Add individual low rank appx

field/src/field.rs

@@ -22,6 +22,7 @@ use bempp_tree::{
     implementations::helpers::find_corners, types::domain::Domain, types::morton::MortonKey,
 };
 
+use crate::types::{SvdFieldTranslationKiFmmIA, SvdM2lOperatorDataIA};
 use crate::{
     array::flip3,
     fft::Fft,
@@ -164,6 +165,116 @@ where
     }
 }
 
+impl<T, U> FieldTranslationData<U> for SvdFieldTranslationKiFmmIA<T, U>
+where
+    T: Float + Default,
+    T: Scalar<Real = T> + Gemm,
+    U: Kernel<T = T> + Default,
+    Array<T, BaseArray<T, VectorContainer<T>, 2>, 2>: MatrixSvd<Item = T>,
+{
+    type TransferVector = Vec<TransferVector>;
+    type M2LOperators = SvdM2lOperatorDataIA<T>;
+    type Domain = Domain<T>;
+
+    fn ncoeffs(&self, order: usize) -> usize {
+        6 * (order - 1).pow(2) + 2
+    }
+
+    fn compute_m2l_operators<'a>(&self, order: usize, domain: Self::Domain) -> Self::M2LOperators {
+        // Compute unique M2L interactions at Level 3 (smallest choice with all vectors)
+
+        // Compute interaction matrices between source and unique targets, defined by unique transfer vectors
+        let nrows = self.ncoeffs(order);
+        let ncols = self.ncoeffs(order);
+
+        let mut u = Vec::new();
+        let mut vt = Vec::new();
+
+        for (_i, t) in self.transfer_vectors.iter().enumerate() {
+            let source_equivalent_surface = t.source.compute_surface(&domain, order, self.alpha);
+            let target_check_surface = t.target.compute_surface(&domain, order, self.alpha);
+
+            let mut tmp_gram_t = rlst_dynamic_array2!(T, [nrows, ncols]);
+
+            self.kernel.assemble_st(
+                EvalType::Value,
+                &target_check_surface[..],
+                &source_equivalent_surface[..],
+                tmp_gram_t.data_mut(),
+            );
+
+            let mut u_i = rlst_dynamic_array2!(T, [nrows, self.k]);
+            let mut sigma_i = vec![T::zero(); self.k];
+            let mut vt_i = rlst_dynamic_array2!(T, [self.k, ncols]);
+
+            tmp_gram_t
+                .into_svd_alloc(
+                    u_i.view_mut(),
+                    vt_i.view_mut(),
+                    &mut sigma_i,
+                    SvdMode::Full,
+                )
+                .unwrap();
+
+            // Retain such that 95% of energy of singular values is retained.
+            let rank = retain_energy(&sigma_i, self.threshold);
+
+            // let rank = 30;
+            // println!("{:?} rank {:?}", _i, rank);
+
+            let mut u_i_compressed = rlst_dynamic_array2!(T, [nrows, rank]);
+            let mut vt_i_compressed_= rlst_dynamic_array2!(T, [rank, ncols]);
+
+            let mut sigma_mat_i_compressed = rlst_dynamic_array2!(T, [rank, rank]);
+
+            u_i_compressed.fill_from(u_i.into_subview([0, 0], [nrows, rank]));
+            vt_i_compressed_.fill_from(vt_i.into_subview([0, 0], [rank, ncols]));
+
+            for (j, s) in sigma_i.iter().enumerate().take(rank) {
+                unsafe {
+                    *sigma_mat_i_compressed.get_unchecked_mut([j, j]) = T::from(*s).unwrap();
+                }
+            }
+
+            let vt_i_compressed = empty_array::<T, 2>()
+            .simple_mult_into_resize(
+                sigma_mat_i_compressed.view(), vt_i_compressed_.view()
+            );
+
+            // Store compressed M2L oeprators
+            u.push(u_i_compressed);
+            vt.push(vt_i_compressed);
+        }
+
+        SvdM2lOperatorDataIA { u, vt }
+    }
+}
+
+fn retain_energy<T: Float + Default + Scalar<Real = T> + Gemm>(
+    singular_values: &Vec<T>,
+    percentage: T,
+) -> usize {
+    // Calculate the total energy.
+    let total_energy: T = singular_values.iter().map(|&s| s * s).sum();
+
+    // Calculate the threshold energy to retain.
+    let threshold_energy = total_energy * (percentage / T::one());
+
+    // Iterate over singular values to find the minimum set that retains the desired energy.
+    let mut cumulative_energy = T::zero();
+    let mut significant_values = Vec::new();
+
+    for (i, &value) in singular_values.iter().enumerate() {
+        cumulative_energy += value * value;
+        significant_values.push(value);
+        if cumulative_energy >= threshold_energy {
+            return i + 1;
+        }
+    }
+
+    significant_values.len()
+}
+
 impl<T, U> SvdFieldTranslationKiFmm<T, U>
 where
     T: Float + Default,
@@ -205,6 +316,40 @@ where
     }
 }
 
+impl<T, U> SvdFieldTranslationKiFmmIA<T, U>
+where
+    T: Float + Default,
+    T: Scalar<Real = T> + rlst_blis::interface::gemm::Gemm,
+    U: Kernel<T = T> + Default,
+    Array<T, BaseArray<T, VectorContainer<T>, 2>, 2>: MatrixSvd<Item = T>,
+{
+    /// Constructor for SVD field translation struct for the kernel independent FMM (KiFMM).
+    ///
+    /// # Arguments
+    /// * `kernel` - The kernel being used, only compatible with homogenous, translationally invariant kernels.
+    /// * `k` - The maximum rank to be used in SVD compression for the translation operators, if none is specified will be taken as  max({50, max_column_rank})
+    /// * `order` - The expansion order for the multipole and local expansions.
+    /// * `domain` - Domain associated with the global point set.
+    /// * `alpha` - The multiplier being used to modify the diameter of the surface grid uniformly along each coordinate axis.
+    pub fn new(kernel: U, threshold: T, order: usize, domain: Domain<T>, alpha: T) -> Self {
+        let mut result = SvdFieldTranslationKiFmmIA {
+            alpha,
+            k: 0,
+            threshold,
+            kernel,
+            operator_data: SvdM2lOperatorDataIA::default(),
+            transfer_vectors: vec![],
+        };
+
+        let ncoeffs = result.ncoeffs(order);
+        result.k = ncoeffs;
+        result.transfer_vectors = compute_transfer_vectors();
+        result.operator_data = result.compute_m2l_operators(order, domain);
+
+        result
+    }
+}
+
 impl<T, U> FieldTranslationData<U> for FftFieldTranslationKiFmm<T, U>
 where
     T: Scalar<Real = T> + Float + Default + Fft,

field/src/types.rs

@@ -63,6 +63,31 @@ where
     pub kernel: U,
 }
 
+/// A type to store the M2L field translation meta-data  and datafor an SVD based sparsification in the kernel independent FMM.
+pub struct SvdFieldTranslationKiFmmIA<T, U>
+where
+    T: Scalar<Real = T> + Float + Default + rlst_blis::interface::gemm::Gemm,
+    U: Kernel<T = T> + Default,
+{
+    /// Amount to dilate inner check surface by when computing operator.
+    pub alpha: T,
+
+    /// Maximum rank taken for SVD compression
+    pub k: usize,
+
+    /// Amount of energy of each M2L operator retained in SVD compression
+    pub threshold: T,
+
+    /// Precomputed data required for SVD compressed M2L interaction.
+    pub operator_data: SvdM2lOperatorDataIA<T>,
+
+    /// Unique transfer vectors to lookup m2l unique kernel interactions.
+    pub transfer_vectors: Vec<TransferVector>,
+
+    /// The associated kernel with this translation operator.
+    pub kernel: U,
+}
+
 /// A type to store a transfer vector between a `source` and `target` Morton key.
 #[derive(Debug)]
 pub struct TransferVector {
@@ -119,3 +144,17 @@ where
         SvdM2lOperatorData { u, st_block, c }
     }
 }
+
+/// Container to store precomputed data required for SVD field translations.
+/// See Fong & Darve (2009) for the definitions of 'fat' and 'thin' M2L matrices.
+#[derive(Default)]
+pub struct SvdM2lOperatorDataIA<T>
+where
+    T: Scalar,
+{
+    /// Left singular vectors from SVD of fat M2L matrix.
+    pub u: Vec<SvdM2lEntry<T>>,
+
+    /// Right singular vectors from SVD of thin M2L matrix, cutoff to a maximum rank of 'k'.
+    pub vt: Vec<SvdM2lEntry<T>>,
+}

fmm/examples/single_node.rs

@@ -4,7 +4,7 @@ use itertools::Itertools;
 
 use rlst_dense::traits::RawAccess;
 
-use bempp_field::types::FftFieldTranslationKiFmm;
+use bempp_field::types::{FftFieldTranslationKiFmm, SvdFieldTranslationKiFmm, SvdFieldTranslationKiFmmIA};
 use bempp_fmm::{
     charge::build_charge_dict,
     types::{FmmDataUniform, KiFmmLinear},
@@ -32,6 +32,7 @@ fn main() {
     let kernel = Laplace3dKernel::default();
     let m2l_data: FftFieldTranslationKiFmm<f32, Laplace3dKernel<f32>> =
         FftFieldTranslationKiFmm::new(kernel.clone(), order, *tree.get_domain(), alpha_inner);
+
     // let m2l_data = SvdFieldTranslationKiFmm::new(
     //     kernel.clone(),
     //     Some(80),
@@ -40,6 +41,14 @@ fn main() {
     //     alpha_inner,
     // );
 
+    // let m2l_data = SvdFieldTranslationKiFmmIA::new(
+    //     kernel.clone(),
+    //     0.9999,
+    //     order,
+    //     *tree.get_domain(),
+    //     alpha_inner
+    // );
+
     let fmm = KiFmmLinear::new(order, alpha_inner, alpha_outer, kernel, tree, m2l_data);
 
     // Form charge dict, matching charges with their associated global indices

fmm/examples/single_node_matrix.rs

@@ -5,7 +5,7 @@ use itertools::Itertools;
 
 use rlst_dense::traits::RawAccess;
 
-use bempp_field::types::SvdFieldTranslationKiFmm;
+use bempp_field::types::{SvdFieldTranslationKiFmm, SvdFieldTranslationKiFmmIA};
 use bempp_fmm::charge::build_charge_dict;
 use bempp_kernel::laplace_3d::Laplace3dKernel;
 use bempp_traits::{fmm::FmmLoop, tree::Tree};
@@ -23,7 +23,7 @@ fn main() {
 
     // Test matrix input
     let points = points_fixture::<f32>(npoints, None, None);
-    let ncharge_vecs = 3;
+    let ncharge_vecs = 5;
     let depth = 4;
 
     let mut charge_mat = vec![vec![0.0; npoints]; ncharge_vecs];
@@ -39,12 +39,19 @@ fn main() {
     let tree = SingleNodeTree::new(points.data(), false, None, Some(depth), &global_idxs, true);
 
     // Precompute the M2L data
-    let m2l_data = SvdFieldTranslationKiFmm::new(
+    // let m2l_data = SvdFieldTranslationKiFmm::new(
+    //     kernel.clone(),
+    //     Some(80),
+    //     order,
+    //     *tree.get_domain(),
+    //     alpha_inner,
+    // );
+    let m2l_data = SvdFieldTranslationKiFmmIA::new(
         kernel.clone(),
-        Some(80),
+        0.9999,
         order,
         *tree.get_domain(),
-        alpha_inner,
+        alpha_inner
     );
 
     let fmm = bempp_fmm::types::KiFmmLinearMatrix::new(