Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

FIX: Accuracy Issue Due to Pseudo Inverse #129

Merged
merged 3 commits into from
Nov 1, 2023
Merged

FIX: Accuracy Issue Due to Pseudo Inverse #129

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1fc7dd0
Fix accuracy issue
skailasa Nov 1, 2023
426d0c2
Fix accuracy issue
skailasa Nov 1, 2023
11cb6c4
Remove unused manifest for dev
skailasa Nov 1, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion field/src/field.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -440,7 +440,7 @@ where
///
/// # Arguments
/// * `order` - The expansion order for the multipole and local expansions.
/// * `convolution_grid` - Cartesian coordinates of points on the convolution grid at a source box, expected in row major order.
/// * `convolution_grid` - Cartesian coordinates of points on the convolution grid at a source box, expected in column major order.
/// * `target_pt` - The point on the target box's surface grid, with which kernels are being evaluated with respect to.
pub fn compute_kernel(
&self,
Expand Down
3 changes: 0 additions & 3 deletions fmm/Cargo.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,3 @@ num = "0.4"
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" }
1 change: 0 additions & 1 deletion fmm/examples/test_fmm.rs
View file
Open in desktop

This file was deleted.

23 changes: 12 additions & 11 deletions fmm/src/field_translation.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,6 @@ where

let nsources = leaf_coordinates.len() / self.fmm.kernel.space_dimension();

// Get into row major order
let leaf_coordinates = unsafe {
rlst_pointer_mat!['a, f64, leaf_coordinates.as_ptr(), (nsources, fmm_arc.kernel.space_dimension()), (fmm_arc.kernel.space_dimension(), 1)]
}.eval();
Expand All @@ -84,7 +83,7 @@ where

let leaf_multipole_owned = (
fmm_arc.kernel.scale(leaf.level())
* fmm_arc.uc2e_inv.dot(&check_potential)
* fmm_arc.uc2e_inv_1.dot(&fmm_arc.uc2e_inv_2.dot(&check_potential))
).eval();

let mut leaf_multipole_lock = leaf_multipole_arc.lock().unwrap();
Expand Down Expand Up @@ -174,7 +173,6 @@ where

let ntargets = target_coordinates.len() / self.fmm.kernel.space_dimension();

// Get into row major order
let target_coordinates = unsafe {
rlst_pointer_mat!['a, f64, target_coordinates.as_ptr(), (ntargets, fmm_arc.kernel.space_dimension()), (fmm_arc.kernel.space_dimension(), 1)]
}.eval();
Expand Down Expand Up @@ -216,7 +214,6 @@ where
.collect_vec();
let ntargets = target_coordinates.len() / self.fmm.kernel.space_dimension();

// Get into row major order
let target_coordinates = unsafe {
rlst_pointer_mat!['a, f64, target_coordinates.as_ptr(), (ntargets, fmm_arc.kernel.space_dimension()), (fmm_arc.kernel.space_dimension(), 1)]
}.eval();
Expand Down Expand Up @@ -264,7 +261,6 @@ where

let nsources = source_coordinates.len() / self.fmm.kernel.space_dimension();

// Get into row major order
let source_coordinates = unsafe {
rlst_pointer_mat!['a, f64, source_coordinates.as_ptr(), (nsources, fmm_arc.kernel.space_dimension()), (fmm_arc.kernel.space_dimension(), 1)]
}.eval();
Expand All @@ -291,7 +287,7 @@ where

let mut target_local_lock = target_local_arc.lock().unwrap();

let target_local_owned = (fmm_arc.kernel.scale(leaf.level()) * fmm_arc.dc2e_inv.dot(&downward_check_potential)).eval();
let target_local_owned = (fmm_arc.kernel.scale(leaf.level()) * fmm_arc.dc2e_inv_1.dot(&fmm_arc.dc2e_inv_2.dot(&downward_check_potential))).eval();

*target_local_lock.deref_mut() = (target_local_lock.deref() + target_local_owned).eval();
}
Expand All @@ -316,7 +312,6 @@ where

let ntargets= target_coordinates.len() / self.fmm.kernel.space_dimension();

// Get into row major order
let target_coordinates = unsafe {
rlst_pointer_mat!['a, f64, target_coordinates.as_ptr(), (ntargets, fmm_arc.kernel.space_dimension()), (fmm_arc.kernel.space_dimension(), 1)]
}.eval();
Expand All @@ -332,7 +327,6 @@ where

let nsources = source_coordinates.len() / self.fmm.kernel.space_dimension();

// Get into row major order
let source_coordinates = unsafe {
rlst_pointer_mat!['a, f64, source_coordinates.as_ptr(), (nsources, fmm_arc.kernel.space_dimension()), (fmm_arc.kernel.space_dimension(), 1)]
}.eval();
Expand Down Expand Up @@ -463,7 +457,10 @@ where
// println!("check potential svd {:?}", check_potential_owned.data());

// Compute local
let locals_owned = (self.fmm.dc2e_inv.dot(&check_potential_owned)
let locals_owned = (self
.fmm
.dc2e_inv_1
.dot(&self.fmm.dc2e_inv_2.dot(&check_potential_owned))
* self.fmm.kernel.scale(level)
* self.m2l_scale(level))
.eval();
Expand Down Expand Up @@ -727,8 +724,12 @@ where
rlst_pointer_mat!['a, f64, check_potentials.as_ptr(), (ncoeffs, ntargets), (1, ncoeffs)]
};

let locals =
(self.fmm.dc2e_inv.dot(&check_potentials) * self.fmm.kernel.scale(level)).eval();
let locals = (self
.fmm
.dc2e_inv_1
.dot(&self.fmm.dc2e_inv_2.dot(&check_potentials))
* self.fmm.kernel.scale(level))
.eval();

for (i, target) in targets.iter().enumerate() {
let target_local_arc = Arc::clone(self.locals.get(target).unwrap());
Expand Down
28 changes: 16 additions & 12 deletions fmm/src/fmm.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,15 +100,17 @@ where
for i in 0..s.len() {
mat_s[[i, i]] = s[i];
}
let uc2e_inv = v.dot(&mat_s).dot(&ut);
let uc2e_inv_1 = v.dot(&mat_s);
let uc2e_inv_2 = ut;

let (s, ut, v) = pinv::<f64>(&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];
}
let dc2e_inv = v.dot(&mat_s).dot(&ut);
let dc2e_inv_1 = v.dot(&mat_s);
let dc2e_inv_2 = ut;

// Calculate M2M/L2L matrices
let children = ROOT.children();
Expand Down Expand Up @@ -139,7 +141,7 @@ where
// Need to transpose so that rows correspond to targets, and columns to sources
let pc2ce = pc2ce.transpose().eval();

m2m.push(uc2e_inv.dot(&pc2ce).eval());
m2m.push(uc2e_inv_1.dot(&uc2e_inv_2.dot(&pc2ce)).eval());

let mut cc2pe = rlst_dynamic_mat![f64, (ncheck_surface, nequiv_surface)];

Expand All @@ -152,13 +154,17 @@ where

// Need to transpose so that rows correspond to targets, and columns to sources
let cc2pe = cc2pe.transpose().eval();
l2l.push((kernel.scale(child.level()) * dc2e_inv.dot(&cc2pe)).eval());
l2l.push(
(kernel.scale(child.level()) * dc2e_inv_1.dot(&dc2e_inv_2.dot(&cc2pe))).eval(),
);
}

Self {
order,
uc2e_inv,
dc2e_inv,
uc2e_inv_1,
uc2e_inv_2,
dc2e_inv_1,
dc2e_inv_2,
alpha_inner,
alpha_outer,
m2m,
Expand Down Expand Up @@ -392,7 +398,7 @@ mod test {
// Setup a FMM experiment
let order = 6;
let alpha_inner = 1.05;
let alpha_outer = 1.99;
let alpha_outer = 2.95;
let adaptive = false;
let k = 1000;
let ncrit = 150;
Expand Down Expand Up @@ -445,7 +451,6 @@ mod test {

let ntargets = leaf_coordinates.len() / datatree.fmm.kernel.space_dimension();

// Get into row major order
let leaf_coordinates = unsafe {
rlst_pointer_mat!['static, f64, leaf_coordinates.as_ptr(), (ntargets, datatree.fmm.kernel.space_dimension()), (datatree.fmm.kernel.space_dimension(), 1)]
}.eval();
Expand Down Expand Up @@ -473,7 +478,7 @@ mod test {
.sum();
let rel_error: f64 = abs_error / (direct.iter().sum::<f64>());

assert!(rel_error <= 1e-5);
assert!(rel_error <= 1e-6);
}

#[test]
Expand All @@ -485,7 +490,7 @@ mod test {

let order = 6;
let alpha_inner = 1.05;
let alpha_outer = 1.95;
let alpha_outer = 2.95;
let adaptive = false;
let ncrit = 150;
let depth = 3;
Expand Down Expand Up @@ -524,7 +529,6 @@ mod test {

let ntargets = leaf_coordinates.len() / datatree.fmm.kernel.space_dimension();

// Get into row major order
let leaf_coordinates = unsafe {
rlst_pointer_mat!['static, f64, leaf_coordinates.as_ptr(), (ntargets, datatree.fmm.kernel.space_dimension()), (datatree.fmm.kernel.space_dimension(), 1)]
}.eval();
Expand All @@ -551,6 +555,6 @@ mod test {
.map(|(a, b)| (a - b).abs())
.sum();
let rel_error: f64 = abs_error / (direct.iter().sum::<f64>());
assert!(rel_error <= 1e-5);
assert!(rel_error <= 1e-6);
}
}
8 changes: 6 additions & 2 deletions fmm/src/types.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,10 +54,14 @@ pub struct KiFmm<T: Tree, U: Kernel, V: FieldTranslationData<U>> {
pub order: usize,

/// The pseudo-inverse of the dense interaction matrix between the upward check and upward equivalent surfaces.
pub uc2e_inv: C2EType,
/// Store in two parts to avoid propagating error from computing pseudo-inverse
pub uc2e_inv_1: C2EType,
pub uc2e_inv_2: C2EType,

/// The pseudo-inverse of the dense interaction matrix between the downward check and downward equivalent surfaces.
pub dc2e_inv: C2EType,
/// Store in two parts to avoid propagating error from computing pseudo-inverse
pub dc2e_inv_1: C2EType,
pub dc2e_inv_2: C2EType,

/// The ratio of the inner check surface diamater in comparison to the surface discretising a box.
pub alpha_inner: f64,
Expand Down
Loading