Add hadamard

field/src/field.rs

@@ -159,7 +159,7 @@ where
             alpha,
             k: 0,
             kernel,
-            m2l: SvdM2lOperatorData::default(),
+            operator_data: SvdM2lOperatorData::default(),
             transfer_vectors: Vec::new(),
         };
 
@@ -176,7 +176,7 @@ where
         }
 
         (result.transfer_vectors, _) = result.compute_transfer_vectors();
-        result.m2l = result.compute_m2l_operators(order, domain);
+        result.operator_data = result.compute_m2l_operators(order, domain);
 
         result
     }
@@ -194,7 +194,12 @@ where
     type TransferVectorMap = HashMap<usize, usize>;
 
     fn compute_m2l_operators(&self, order: usize, domain: Self::Domain) -> Self::M2LOperators {
-        let mut result = Vec::new();
+        // let mut result = Vec::new();
+
+        let n = 2 * order - 1;
+        let n_p = n + 1;
+        let size_real = n_p * n_p * (n_p / 2 + 1);
+        let mut kernel_data = rlst_mat![c64, (self.transfer_vectors.len() * size_real, 1)];
 
         let mut surface_maps = Vec::new();
         let mut inv_surface_maps = Vec::new();
@@ -204,7 +209,7 @@ where
         // Create a map between corner indices and surface indices
         let (map_corner_to_surface, inv_map_corner_to_surface) = map_corners_to_surface(order);
 
-        for t in self.transfer_vectors.iter() {
+        for (t_idx, t) in self.transfer_vectors.iter().enumerate() {
             let source_equivalent_surface = t.source.compute_surface(&domain, order, self.alpha);
             let nsources = source_equivalent_surface.len() / 3;
 
@@ -412,7 +417,8 @@ where
             );
 
             // Store FFT of kernel for this transfer vector
-            result.push(padded_reflected_kernel_hat);
+            kernel_data.data_mut()[t_idx * size_real..(t_idx + 1) * size_real]
+                .copy_from_slice(padded_reflected_kernel_hat.get_data());
 
             // Save surface and convolution maps
             surface_maps.push(map_surface);
@@ -422,7 +428,7 @@ where
         }
 
         FftM2lOperatorData {
-            m2l: result,
+            kernel_data,
             surface_map: surface_maps,
             inv_surface_map: inv_surface_maps,
             conv_map: conv_maps,
@@ -456,7 +462,7 @@ where
             kernel,
             surf_to_conv_map: HashMap::default(),
             conv_to_surf_map: HashMap::default(),
-            m2l: FftM2lOperatorData::default(),
+            operator_data: FftM2lOperatorData::default(),
             transfer_vectors: Vec::default(),
             transfer_vector_map: HashMap::default(),
         };
@@ -469,7 +475,7 @@ where
         result.conv_to_surf_map = conv_to_surf;
         (result.transfer_vectors, result.transfer_vector_map) = result.compute_transfer_vectors();
 
-        result.m2l = result.compute_m2l_operators(order, domain);
+        result.operator_data = result.compute_m2l_operators(order, domain);
 
         result
     }
@@ -687,6 +693,11 @@ mod test {
         let m2l = fft.compute_m2l_operators(order, domain);
 
         // Test that the number of precomputed kernel interactions matches the number of transfer vectors
-        assert_eq!(m2l.m2l.len(), 16);
+        let n = 2 * order - 1;
+        let n_p = n + 1;
+        let size_real = n_p * n_p * (n_p / 2 + 1);
+        let shp = m2l.kernel_data.shape();
+        let expected = shp.0 / size_real;
+        assert_eq!(expected, 16);
     }
 }

field/src/hadamard.rs

@@ -0,0 +1,92 @@
+//! Fast Hadamard product kernels for coefficient data on octrees.
+
+use num::complex::Complex64;
+
+/// Compute Hadamard product between the Fourier transform of a set of coefficients corresponding to the coefficients of 8
+/// octree siblings, and the 16 unique sets of kernel evaluations corresponding to unique field translation expected when using
+/// translationally invariant, homogenous kernels that have been reflected into the reference cone.
+///
+/// # Arguments
+/// * `order` - The expansion order for the multipole and local expansions.
+/// * `sibling_coefficients` - A set of Fourier transforms of the multipole coefficients arranged on a convolution grid for
+///  a set of siblings, arranged in Morton order.
+/// * `kernel_evaluations` - The set of 16 unique kernel evaluations corresponding to unique transfer vectors in the reference cone for
+/// translationally invariant and homogenous kernels.
+/// * `result` - Mutable vector to store results
+pub fn hadamard_product_sibling(
+    order: usize,
+    sibling_coefficients: &[Complex64],
+    kernel_evaluations: &[Complex64],
+    result: &mut [Complex64],
+) {
+    let n = 2 * order - 1;
+    let &(m, n, o) = &(n, n, n);
+
+    let p = m + 1;
+    let q = n + 1;
+    let r = o + 1;
+    let size_real = p * q * (r / 2 + 1);
+
+    let nsiblings = 8;
+    let nconvolutions = 16;
+
+    for i in 0..nconvolutions {
+        let offset_i = i * size_real;
+
+        // Load each kernel matrix into cache, the most expensive operation
+        let kernel = &kernel_evaluations[offset_i..offset_i + size_real];
+
+        // Iterate over siblings in inner loop, applying same convolution to each one.
+        for j in 0..nsiblings {
+            let offset_j = j * size_real;
+            let signal = &sibling_coefficients[offset_j..offset_j + size_real];
+
+            for k in 0..size_real {
+                result[j * size_real * nconvolutions + i * size_real + k] += kernel[k] * signal[k];
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod test {
+
+    use super::*;
+    use cauchy::c64;
+    use rlst::dense::{rlst_mat, rlst_rand_mat, RawAccess, RawAccessMut};
+
+    #[test]
+    fn test_hadamard_product_sibling() {
+        let nconv = 16;
+        let nsiblings = 8;
+        let order = 5;
+        let n = 2 * order - 1;
+        let &(m, n, o) = &(n, n, n);
+
+        let p = m + 1;
+        let q = n + 1;
+        let r = o + 1;
+        let size_real = p * q * (r / 2 + 1);
+
+        let sibling_coefficients = rlst_rand_mat![c64, (nsiblings * size_real, 1)];
+        let kernel_evaluations = rlst_rand_mat![c64, (nconv * size_real, 1)];
+        let mut result = rlst_mat![c64, (size_real * nsiblings * nconv, 1)];
+
+        hadamard_product_sibling(
+            order,
+            sibling_coefficients.data(),
+            kernel_evaluations.data(),
+            result.data_mut(),
+        );
+
+        // Test just the convolutions being applied to the first child. Also provides an example of how to access convolutions.
+        for i in 0..16 {
+            for j in 0..size_real {
+                let expected =
+                    sibling_coefficients.data()[j] * kernel_evaluations.data()[i * size_real + j];
+                let res = result.data()[j + i * size_real];
+                assert_eq!(res, expected)
+            }
+        }
+    }
+}

field/src/lib.rs

@@ -4,6 +4,7 @@
 pub mod array;
 pub mod fft;
 pub mod field;
+pub mod hadamard;
 pub mod surface;
 pub mod transfer_vector;
 pub mod types;

field/src/types.rs

@@ -1,8 +1,8 @@
 //! Types for storing field translation data.
 use std::collections::HashMap;
 
-use bempp_tools::Array3D;
-use num::Complex;
+use cauchy::c64;
+
 use rlst::{
     common::traits::{Eval, NewLikeSelf},
     dense::{
@@ -18,7 +18,9 @@ pub type SvdM2lEntry =
     Matrix<f64, BaseMatrix<f64, VectorContainer<f64>, Dynamic, Dynamic>, Dynamic, Dynamic>;
 
 /// Simple alias for an Array3D<Complexf64>
-pub type FftM2lEntry = Array3D<Complex<f64>>;
+// pub type FftM2lEntry = Array3D<Complex<f64>>;
+pub type FftM2lEntry =
+    Matrix<c64, BaseMatrix<c64, VectorContainer<c64>, Dynamic, Dynamic>, Dynamic, Dynamic>;
 
 /// A type to store the M2L field translation meta-data and data for an FFT based sparsification in the kernel independent FMM.
 pub struct FftFieldTranslationKiFmm<T>
@@ -35,7 +37,7 @@ where
     pub conv_to_surf_map: HashMap<usize, usize>,
 
     /// Precomputed data required for FFT compressed M2L interaction.
-    pub m2l: FftM2lOperatorData,
+    pub operator_data: FftM2lOperatorData,
 
     /// Unique transfer vectors to lookup m2l unique kernel interactions
     pub transfer_vectors: Vec<TransferVector>,
@@ -59,7 +61,7 @@ where
     pub k: usize,
 
     /// Precomputed data required for SVD compressed M2L interaction.
-    pub m2l: SvdM2lOperatorData,
+    pub operator_data: SvdM2lOperatorData,
 
     /// Unique transfer vectors to lookup m2l unique kernel interactions.
     pub transfer_vectors: Vec<TransferVector>,
@@ -85,10 +87,9 @@ pub struct TransferVector {
 }
 
 /// Container to store precomputed data required for FFT field translations.
-#[derive(Default)]
 pub struct FftM2lOperatorData {
     /// The FFT of unique kernel evaluations for each transfer vector
-    pub m2l: Vec<FftM2lEntry>,
+    pub kernel_data: FftM2lEntry,
 
     /// Map between surface grid indices on unreflected/reflected surfaces.
     /// Each index in the Vec corresponds to an un-reflected transfer vector.
@@ -134,3 +135,17 @@ impl Default for SvdM2lOperatorData {
         }
     }
 }
+
+impl Default for FftM2lOperatorData {
+    fn default() -> Self {
+        let tmp = rlst_mat![c64, (1, 1)];
+
+        FftM2lOperatorData {
+            kernel_data: tmp.new_like_self().eval(),
+            surface_map: Vec::default(),
+            inv_surface_map: Vec::default(),
+            conv_map: Vec::default(),
+            inv_conv_map: Vec::default(),
+        }
+    }
+}