diff --git a/fmm/Cargo.toml b/fmm/Cargo.toml index 7f7be620..96f3a658 100644 --- a/fmm/Cargo.toml +++ b/fmm/Cargo.toml @@ -35,6 +35,9 @@ rand = "0.8.*" float-cmp = "0.9.0" num_cpus = "1" num = "0.4" -rlst = {git = "https://github.com/skailasa/rlst.git", branch = "enh/moore-penrose-pseudo-inverse" } +rlst = { git = "https://github.com/linalg-rs/rlst.git" } fftw = {git = "https://github.com/skailasa/fftw.git" } rayon = "1.7" + +[dev.dependencies] +rlst = { git = "https://github.com/linalg-rs/rlst.git" } \ No newline at end of file diff --git a/fmm/src/field_translation.rs b/fmm/src/field_translation.rs index 02779ac2..e1336e97 100644 --- a/fmm/src/field_translation.rs +++ b/fmm/src/field_translation.rs @@ -28,7 +28,7 @@ use bempp_tree::types::{morton::MortonKey, single_node::SingleNodeTree}; use rlst::{ common::traits::*, - dense::{rlst_col_vec, rlst_mat, rlst_pointer_mat, traits::*, Dot, Shape}, + dense::{rlst_col_vec, rlst_dynamic_mat, rlst_pointer_mat, traits::*, Dot, Shape}, }; use crate::types::{FmmData, KiFmm}; @@ -426,7 +426,7 @@ where let c_sub = self.fmm.m2l.operator_data.c.block(top_left, dim); let m2l_rw = m2l_arc.read().unwrap(); - let mut multipoles = rlst_mat![f64, (self.fmm.m2l.k, m2l_rw.len())]; + let mut multipoles = rlst_dynamic_mat![f64, (self.fmm.m2l.k, m2l_rw.len())]; for (i, (source, _)) in m2l_rw.iter().enumerate() { let source_multipole_arc = Arc::clone(self.multipoles.get(source).unwrap()); diff --git a/fmm/src/fmm.rs b/fmm/src/fmm.rs index e715d214..972b60a8 100644 --- a/fmm/src/fmm.rs +++ b/fmm/src/fmm.rs @@ -7,9 +7,8 @@ use std::{ }; use rlst::{ - algorithms::{linalg::LinAlg, traits::pseudo_inverse::Pinv}, common::traits::{Eval, NewLikeSelf, Transpose}, - dense::{rlst_col_vec, rlst_mat, rlst_pointer_mat, traits::*, Dot}, + dense::{rlst_col_vec, rlst_dynamic_mat, rlst_pointer_mat, traits::*, Dot}, }; use bempp_traits::{ @@ -21,6 +20,7 @@ use bempp_traits::{ }; use bempp_tree::{constants::ROOT, types::single_node::SingleNodeTree}; +use crate::pinv::pinv; use crate::types::{C2EType, ChargeDict, FmmData, KiFmm}; /// Implementation of constructor for single node KiFMM @@ -72,7 +72,7 @@ where // Compute upward check to equivalent, and downward check to equivalent Gram matrices // as well as their inverses using DGESVD. - let mut uc2e = rlst_mat![f64, (ncheck_surface, nequiv_surface)]; + let mut uc2e = rlst_dynamic_mat![f64, (ncheck_surface, nequiv_surface)]; kernel.assemble_st( EvalType::Value, upward_equivalent_surface.data(), @@ -83,7 +83,7 @@ where // Need to tranapose so that rows correspond to targets and columns to sources let uc2e = uc2e.transpose().eval(); - let mut dc2e = rlst_mat![f64, (ncheck_surface, nequiv_surface)]; + let mut dc2e = rlst_dynamic_mat![f64, (ncheck_surface, nequiv_surface)]; kernel.assemble_st( EvalType::Value, downward_equivalent_surface.data(), @@ -94,21 +94,17 @@ where // Need to tranapose so that rows correspond to targets and columns to sources let dc2e = dc2e.transpose().eval(); - let (s, ut, v) = uc2e.linalg().pinv(None).unwrap(); - let s = s.unwrap(); - let ut = ut.unwrap(); - let v = v.unwrap(); - let mut mat_s = rlst_mat![f64, (s.len(), s.len())]; + let (s, ut, v) = pinv::(&uc2e, None, None).unwrap(); + + let mut mat_s = rlst_dynamic_mat![f64, (s.len(), s.len())]; for i in 0..s.len() { mat_s[[i, i]] = s[i]; } let uc2e_inv = v.dot(&mat_s).dot(&ut); - let (s, ut, v) = dc2e.linalg().pinv(None).unwrap(); - let s = s.unwrap(); - let ut = ut.unwrap(); - let v = v.unwrap(); - let mut mat_s = rlst_mat![f64, (s.len(), s.len())]; + let (s, ut, v) = pinv::(&dc2e, None, None).unwrap(); + + let mut mat_s = rlst_dynamic_mat![f64, (s.len(), s.len())]; for i in 0..s.len() { mat_s[[i, i]] = s[i]; } @@ -131,7 +127,7 @@ where rlst_pointer_mat!['a, f64, child_downward_check_surface.as_ptr(), (ncheck_surface, kernel.space_dimension()), (1, ncheck_surface)] }; - let mut pc2ce = rlst_mat![f64, (ncheck_surface, nequiv_surface)]; + let mut pc2ce = rlst_dynamic_mat![f64, (ncheck_surface, nequiv_surface)]; kernel.assemble_st( EvalType::Value, @@ -145,7 +141,7 @@ where m2m.push(uc2e_inv.dot(&pc2ce).eval()); - let mut cc2pe = rlst_mat![f64, (ncheck_surface, nequiv_surface)]; + let mut cc2pe = rlst_dynamic_mat![f64, (ncheck_surface, nequiv_surface)]; kernel.assemble_st( EvalType::Value, @@ -394,7 +390,7 @@ mod test { let charges = vec![1.0; npoints]; // Setup a FMM experiment - let order = 9; + let order = 6; let alpha_inner = 1.05; let alpha_outer = 1.99; let adaptive = false; @@ -477,8 +473,7 @@ mod test { .sum(); let rel_error: f64 = abs_error / (direct.iter().sum::()); - println!("{:?}", rel_error); - assert!(rel_error <= 1e-6); + assert!(rel_error <= 1e-5); } #[test] @@ -488,7 +483,7 @@ mod test { let global_idxs = (0..npoints).collect_vec(); let charges = vec![1.0; npoints]; - let order = 9; + let order = 6; let alpha_inner = 1.05; let alpha_outer = 1.95; let adaptive = false; @@ -514,9 +509,7 @@ mod test { let datatree = FmmData::new(fmm, &charge_dict); - let times = datatree.run(Some(true)); - - println!("FFT times {:?}", times); + datatree.run(Some(true)); let leaf = &datatree.fmm.tree.get_leaves().unwrap()[0]; @@ -558,7 +551,6 @@ mod test { .map(|(a, b)| (a - b).abs()) .sum(); let rel_error: f64 = abs_error / (direct.iter().sum::()); - println!("rel error {:?}", rel_error); - assert!(rel_error <= 1e-6); + assert!(rel_error <= 1e-5); } } diff --git a/fmm/src/lib.rs b/fmm/src/lib.rs index fb715141..fbc31a22 100644 --- a/fmm/src/lib.rs +++ b/fmm/src/lib.rs @@ -4,4 +4,5 @@ pub mod constants; pub mod field_translation; pub mod fmm; pub mod interaction_lists; +pub mod pinv; pub mod types; diff --git a/fmm/src/pinv.rs b/fmm/src/pinv.rs new file mode 100644 index 00000000..a0b00498 --- /dev/null +++ b/fmm/src/pinv.rs @@ -0,0 +1,117 @@ +//! Implementation of Moore-Penrose PseudoInverse +use num::{Float, Zero}; +use rlst::algorithms::linalg::LinAlg; +use rlst::algorithms::traits::svd::{Mode, Svd}; +use rlst::dense::{ + base_matrix::BaseMatrix, data_container::VectorContainer, matrix::Matrix, Dynamic, Shape, +}; +// use rlst_common::traits::*; +use rlst::common::traits::{Eval, Transpose}; +use rlst::common::types::{RlstError, RlstResult, Scalar}; +use rlst::dense::MatrixD; + +pub type PinvMatrix = Matrix, Dynamic>, Dynamic>; + +type PinvReturnType = RlstResult<( + Vec<::Real>, + MatrixD<::Real>, + MatrixD<::Real>, +)>; + +/// Compute the (Moore-Penrose) pseudo-inverse of a matrix. +/// +/// Calculate a generalised inverse using its singular value decomposition `U @ S @ V*`. +/// If `s` is the maximum singular value, then the signifance cut-off value is determined by +/// `atol + rtol * s`. Any singular value below this is assumed insignificant. +/// +/// # Arguments +/// * `mat` - (M, N) matrix to be inverted. +/// * `atol` - Absolute threshold term, default is 0. +/// * `rtol` - Relative threshold term, default value is max(M, N) * eps +pub fn pinv + Float>( + mat: &PinvMatrix, + atol: Option, + rtol: Option, +) -> PinvReturnType { + let shape = mat.shape(); + + if shape.0 == 0 || shape.1 == 0 { + return Err(RlstError::MatrixIsEmpty(shape)); + } + + // If we have a vector return error + if shape.0 == 1 || shape.1 == 1 { + Err(RlstError::SingleDimensionError { + expected: 2, + actual: 1, + }) + } else { + // For matrices compute the full SVD + let (mut s, u, vt) = mat.linalg().svd(Mode::All, Mode::All)?; + let u = u.unwrap(); + let vt = vt.unwrap(); + + let eps = T::real(T::epsilon()); + let max_dim = T::real(std::cmp::max(shape.0, shape.1)); + + let atol = atol.unwrap_or(T::Real::zero()); + let rtol = rtol.unwrap_or(max_dim * eps); + + let max_s = T::real(s[0]); + let threshold = atol + rtol * max_s; + // Filter singular values below this threshold + for s in s.iter_mut() { + if *s > threshold { + *s = T::real(1.0) / *s; + } else { + *s = T::real(0.) + } + } + + // Return pseudo-inverse in component form + let v = vt.transpose().eval(); + let ut = u.transpose().eval(); + + Ok((s, ut, v)) + } +} + +#[cfg(test)] +mod test { + + use super::*; + use approx::assert_relative_eq; + use rlst::common::traits::ColumnMajorIterator; + use rlst::common::traits::NewLikeSelf; + use rlst::dense::{rlst_dynamic_mat, rlst_rand_mat, Dot}; + + #[test] + fn test_pinv() { + let dim: usize = 5; + let mat = rlst_rand_mat![f64, (dim, dim)]; + + let (s, ut, v) = pinv::(&mat, None, None).unwrap(); + + let ut = ut; + let v = v; + + let mut mat_s = rlst_dynamic_mat![f64, (s.len(), s.len())]; + for i in 0..s.len() { + mat_s[[i, i]] = s[i]; + } + + let inv = v.dot(&mat_s).dot(&ut); + + let actual = inv.dot(&mat); + + // Expect the identity matrix + let mut expected = actual.new_like_self(); + for i in 0..dim { + expected[[i, i]] = 1.0 + } + + for (a, e) in actual.iter_col_major().zip(expected.iter_col_major()) { + assert_relative_eq!(a, e, epsilon = 1E-13); + } + } +} diff --git a/fmm/src/types.rs b/fmm/src/types.rs index 3232ffc5..08f6b85e 100644 --- a/fmm/src/types.rs +++ b/fmm/src/types.rs @@ -19,16 +19,13 @@ pub type GlobalIdx = usize; pub type ChargeDict = HashMap; /// Type alias for multipole/local expansion containers. -pub type Expansions = - Matrix, Dynamic, Dynamic>, Dynamic, Dynamic>; +pub type Expansions = Matrix, Dynamic>, Dynamic>; /// Type alias for potential containers. -pub type Potentials = - Matrix, Dynamic, Dynamic>, Dynamic, Dynamic>; +pub type Potentials = Matrix, Dynamic>, Dynamic>; /// Type alias for approximation of FMM operator matrices. -pub type C2EType = - Matrix, Dynamic, Dynamic>, Dynamic, Dynamic>; +pub type C2EType = Matrix, Dynamic>, Dynamic>; /// Type to store data associated with an FMM in. pub struct FmmData {