From 30c563c19d95cc56b096f7deb2804fb15d70ea3f Mon Sep 17 00:00:00 2001
From: Magnus Ulimoen <magnus@ulimoen.dev>
Date: Thu, 28 Jan 2021 21:27:25 +0100
Subject: [PATCH] working d1 SBP4

---
 sbp/src/lib.rs                    |  1 +
 sbp/src/operators/algos.rs        | 12 ++++++------
 sbp/src/operators/traditional4.rs | 30 +++++++++++++++++++++++++++++-
 3 files changed, 36 insertions(+), 7 deletions(-)

diff --git a/sbp/src/lib.rs b/sbp/src/lib.rs
index 5b3eaa2..1d09854 100644
--- a/sbp/src/lib.rs
+++ b/sbp/src/lib.rs
@@ -1,5 +1,6 @@
 #![feature(core_intrinsics)]
 #![feature(array_windows)]
+#![feature(array_chunks)]
 
 /// 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 93657e3..3fd83bf 100644
--- a/sbp/src/operators/algos.rs
+++ b/sbp/src/operators/algos.rs
@@ -114,7 +114,7 @@ pub(crate) mod constmatrix {
             let mut v = Self::default();
             for i in 0..M {
                 for j in 0..N {
-                    v[(i, j)] = self[(N - 1 - i, M - 1 - j)].clone()
+                    v[(i, j)] = self[(M - 1 - i, N - 1 - j)].clone()
                 }
             }
             v
@@ -276,12 +276,11 @@ pub(crate) fn diff_op_1d_slice_matrix<const M: usize, const N: usize, const D: u
     for (window, f) in prev
         .array_windows::<D>()
         .skip(window_elems_to_skip)
-        .zip(fut.iter_mut().skip(M))
+        .zip(fut.array_chunks_mut::<1>().skip(M))
         .take(nx - 2 * M)
     {
-        let fut = ColVector::<Float, 1>::map_to_col_mut(unsafe {
-            std::mem::transmute::<&mut Float, &mut [Float; 1]>(f)
-        });
+        // impl From here?
+        let fut = ColVector::<Float, 1>::map_to_col_mut(f);
         let prev = ColVector::<_, D>::map_to_col(window);
 
         diag.matmul_into(prev, fut);
@@ -297,7 +296,8 @@ pub(crate) fn diff_op_1d_slice_matrix<const M: usize, const N: usize, const D: u
     };
 
     {
-        let prev = ColVector::<_, N>::map_to_col((&prev[nx - N..]).try_into().unwrap());
+        let prev = prev.array_windows::<N>().last().unwrap();
+        let prev = ColVector::<_, N>::map_to_col(prev);
         let fut = ColVector::<_, M>::map_to_col_mut((&mut fut[nx - M..]).try_into().unwrap());
 
         flipped.matmul_into(prev, fut);
diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs
index 649c43c..c29944d 100644
--- a/sbp/src/operators/traditional4.rs
+++ b/sbp/src/operators/traditional4.rs
@@ -10,10 +10,14 @@ impl SBP4 {
     const HBLOCK: &'static [Float] = &[
         17.0 / 48.0, 59.0 / 48.0, 43.0 / 48.0, 49.0 / 48.0,
     ];
+
     #[rustfmt::skip]
     const DIAG: &'static [Float] = &[
         1.0 / 12.0, -2.0 / 3.0, 0.0, 2.0 / 3.0, -1.0 / 12.0,
     ];
+
+    const DIAG_MATRIX: super::RowVector<Float, 5> =
+        super::RowVector::new([[1.0 / 12.0, -2.0 / 3.0, 0.0, 2.0 / 3.0, -1.0 / 12.0]]);
     #[rustfmt::skip]
     const BLOCK: &'static [&'static [Float]] = &[
         &[-24.0/17.0, 59.0/34.0, -4.0/17.0, -3.0/34.0],
@@ -21,6 +25,13 @@ impl SBP4 {
         &[4.0/43.0, -59.0/86.0, 0.0, 59.0/86.0, -4.0/43.0],
         &[3.0/98.0, 0.0, -59.0/98.0, 0.0, 32.0/49.0, -4.0/49.0]
     ];
+    #[rustfmt::skip]
+    const BLOCK_MATRIX: super::Matrix<Float, 4, 6> = super::Matrix::new([
+        [-24.0/17.0, 59.0/34.0, -4.0/17.0, -3.0/34.0, 0.0, 0.0],
+        [-1.0/2.0, 0.0, 1.0/2.0, 0.0, 0.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]
+    ]);
 
     #[rustfmt::skip]
     const D2DIAG: &'static [Float] = &[
@@ -71,6 +82,22 @@ impl SbpOperator1d for SBP4 {
     }
 }
 
+fn diff_op_row_local(prev: ndarray::ArrayView2<Float>, mut fut: ndarray::ArrayViewMut2<Float>) {
+    for (p, mut f) in prev
+        .axis_iter(ndarray::Axis(0))
+        .zip(fut.axis_iter_mut(ndarray::Axis(0)))
+    {
+        super::diff_op_1d_slice_matrix(
+            &SBP4::BLOCK_MATRIX,
+            &SBP4::DIAG_MATRIX,
+            super::Symmetry::AntiSymmetric,
+            super::OperatorType::Normal,
+            p.as_slice().unwrap(),
+            f.as_slice_mut().unwrap(),
+        )
+    }
+}
+
 impl SbpOperator2d for SBP4 {
     fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
         assert_eq!(prev.shape(), fut.shape());
@@ -81,7 +108,8 @@ impl SbpOperator2d for SBP4 {
 
         match (prev.strides(), fut.strides()) {
             ([_, 1], [_, 1]) => {
-                diff_op_row(SBP4::BLOCK, SBP4::DIAG, symmetry, optype)(prev, fut);
+                //diff_op_row(SBP4::BLOCK, SBP4::DIAG, symmetry, optype)(prev, fut);
+                diff_op_row_local(prev, fut)
             }
             ([1, _], [1, _]) => {
                 diff_op_col(SBP4::BLOCK, SBP4::DIAG, symmetry, optype)(prev, fut);