Remove an expensive heap allocation (#150)

Co-authored-by: Erik Tollerud <[email protected]>

Cargo.toml

@@ -23,3 +23,6 @@ float_eq = "1.0.0"
 [lints.clippy]
 needless_return = "allow"
 implicit_return = "warn"
+
+[profile.release]
+codegen-units = 1

benchmarks/benchmark.py

@@ -18,7 +18,7 @@
 grid_spacing = [1, 2, 3]
 grid = VectorGrid(field, grid_spacing)
 
-seedlist = 2 ** np.arange(12)
+seedlist = 2 ** np.arange(17)
 times = []
 for nseeds in seedlist:
     dts = []

src/interp.rs

@@ -1,6 +1,6 @@
 //! Helper functions for interpolation.
 
-use numpy::ndarray::{Array, Array1, ArrayBase, Data, Ix3};
+use numpy::ndarray::{Array1, ArrayBase, Data, Ix3};
 
 /// Trilinear-interpolation of a scalar defined on
 /// the eight corners of a cuboid.
@@ -18,20 +18,21 @@ where
     }
     let m_x = 1. - x;
 
-    let mut c = Array::zeros((2, 2));
+    let mut c: [f64; 4] = [0.0; 4];
     // Interpolate over x
     for iy in 0..2 {
         for iz in 0..2 {
-            c[[iy, iz]] = values[[0, iy, iz]] * m_x[[0]] + values[[1, iy, iz]] * x[[0]];
+            let iix = (2 * iy) + iz;
+            c[iix] = values[[0, iy, iz]] * m_x[[0]] + values[[1, iy, iz]] * x[[0]];
         }
     }
 
     // Interpolate over y
-    let mut c1 = Array::zeros(2);
+    let mut c1: [f64; 2] = [0.0; 2];
     for iz in 0..2 {
-        c1[[iz]] = c[[0, iz]] * m_x[[1]] + c[[1, iz]] * x[[1]];
+        c1[iz] = c[iz] * m_x[[1]] + c[iz + 2] * x[[1]];
     }
 
     // Interpolate over z
-    return c1[[0]] * m_x[[2]] + c1[[1]] * x[[2]];
+    return c1[0] * m_x[[2]] + c1[1] * x[[2]];
 }