diff --git a/sbp/src/lib.rs b/sbp/src/lib.rs index 62cb872..52a1afa 100644 --- a/sbp/src/lib.rs +++ b/sbp/src/lib.rs @@ -1,6 +1,7 @@ #![cfg_attr(feature = "fast-float", feature(core_intrinsics))] #![feature(array_windows)] #![feature(array_chunks)] +#![feature(const_fn_floating_point_arithmetic)] /// Type used for floats, configure with the `f32` feature #[cfg(feature = "f32")] diff --git a/sbp/src/operators/algos.rs b/sbp/src/operators/algos.rs index 11b9923..7fe8dd6 100644 --- a/sbp/src/operators/algos.rs +++ b/sbp/src/operators/algos.rs @@ -6,7 +6,7 @@ pub(crate) mod constmatrix { #[derive(Debug, Clone, PartialEq, Eq, Hash)] #[repr(C)] pub struct Matrix { - data: [[T; N]; M], + pub data: [[T; N]; M], } pub type RowVector = Matrix; pub type ColVector = Matrix; @@ -75,19 +75,6 @@ pub(crate) mod constmatrix { ) -> impl ExactSizeIterator + DoubleEndedIterator { 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 ColVector { @@ -244,9 +231,77 @@ pub(crate) mod constmatrix { } } } + + pub(crate) const fn flip_ud( + mut m: Matrix, + ) -> Matrix { + 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) const fn flip_lr( + mut m: Matrix, + ) -> Matrix { + 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( + mut m: Matrix, + ) -> Matrix { + 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::{ColVector, Matrix, RowVector}; +pub(crate) use constmatrix::{flip_lr, flip_sign, flip_ud, ColVector, Matrix, RowVector}; #[inline(always)] pub(crate) fn diff_op_1d_matrix( diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs index f480cda..417b799 100644 --- a/sbp/src/operators/traditional4.rs +++ b/sbp/src/operators/traditional4.rs @@ -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], [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 = super::Matrix::new([ - [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], - ]); + const BLOCKEND_MATRIX: super::Matrix = + super::flip_sign(super::flip_ud(super::flip_lr(Self::BLOCK_MATRIX))); #[rustfmt::skip] const D2DIAG: &'static [Float] = &[ @@ -234,11 +229,3 @@ fn test_trad4() { 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); -} diff --git a/sbp/src/operators/traditional8.rs b/sbp/src/operators/traditional8.rs index 48891f8..32b392c 100644 --- a/sbp/src/operators/traditional8.rs +++ b/sbp/src/operators/traditional8.rs @@ -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], [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: Matrix = Matrix::new([ - [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], - ]); + const BLOCKEND_MATRIX: super::Matrix = + super::flip_sign(super::flip_ud(super::flip_lr(Self::BLOCK_MATRIX))); } impl SbpOperator1d for SBP8 { @@ -215,11 +206,3 @@ fn test_trad8() { 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); -} diff --git a/sbp/src/operators/upwind4.rs b/sbp/src/operators/upwind4.rs index 66b9852..d2a4eb4 100644 --- a/sbp/src/operators/upwind4.rs +++ b/sbp/src/operators/upwind4.rs @@ -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], [ 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 = Matrix::new([ - [ -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], - ]); + const BLOCKEND_MATRIX: Matrix = + super::flip_sign(super::flip_ud(super::flip_lr(Self::BLOCK_MATRIX))); #[rustfmt::skip] const DISS_BLOCK: &'static [&'static [Float]] = &[ @@ -495,11 +490,3 @@ fn upwind4_test2() { 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); -}