mulimoen
/
SummationByParts
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

make flip const functions

This commit is contained in:
Magnus Ulimoen 2021-02-01 18:37:59 +01:00
parent 1f3aa2c116
commit f7f8a7ffff
5 changed files with 77 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
sbp/src/lib.rs
View File

@ -1,6 +1,7 @@
#![cfg_attr(feature = "fast-float", feature(core_intrinsics))] #![cfg_attr(feature = "fast-float", feature(core_intrinsics))]
#![feature(array_windows)] #![feature(array_windows)]
#![feature(array_chunks)] #![feature(array_chunks)]
#![feature(const_fn_floating_point_arithmetic)]
/// Type used for floats, configure with the `f32` feature /// Type used for floats, configure with the `f32` feature
#[cfg(feature = "f32")] #[cfg(feature = "f32")]

85
sbp/src/operators/algos.rs
View File

@ -6,7 +6,7 @@ pub(crate) mod constmatrix {
#[derive(Debug, Clone, PartialEq, Eq, Hash)] #[derive(Debug, Clone, PartialEq, Eq, Hash)]
#[repr(C)] #[repr(C)]
pub struct Matrix<T, const M: usize, const N: usize> { pub struct Matrix<T, const M: usize, const N: usize> {
data: [[T; N]; M], pub data: [[T; N]; M],
} }
pub type RowVector<T, const N: usize> = Matrix<T, 1, N>; pub type RowVector<T, const N: usize> = Matrix<T, 1, N>;
pub type ColVector<T, const N: usize> = Matrix<T, N, 1>; pub type ColVector<T, const N: usize> = Matrix<T, N, 1>;
@ -75,19 +75,6 @@ pub(crate) mod constmatrix {
) -> impl ExactSizeIterator<Item = &[T; N]> + DoubleEndedIterator<Item = &[T; N]> { ) -> impl ExactSizeIterator<Item = &[T; N]> + DoubleEndedIterator<Item = &[T; N]> {
self.data.iter() self.data.iter()
} }
pub fn flip(&self) -> Self
where
T: Default + Copy,
{
let mut v = Self::default();
for i in 0..M {
for j in 0..N {
v[(i, j)] = self[(M - 1 - i, N - 1 - j)]
}
}
v
}
} }
impl<T, const N: usize> ColVector<T, N> { impl<T, const N: usize> ColVector<T, N> {
@ -244,9 +231,77 @@ pub(crate) mod constmatrix {
} }
} }
} }
pub(crate) const fn flip_ud<const M: usize, const N: usize>(
mut m: Matrix<super::Float, M, N>,
) -> Matrix<super::Float, M, N> {
let mut i = 0;
while i < M / 2 {
let tmp = m.data[i];
m.data[i] = m.data[M - 1 - i];
m.data[M - 1 - i] = tmp;
i += 1;
}
m
} }
pub(crate) use constmatrix::{ColVector, Matrix, RowVector}; pub(crate) const fn flip_lr<const M: usize, const N: usize>(
mut m: Matrix<super::Float, M, N>,
) -> Matrix<super::Float, M, N> {
let mut i = 0;
while i < M {
let mut j = 0;
while j < N / 2 {
let tmp = m.data[i][j];
m.data[i][j] = m.data[i][N - 1 - j];
m.data[i][N - 1 - j] = tmp;
j += 1;
}
i += 1;
}
m
}
/// Flip all sign bits
pub(crate) const fn flip_sign<const M: usize, const N: usize>(
mut m: Matrix<super::Float, M, N>,
) -> Matrix<super::Float, M, N> {
let mut i = 0;
while i < M {
let mut j = 0;
while j < N {
m.data[i][j] = -m.data[i][j];
j += 1;
}
i += 1;
}
m
}
mod flipping {
use super::*;
#[test]
fn flip_lr_test() {
let m = Matrix::new([[1.0, 2.0, 3.0, 4.0]]);
let flipped = flip_lr(m);
assert_eq!(flipped, Matrix::new([[4.0, 3.0, 2.0, 1.0]]));
let m = Matrix::new([[1.0, 2.0, 3.0, 4.0, 5.0]]);
let flipped = flip_lr(m);
assert_eq!(flipped, Matrix::new([[5.0, 4.0, 3.0, 2.0, 1.0]]));
}
#[test]
fn flip_ud_test() {
let m = Matrix::new([[1.0], [2.0], [3.0], [4.0]]);
let flipped = flip_ud(m);
assert_eq!(flipped, Matrix::new([[4.0], [3.0], [2.0], [1.0]]));
let m = Matrix::new([[1.0], [2.0], [3.0], [4.0], [5.0]]);
let flipped = flip_ud(m);
assert_eq!(flipped, Matrix::new([[5.0], [4.0], [3.0], [2.0], [1.0]]));
}
}
}
pub(crate) use constmatrix::{flip_lr, flip_sign, flip_ud, ColVector, Matrix, RowVector};
#[inline(always)] #[inline(always)]
pub(crate) fn diff_op_1d_matrix<const M: usize, const N: usize, const D: usize>( pub(crate) fn diff_op_1d_matrix<const M: usize, const N: usize, const D: usize>(

17
sbp/src/operators/traditional4.rs
View File

@ -32,13 +32,8 @@ impl SBP4 {
[4.0/43.0, -59.0/86.0, 0.0, 59.0/86.0, -4.0/43.0, 0.0], [4.0/43.0, -59.0/86.0, 0.0, 59.0/86.0, -4.0/43.0, 0.0],
[3.0/98.0, 0.0, -59.0/98.0, 0.0, 32.0/49.0, -4.0/49.0] [3.0/98.0, 0.0, -59.0/98.0, 0.0, 32.0/49.0, -4.0/49.0]
]); ]);
#[rustfmt::skip] const BLOCKEND_MATRIX: super::Matrix<Float, 4, 6> =
const BLOCKEND_MATRIX: super::Matrix<Float, 4, 6> = super::Matrix::new([ super::flip_sign(super::flip_ud(super::flip_lr(Self::BLOCK_MATRIX)));
[4.0/49.0, -32.0/49.0, 0.0, 59.0/98.0, 0.0, -3.0/98.0 ],
[0.0, 4.0/43.0, -59.0/86.0, 0.0, 59.0/86.0, -4.0/43.0],
[0.0, 0.0, 0.0, -1.0/2.0, 0.0, 1.0/2.0],
[0.0, 0.0, 3.0/34.0, 4.0/17.0, -59.0/34.0, 24.0/17.0],
]);
#[rustfmt::skip] #[rustfmt::skip]
const D2DIAG: &'static [Float] = &[ const D2DIAG: &'static [Float] = &[
@ -234,11 +229,3 @@ fn test_trad4() {
1e-1, 1e-1,
); );
} }
#[test]
fn block_equality() {
let mut flipped_inverted = SBP4::BLOCK_MATRIX.flip();
flipped_inverted *= -1.0;
approx::assert_ulps_eq!(SBP4::BLOCKEND_MATRIX, flipped_inverted, max_ulps = 1);
}

21
sbp/src/operators/traditional8.rs
View File

@ -40,17 +40,8 @@ impl SBP8 {
[2.75587615266177e-02, -8.71295642560637e-02, 5.01135077563584e-02, 7.68229253600969e-02, 4.60181213406519e-02, -7.55873581663580e-01, 0.00000000000000e+00, 8.21713248844682e-01, -2.16615355227872e-01, 4.12600676624518e-02, -3.86813134335486e-03, 0.0], [2.75587615266177e-02, -8.71295642560637e-02, 5.01135077563584e-02, 7.68229253600969e-02, 4.60181213406519e-02, -7.55873581663580e-01, 0.00000000000000e+00, 8.21713248844682e-01, -2.16615355227872e-01, 4.12600676624518e-02, -3.86813134335486e-03, 0.0],
[5.84767272160451e-03, -1.08336661209337e-02, -1.42810403117803e-02, 3.50919361287023e-02, -1.22244235731112e-02, 1.19411743193552e-01, -7.51668243727123e-01, 0.00000000000000e+00, 7.92601963555477e-01, -1.98150490888869e-01, 3.77429506454989e-02, -3.53840162301552e-03], [5.84767272160451e-03, -1.08336661209337e-02, -1.42810403117803e-02, 3.50919361287023e-02, -1.22244235731112e-02, 1.19411743193552e-01, -7.51668243727123e-01, 0.00000000000000e+00, 7.92601963555477e-01, -1.98150490888869e-01, 3.77429506454989e-02, -3.53840162301552e-03],
]); ]);
#[rustfmt::skip] const BLOCKEND_MATRIX: super::Matrix<Float, 8, 12> =
const BLOCKEND_MATRIX: Matrix<Float, 8, 12> = Matrix::new([ super::flip_sign(super::flip_ud(super::flip_lr(Self::BLOCK_MATRIX)));
[0.0035384016230155199, -0.037742950645498902, 0.19815049088886899, -0.79260196355547696, -0.0, 0.75166824372712304, -0.11941174319355199, 0.012224423573111201, -0.035091936128702303, 0.0142810403117803, 0.010833666120933699, -0.0058476727216045096],
[-0.0, 0.0038681313433548601, -0.041260067662451799, 0.21661535522787201, -0.821713248844682, -0.0, 0.75587358166357999, -0.046018121340651898, -0.076822925360096897, -0.050113507756358401, 0.087129564256063705, -0.027558761526617698],
[-0.0, -0.0, 0.0027934857464329702, -0.029797181295285, 0.094272792663829805, -0.54587651996612696, -0.0, 0.205756876210586, 0.39561116774840099, 0.020799982439143602, -0.22287532317150199, 0.079314719624520302],
[-0.0, -0.0, -0.0, 0.0086536439119011006, -0.029896495621603899, 0.102950081118829, -0.63739245543855805, -0.0, 0.54336288636530095, 0.085109057027750595, -0.090998664615454999, 0.01821194725190089],
[-0.0, -0.0, -0.0, -0.0, 0.019698131073842998, 0.039447042394264102, -0.28128521251909999, -0.124714160903055, -0.0, 0.055770702144577897, 0.41917213003600801, -0.128088632226564],
[-0.0, -0.0, -0.0, -0.0, -0.055988255885229599, 0.179720579323281, -0.10329058280084499, -0.13643348652854601, -0.389516189693211, -0.0, 0.54432974445498405, -0.038821808870425301],
[-0.0, -0.0, -0.0, -0.0, -0.0071669648308011497, -0.052726783847581303, 0.18675940343293501, 0.024615176293723499, -0.49400862680798402, -0.091851192507295606, -0.0, 0.43437898826698501],
[-0.0, -0.0, -0.0, -0.0, 0.020015058331576099, 0.086285816263333495, -0.34386522738887298, -0.0254881486107905, 0.78102816812674902, 0.033893192260150001, -2.2474124634140402, 1.6955436044319001],
]);
} }
impl SbpOperator1d for SBP8 { impl SbpOperator1d for SBP8 {
@ -215,11 +206,3 @@ fn test_trad8() {
1e-4, 1e-4,
); );
} }
#[test]
fn block_equality() {
let mut flipped_inverted = SBP8::BLOCK_MATRIX.flip();
flipped_inverted *= -1.0;
approx::assert_ulps_eq!(SBP8::BLOCKEND_MATRIX, flipped_inverted, max_ulps = 1);
}

17
sbp/src/operators/upwind4.rs
View File

@ -192,13 +192,8 @@ impl Upwind4 {
[ 20.0 / 123.0, -69.0 / 82.0, 0.0, 227.0 / 246.0, -12.0 / 41.0, 2.0 / 41.0, 0.0], [ 20.0 / 123.0, -69.0 / 82.0, 0.0, 227.0 / 246.0, -12.0 / 41.0, 2.0 / 41.0, 0.0],
[ 3.0 / 298.0, 16.0 / 149.0, -227.0 / 298.0, 0.0, 126.0 / 149.0, -36.0 / 149.0, 6.0 / 149.0], [ 3.0 / 298.0, 16.0 / 149.0, -227.0 / 298.0, 0.0, 126.0 / 149.0, -36.0 / 149.0, 6.0 / 149.0],
]); ]);
#[rustfmt::skip] const BLOCKEND_MATRIX: Matrix<Float, 4, 7> =
const BLOCKEND_MATRIX: Matrix<Float, 4, 7> = Matrix::new([ super::flip_sign(super::flip_ud(super::flip_lr(Self::BLOCK_MATRIX)));
[ -6.0 / 149.0, 36.0 / 149.0, -126.0 / 149.0, 0.0, 227.0 / 298.0, -16.0 / 149.0, -3.0 / 298.0],
[ 0.0, -2.0/41.0, 12.0/41.0, -227.0/246.0, 0.0, 69.0/82.0, -20.0/123.0],
[ 0.0, 0.0, -2.0/61.0, 16.0/183.0, -69.0/122.0, 0.0, 187.0/366.0],
[ 0.0, 0.0, 0.0, 3.0 / 98.0, 20.0 / 49.0, -187.0 / 98.0, 72.0/49.0],
]);
#[rustfmt::skip] #[rustfmt::skip]
const DISS_BLOCK: &'static [&'static [Float]] = &[ const DISS_BLOCK: &'static [&'static [Float]] = &[
@ -495,11 +490,3 @@ fn upwind4_test2() {
1e-1, 1e-1,
); );
} }
#[test]
fn block_equality() {
let mut flipped_inverted = Upwind4::BLOCK_MATRIX.flip();
flipped_inverted *= -1.0;
approx::assert_ulps_eq!(Upwind4::BLOCKEND_MATRIX, flipped_inverted, max_ulps = 1);
}