diff --git a/fmm/src/field_translation/source.rs b/fmm/src/field_translation/source.rs index b8eb5004..adf9f171 100644 --- a/fmm/src/field_translation/source.rs +++ b/fmm/src/field_translation/source.rs @@ -17,7 +17,7 @@ use bempp_tree::types::{morton::MortonKey, single_node::SingleNodeTree}; use crate::{ constants::{M2M_MAX_CHUNK_SIZE, P2M_MAX_CHUNK_SIZE}, helpers::find_chunk_size, - types::{FmmDataAdaptive, FmmDataUniform, KiFmmLinear, FmmDataUniformMatrix}, + types::{FmmDataAdaptive, FmmDataUniform, KiFmmLinear, FmmDataUniformMatrix, KiFmmLinearMatrix}, }; use rlst::{ common::traits::*, @@ -292,7 +292,7 @@ where } -impl SourceTranslation for FmmDataUniformMatrix, T, U, V>, V> +impl SourceTranslation for FmmDataUniformMatrix, T, U, V>, V> where T: Kernel + ScaleInvariantKernel + std::marker::Send + std::marker::Sync, U: FieldTranslationData + std::marker::Sync + std::marker::Send, @@ -381,39 +381,56 @@ where /// Multipole to multipole translations, multithreaded over all boxes at a given level. fn m2m<'a>(&self, level: u64) { if let Some(child_sources) = self.fmm.tree().get_keys(level) { - // let ncoeffs = self.fmm.m2l.ncoeffs(self.fmm.order); - // let nsiblings = 8; - - // // 1. Lookup parents and corresponding children that exist for this set of sources - // // Must explicitly lookup as boxes may be empty at this level, and the next. - // let parent_targets: HashSet = - // child_sources.iter().map(|source| source.parent()).collect(); - // let mut parent_targets = parent_targets.into_iter().collect_vec(); - // parent_targets.sort(); - // let nparents = parent_targets.len(); - // let mut parent_multipoles = Vec::new(); - // for parent in parent_targets.iter() { - // let parent_index_pointer = *self.level_index_pointer[(level - 1) as usize] - // .get(parent) - // .unwrap(); - // let parent_multipole = - // self.level_multipoles[(level - 1) as usize][parent_index_pointer]; - // parent_multipoles.push(parent_multipole); - // } - - // let n_child_sources = child_sources.len(); - // let min: &MortonKey = &child_sources[0]; - // let max = &child_sources[n_child_sources - 1]; - // let min_idx = self.fmm.tree().key_to_index.get(min).unwrap(); - // let max_idx = self.fmm.tree().key_to_index.get(max).unwrap(); - - // let child_multipoles = &self.multipoles[min_idx * ncoeffs..(max_idx + 1) * ncoeffs]; - - // let mut max_chunk_size = nparents; - // if max_chunk_size > M2M_MAX_CHUNK_SIZE { - // max_chunk_size = M2M_MAX_CHUNK_SIZE - // } - // let chunk_size = find_chunk_size(nparents, max_chunk_size); + let nsiblings = 8; + + // 1. Lookup parents and corresponding children that exist for this set of sources + // Must explicitly lookup as boxes may be empty at this level, and the next. + let parent_targets: HashSet = + child_sources.iter().map(|source| source.parent()).collect(); + let mut parent_targets = parent_targets.into_iter().collect_vec(); + parent_targets.sort(); + let nparents = parent_targets.len(); + let mut parent_multipoles = vec![Vec::new(); nparents]; + + for (parent_idx, parent) in parent_targets.iter().enumerate() { + for charge_vec_idx in 0..self.ncharge_vectors { + let parent_index_pointer = *self.level_index_pointer[(level - 1) as usize] + .get(parent) + .unwrap(); + let parent_multipole = + self.level_multipoles[(level - 1) as usize][parent_index_pointer][charge_vec_idx]; + parent_multipoles[parent_idx].push(parent_multipole); + } + } + + let n_child_sources = child_sources.len(); + let min: &MortonKey = &child_sources[0]; + let max = &child_sources[n_child_sources - 1]; + let min_idx = self.fmm.tree.get_index(min).unwrap(); + let min_key_displacement = min_idx * self.ncoeffs * self.ncharge_vectors; + let max_idx = self.fmm.tree().get_index(max).unwrap(); + let max_key_displacement = (max_idx + 1) * self.ncoeffs * self.ncharge_vectors; + let child_multipoles = &self.multipoles[min_key_displacement..max_key_displacement]; + + child_multipoles + .par_chunks_exact(self.ncharge_vectors * self.ncoeffs * nsiblings) + .zip(parent_multipoles.into_par_iter()) + .for_each(|(child_multipoles, parent_multipole_pointers)| { + + for i in 0..nsiblings { + let sibling_displacement = i * self.ncoeffs * self.ncharge_vectors; + let ptr = unsafe { child_multipoles.as_ptr().add(sibling_displacement) }; + let child_multipoles_i = unsafe { rlst_pointer_mat!['a, V, ptr, (self.ncoeffs, self.ncharge_vectors), (1, self.ncoeffs)] }; + let parent_multipole_i = self.fmm.m2m[i].dot(&child_multipoles_i).eval(); + + for j in 0..self.ncharge_vectors { + let raw = parent_multipole_pointers[j].raw; + let parent_multipole_j = unsafe { std::slice::from_raw_parts_mut(raw, self.ncoeffs) }; + let parent_multipole_ij = &parent_multipole_i.data()[j*self.ncoeffs..(j+1)*self.ncoeffs]; + parent_multipole_j.iter_mut().zip(parent_multipole_ij.iter()).for_each(|(p, r)| *p += *r); + } + } + }); // // 3. Compute M2M kernel over sets of siblings // child_multipoles @@ -728,7 +745,6 @@ mod test { let order = 8; let alpha_inner = 1.05; let alpha_outer = 2.95; - let ncrit = 150; let depth = 3; let adaptive = false; @@ -738,8 +754,8 @@ mod test { let tree = SingleNodeTree::new( points.data(), adaptive, - Some(ncrit), None, + Some(depth), &global_idxs[..], false, ); @@ -747,7 +763,7 @@ mod test { // Precompute the M2L data let m2l_data = FftFieldTranslationKiFmm::new(kernel.clone(), order, *tree.get_domain(), alpha_inner); - let fmm = KiFmmLinear::new(order, alpha_inner, alpha_outer, kernel, tree, m2l_data); + let fmm = KiFmmLinearMatrix::new(order, alpha_inner, alpha_outer, kernel, tree, m2l_data); // Form charge dict, matching charges with their associated global indices let mut charge_dicts = Vec::new(); @@ -767,6 +783,7 @@ mod test { break; } } + let leaf = &datatree.fmm.tree().get_all_leaves().unwrap()[test_idx]; let &leaf_idx = datatree.fmm.tree.get_leaf_index(leaf).unwrap(); let multipoles = &datatree.leaf_multipoles[leaf_idx]; @@ -788,12 +805,11 @@ mod test { rlst_pointer_mat!['static, f64, leaf_coordinates.as_ptr(), (nsources, datatree.fmm.kernel.space_dimension()), (datatree.fmm.kernel.space_dimension(), 1)] }.eval(); - let kernel = Laplace3dKernel::::default(); for i in 0..ncharge_vecs { let charge_vec_displacement = i * ncoordinates; let charges = &datatree.charges[charge_vec_displacement+l..charge_vec_displacement+r]; - kernel.evaluate_st( + datatree.fmm.kernel.evaluate_st( EvalType::Value, leaf_coordinates.data(), &test_point, @@ -804,7 +820,7 @@ mod test { for i in 0..ncharge_vecs { let multipole = unsafe { std::slice::from_raw_parts(multipoles[i].raw, datatree.ncoeffs) }; - kernel.evaluate_st( + datatree.fmm.kernel.evaluate_st( EvalType::Value, &upward_equivalent_surface, &test_point, @@ -818,6 +834,94 @@ mod test { } } + + #[test] + fn test_upward_pass_uniform_matrix() { + let npoints = 10000; + let ncharge_vecs = 10; + let points = points_fixture::(npoints, None, None); + let global_idxs = (0..npoints).collect_vec(); + let mut charge_mat = vec![vec![0.0; npoints]; ncharge_vecs]; + for i in 0..ncharge_vecs { + charge_mat[i] = vec![i as f64 + 1.0; npoints] + } + + let order = 8; + let alpha_inner = 1.05; + let alpha_outer = 2.95; + let depth = 3; + let adaptive = false; + + let kernel = Laplace3dKernel::default(); + + // Create a tree + let tree = SingleNodeTree::new( + points.data(), + adaptive, + None, + Some(depth), + &global_idxs[..], + false, + ); + + // Precompute the M2L data + let m2l_data = + FftFieldTranslationKiFmm::new(kernel.clone(), order, *tree.get_domain(), alpha_inner); + let fmm = KiFmmLinearMatrix::new(order, alpha_inner, alpha_outer, kernel, tree, m2l_data); + + // Form charge dict, matching charges with their associated global indices + let mut charge_dicts = Vec::new(); + for i in 0..ncharge_vecs { + charge_dicts.push(build_charge_dict(&global_idxs, &charge_mat[i])) + } + // Associate data with the FMM + let datatree = FmmDataUniformMatrix::new(fmm, &charge_dicts).unwrap(); + + // Upward pass + { + datatree.p2m(); + + for level in (1..=depth).rev() { + datatree.m2m(level); + } + } + + let &midx = datatree.fmm.tree().get_index(&ROOT).unwrap(); + let multipoles = &datatree.level_multipoles[ROOT.level() as usize][0]; + + let upward_equivalent_surface = + ROOT.compute_surface(&datatree.fmm.tree().domain, order, datatree.fmm.alpha_inner); + + let test_point = vec![100000., 0., 0.]; + + let mut expected = vec![0.; datatree.ncharge_vectors]; + let mut found = vec![0.; datatree.ncharge_vectors]; + + for i in 0..ncharge_vecs { + datatree.fmm.kernel().evaluate_st( + EvalType::Value, + points.data(), + &test_point, + &charge_mat[i], + &mut expected[i..i+1], + ); + } + + for i in 0..ncharge_vecs { + let multipole = unsafe { std::slice::from_raw_parts(multipoles[i].raw, datatree.ncoeffs) }; + datatree.fmm.kernel().evaluate_st( + EvalType::Value, + &upward_equivalent_surface, + &test_point, + multipole, + &mut found[i..i+1], + ); + } + + for (&a, &b) in expected.iter().zip(found.iter()) { + assert_approx_eq!(f64, a, b, epsilon = 1e-5); + } + } #[test] fn test_upward_pass_sphere_adaptive() { diff --git a/fmm/src/fmm.rs b/fmm/src/fmm.rs index 2ffc9076..a9a3ebfa 100644 --- a/fmm/src/fmm.rs +++ b/fmm/src/fmm.rs @@ -7,7 +7,7 @@ use std::time::Instant; use rlst::{ algorithms::{linalg::DenseMatrixLinAlgBuilder, traits::svd::Svd}, common::traits::{Eval, Transpose}, - dense::{rlst_dynamic_mat, rlst_pointer_mat, traits::*, Dot, MultiplyAdd, VectorContainer}, + dense::{rlst_dynamic_mat, rlst_pointer_mat, traits::*, Dot, MultiplyAdd, VectorContainer, DataContainer}, }; use bempp_traits::{ @@ -19,7 +19,7 @@ use bempp_traits::{ }; use bempp_tree::{constants::ROOT, types::single_node::SingleNodeTree}; -use crate::types::{FmmDataAdaptive, FmmDataUniform, FmmDataUniformMatrix}; +use crate::types::{FmmDataAdaptive, FmmDataUniform, FmmDataUniformMatrix, KiFmmLinearMatrix}; use crate::{ pinv::{pinv, SvdScalar}, types::KiFmmLinear, @@ -291,6 +291,270 @@ where } } +/// Implementation of constructor for single node KiFMM +impl<'a, T, U, V> KiFmmLinearMatrix, T, U, V> +where + T: Kernel + ScaleInvariantKernel, + U: FieldTranslationData, + V: Scalar + Default + Float, + SvdScalar: PartialOrd, + SvdScalar: Scalar + Float + ToPrimitive, + DenseMatrixLinAlgBuilder: Svd, + V: MultiplyAdd< + V, + VectorContainer, + VectorContainer, + VectorContainer, + Dynamic, + Dynamic, + Dynamic, + >, + SvdScalar: MultiplyAdd< + SvdScalar, + VectorContainer>, + VectorContainer>, + VectorContainer>, + Dynamic, + Dynamic, + Dynamic, + >, +{ + /// Constructor for single node kernel independent FMM (KiFMM). This object contains all the precomputed operator matrices and metadata, as well as references to + /// the associated single node octree, and the associated kernel function. + /// + /// # Arguments + /// * `order` - The expansion order for the multipole and local expansions. + /// * `alpha_inner` - The ratio of the inner check surface diamater in comparison to the surface discretising a box. + /// * `alpha_order` - The ratio of the outer check surface diamater in comparison to the surface discretising a box. + /// * `kernel` - The kernel function for this FMM. + /// * `tree` - The type of tree associated with this FMM, can be single or multi node. + /// * `m2l` - The M2L operator matrices, as well as metadata associated with this FMM. + pub fn new( + order: usize, + alpha_inner: V, + alpha_outer: V, + kernel: T, + tree: SingleNodeTree, + m2l: U, + ) -> Self { + let upward_equivalent_surface = ROOT.compute_surface(tree.get_domain(), order, alpha_inner); + let upward_check_surface = ROOT.compute_surface(tree.get_domain(), order, alpha_outer); + let downward_equivalent_surface = + ROOT.compute_surface(tree.get_domain(), order, alpha_outer); + let downward_check_surface = ROOT.compute_surface(tree.get_domain(), order, alpha_inner); + + let nequiv_surface = upward_equivalent_surface.len() / kernel.space_dimension(); + let ncheck_surface = upward_check_surface.len() / kernel.space_dimension(); + + // Store in RLST matrices + let upward_equivalent_surface = unsafe { + rlst_pointer_mat!['a, ::Real, upward_equivalent_surface.as_ptr(), (nequiv_surface, kernel.space_dimension()), (1, nequiv_surface)] + }; + let upward_check_surface = unsafe { + rlst_pointer_mat!['a, ::Real, upward_check_surface.as_ptr(), (ncheck_surface, kernel.space_dimension()), (1, ncheck_surface)] + }; + let downward_equivalent_surface = unsafe { + rlst_pointer_mat!['a, ::Real, downward_equivalent_surface.as_ptr(), (nequiv_surface, kernel.space_dimension()), (1, nequiv_surface)] + }; + let downward_check_surface = unsafe { + rlst_pointer_mat!['a, ::Real, downward_check_surface.as_ptr(), (ncheck_surface, kernel.space_dimension()), (1, ncheck_surface)] + }; + + // Compute upward check to equivalent, and downward check to equivalent Gram matrices + // as well as their inverses using DGESVD. + let mut uc2e = rlst_dynamic_mat![V, (ncheck_surface, nequiv_surface)]; + kernel.assemble_st( + EvalType::Value, + upward_equivalent_surface.data(), + upward_check_surface.data(), + uc2e.data_mut(), + ); + + // Need to tranapose so that rows correspond to targets and columns to sources + let uc2e = uc2e.transpose().eval(); + + let mut dc2e = rlst_dynamic_mat![V, (ncheck_surface, nequiv_surface)]; + kernel.assemble_st( + EvalType::Value, + downward_equivalent_surface.data(), + downward_check_surface.data(), + dc2e.data_mut(), + ); + + // Need to tranapose so that rows correspond to targets and columns to sources + let dc2e = dc2e.transpose().eval(); + + let (s, ut, v) = pinv::(&uc2e, None, None).unwrap(); + + let mut mat_s = rlst_dynamic_mat![SvdScalar, (s.len(), s.len())]; + for i in 0..s.len() { + mat_s[[i, i]] = SvdScalar::::from_real(s[i]); + } + let uc2e_inv_1 = v.dot(&mat_s); + let uc2e_inv_2 = ut; + + let uc2e_inv_1_shape = uc2e_inv_1.shape(); + let uc2e_inv_2_shape = uc2e_inv_2.shape(); + + let uc2e_inv_1 = uc2e_inv_1 + .data() + .iter() + .map(|x| V::from(*x).unwrap()) + .collect_vec(); + let uc2e_inv_1 = unsafe { + rlst_pointer_mat!['a, V, uc2e_inv_1.as_ptr(), uc2e_inv_1_shape, (1, uc2e_inv_1_shape.0)] + } + .eval(); + let uc2e_inv_2 = uc2e_inv_2 + .data() + .iter() + .map(|x| V::from(*x).unwrap()) + .collect_vec(); + let uc2e_inv_2 = unsafe { + rlst_pointer_mat!['a, V, uc2e_inv_2.as_ptr(), uc2e_inv_2_shape, (1, uc2e_inv_2_shape.0)] + } + .eval(); + + let (s, ut, v) = pinv::(&dc2e, None, None).unwrap(); + + let mut mat_s = rlst_dynamic_mat![SvdScalar, (s.len(), s.len())]; + for i in 0..s.len() { + mat_s[[i, i]] = SvdScalar::::from_real(s[i]); + } + + let dc2e_inv_1 = v.dot(&mat_s); + let dc2e_inv_2 = ut; + + let dc2e_inv_1_shape = dc2e_inv_1.shape(); + let dc2e_inv_2_shape = dc2e_inv_2.shape(); + + let dc2e_inv_1 = dc2e_inv_1 + .data() + .iter() + .map(|x| V::from(*x).unwrap()) + .collect_vec(); + let dc2e_inv_1 = unsafe { + rlst_pointer_mat!['a, V, dc2e_inv_1.as_ptr(), dc2e_inv_1_shape, (1, dc2e_inv_1_shape.0)] + } + .eval(); + let dc2e_inv_2 = dc2e_inv_2 + .data() + .iter() + .map(|x| V::from(*x).unwrap()) + .collect_vec(); + let dc2e_inv_2 = unsafe { + rlst_pointer_mat!['a, V, dc2e_inv_2.as_ptr(), dc2e_inv_2_shape, (1, dc2e_inv_2_shape.0)] + } + .eval(); + + // Calculate M2M/L2L matrices + let children = ROOT.children(); + let mut m2m = Vec::new(); + let mut l2l = Vec::new(); + + for (i, child) in children.iter().enumerate() { + let child_upward_equivalent_surface = + child.compute_surface(tree.get_domain(), order, alpha_inner); + let child_downward_check_surface = + child.compute_surface(tree.get_domain(), order, alpha_inner); + let child_upward_equivalent_surface = unsafe { + rlst_pointer_mat!['a, ::Real, child_upward_equivalent_surface.as_ptr(), (nequiv_surface, kernel.space_dimension()), (1, nequiv_surface)] + }; + let child_downward_check_surface = unsafe { + rlst_pointer_mat!['a, ::Real, child_downward_check_surface.as_ptr(), (ncheck_surface, kernel.space_dimension()), (1, ncheck_surface)] + }; + + let mut pc2ce = rlst_dynamic_mat![V, (ncheck_surface, nequiv_surface)]; + + kernel.assemble_st( + EvalType::Value, + child_upward_equivalent_surface.data(), + upward_check_surface.data(), + pc2ce.data_mut(), + ); + + // Need to transpose so that rows correspond to targets, and columns to sources + let pc2ce = pc2ce.transpose().eval(); + + let tmp = uc2e_inv_1.dot(&uc2e_inv_2.dot(&pc2ce)).eval(); + m2m.push(tmp); + + let mut cc2pe = rlst_dynamic_mat![V, (ncheck_surface, nequiv_surface)]; + + kernel.assemble_st( + EvalType::Value, + downward_equivalent_surface.data(), + child_downward_check_surface.data(), + cc2pe.data_mut(), + ); + + // Need to transpose so that rows correspond to targets, and columns to sources + let cc2pe = cc2pe.transpose().eval(); + let mut tmp = dc2e_inv_1.dot(&dc2e_inv_2.dot(&cc2pe)).eval(); + tmp.data_mut() + .iter_mut() + .for_each(|d| *d *= kernel.scale(child.level())); + + l2l.push(tmp); + } + + Self { + order, + uc2e_inv_1, + uc2e_inv_2, + dc2e_inv_1, + dc2e_inv_2, + alpha_inner, + alpha_outer, + m2m, + l2l, + kernel, + tree, + m2l, + } + } +} + +impl KiFmmTrait for KiFmmLinearMatrix +where + T: Tree, + U: Kernel, + V: FieldTranslationData, + W: Scalar + Float + Default, +{ + fn alpha_inner(&self) -> <::Kernel as Kernel>::T { + self.alpha_inner + } + + fn alpha_outer(&self) -> <::Kernel as Kernel>::T { + self.alpha_outer + } +} + +impl Fmm for KiFmmLinearMatrix +where + T: Tree, + U: Kernel, + V: FieldTranslationData, + W: Scalar + Float + Default, +{ + type Tree = T; + type Kernel = U; + + fn order(&self) -> usize { + self.order + } + + fn kernel(&self) -> &Self::Kernel { + &self.kernel + } + + fn tree(&self) -> &Self::Tree { + &self.tree + } +} + + impl FmmLoop for FmmDataAdaptive where T: Fmm, diff --git a/fmm/src/types.rs b/fmm/src/types.rs index 32a40a0a..b826b0b3 100644 --- a/fmm/src/types.rs +++ b/fmm/src/types.rs @@ -229,11 +229,17 @@ where /// The pseudo-inverse of the dense interaction matrix between the upward check and upward equivalent surfaces. /// Store in two parts to avoid propagating error from computing pseudo-inverse pub uc2e_inv_1: C2EType, + + /// The pseudo-inverse of the dense interaction matrix between the upward check and upward equivalent surfaces. + /// Store in two parts to avoid propagating error from computing pseudo-inverse pub uc2e_inv_2: C2EType, /// The pseudo-inverse of the dense interaction matrix between the downward check and downward equivalent surfaces. /// Store in two parts to avoid propagating error from computing pseudo-inverse pub dc2e_inv_1: C2EType, + + /// The pseudo-inverse of the dense interaction matrix between the downward check and downward equivalent surfaces. + /// Store in two parts to avoid propagating error from computing pseudo-inverse pub dc2e_inv_2: C2EType, /// The ratio of the inner check surface diamater in comparison to the surface discretising a box. @@ -258,6 +264,57 @@ where pub m2l: V, } + +/// Type to store data associated with the kernel independent (KiFMM) in for matrix input FMMs. +pub struct KiFmmLinearMatrix +where + T: Tree, + U: Kernel, + V: FieldTranslationData, + W: Scalar + Float + Default, +{ + /// The expansion order + pub order: usize, + + /// The pseudo-inverse of the dense interaction matrix between the upward check and upward equivalent surfaces. + /// Store in two parts to avoid propagating error from computing pseudo-inverse + pub uc2e_inv_1: C2EType, + + /// The pseudo-inverse of the dense interaction matrix between the upward check and upward equivalent surfaces. + /// Store in two parts to avoid propagating error from computing pseudo-inverse + pub uc2e_inv_2: C2EType, + + /// The pseudo-inverse of the dense interaction matrix between the downward check and downward equivalent surfaces. + /// Store in two parts to avoid propagating error from computing pseudo-inverse + pub dc2e_inv_1: C2EType, + + /// The pseudo-inverse of the dense interaction matrix between the downward check and downward equivalent surfaces. + /// Store in two parts to avoid propagating error from computing pseudo-inverse + pub dc2e_inv_2: C2EType, + + /// The ratio of the inner check surface diamater in comparison to the surface discretising a box. + pub alpha_inner: W, + + /// The ratio of the outer check surface diamater in comparison to the surface discretising a box. + pub alpha_outer: W, + + /// The multipole to multipole operator matrices, each index is associated with a child box (in sequential Morton order), + pub m2m: Vec>, + + /// The local to local operator matrices, each index is associated with a child box (in sequential Morton order). + pub l2l: Vec>, + + /// The tree (single or multi node) associated with this FMM + pub tree: T, + + /// The kernel associated with this FMM. + pub kernel: U, + + /// The M2L operator matrices, as well as metadata associated with this FMM. + pub m2l: V, +} + + /// A threadsafe mutable raw pointer #[derive(Clone, Debug, Copy)] pub struct SendPtrMut { @@ -466,7 +523,7 @@ where } /// Implementation of the data structure to store the data for the single node KiFMM. -impl FmmDataUniformMatrix, T, U, V>, V> +impl FmmDataUniformMatrix, T, U, V>, V> where T: Kernel + ScaleInvariantKernel, U: FieldTranslationData, @@ -478,7 +535,7 @@ where /// `fmm` - A single node KiFMM object. /// `global_charges` - The charge data associated to the point data via unique global indices. pub fn new( - fmm: KiFmmLinear, T, U, V>, + fmm: KiFmmLinearMatrix, T, U, V>, global_charges: &[ChargeDict], ) -> Result { if let Some(keys) = fmm.tree().get_all_keys() { @@ -903,7 +960,7 @@ where mod test { use crate::{ charge::build_charge_dict, - types::{FmmDataUniform, FmmDataUniformMatrix, KiFmmLinear}, + types::{FmmDataUniform, FmmDataUniformMatrix, KiFmmLinear, KiFmmLinearMatrix}, }; use bempp_field::types::FftFieldTranslationKiFmm; use bempp_kernel::laplace_3d::Laplace3dKernel; @@ -954,7 +1011,7 @@ mod test { alpha_inner, ); - let fmm = KiFmmLinear::new(order, alpha_inner, alpha_outer, kernel, tree, m2l_data); + let fmm = KiFmmLinearMatrix::new(order, alpha_inner, alpha_outer, kernel, tree, m2l_data); // Form charge dicts, matching all charges with their associated global indices let mut charge_dicts = Vec::new();