From 23b3a10c1c7251198cad5a8bbadc8f048752726b Mon Sep 17 00:00:00 2001 From: William Gurecky Date: Wed, 7 Aug 2024 23:31:02 -0400 Subject: [PATCH] migrate to faer LinOp for precond --- Cargo.toml | 2 +- examples/ex_2.rs | 60 ++++++++++++++++++++++++++++++++++ readme.md | 2 +- src/lib.rs | 85 +++++++++++++++++++++++++++++++++++++++--------- 4 files changed, 132 insertions(+), 17 deletions(-) create mode 100644 examples/ex_2.rs diff --git a/Cargo.toml b/Cargo.toml index 04f25cb..448c2be 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,5 +17,5 @@ path = 'src/lib.rs' [dependencies] assert_approx_eq = "1.1.0" num-traits = "0.2.18" -faer = {version = "0.19.0"} +faer = {version = "0.19.1", features=["unstable"]} reborrow = "*" diff --git a/examples/ex_2.rs b/examples/ex_2.rs new file mode 100644 index 0000000..8b0f799 --- /dev/null +++ b/examples/ex_2.rs @@ -0,0 +1,60 @@ +use faer_gmres::{gmres, restarted_gmres, JacobiPreconLinOp}; +use faer::prelude::*; +use faer::sparse::*; +use std::fs::read_to_string; +use std::time::Instant; + + +fn main() { + let mut a_triplets = vec![]; + let mut b_rhs = vec![]; + + // read A matrix from file + for line in read_to_string("./examples/data/e40r0100.txt").unwrap().lines() { + let mut tmp_vec_line = vec![]; + let iter = line.split_whitespace(); + for word in iter { + tmp_vec_line.push(word); + } + let idx_i = tmp_vec_line[0].parse::().unwrap()-1; + let idx_j = tmp_vec_line[1].parse::().unwrap()-1; + let val = tmp_vec_line[tmp_vec_line.len()-1].parse::().unwrap(); + a_triplets.push((idx_i, idx_j, val)); + } + + // read rhs from file + for line in read_to_string("./examples/data/e40r0100_rhs1.txt").unwrap().lines() { + let mut tmp_vec_line = vec![]; + let iter = line.split_whitespace(); + for word in iter { + tmp_vec_line.push(word); + } + let val = tmp_vec_line[tmp_vec_line.len()-1].parse::().unwrap(); + b_rhs.push(val); + } + + // create sparse mat + let a_test = SparseColMat::::try_new_from_triplets( + 17281, 17281, + &a_triplets).unwrap(); + + // create rhs + let mut rhs = faer::Mat::zeros(17281, 1); + for (i, rhs_val) in b_rhs.into_iter().enumerate() { + rhs.write(i, 0, rhs_val); + } + + // init guess + let mut x0 = faer::Mat::zeros(17281, 1); + + // solve the system + let jacobi_pre = JacobiPreconLinOp::new(a_test.as_ref()); + let now = Instant::now(); + let (err, iters) = gmres(a_test.as_ref(), rhs.as_ref(), x0.as_mut(), 5000, 1e-6, Some(&jacobi_pre)).unwrap(); + let dt = now.elapsed(); + println!("Result x: {:?}", x0); + println!("Error x: {:?}", err); + println!("Iters : {:?}", iters); + println!("Solve time: {:?} ", dt); + +} diff --git a/readme.md b/readme.md index 5c725fe..d6e4d8a 100644 --- a/readme.md +++ b/readme.md @@ -53,7 +53,7 @@ Example use: A preconditioner can be supplied: // continued from above... - use faer_gmres::{JacobiPreconLinOp, LinOp}; + use faer_gmres::{JacobiPreconLinOp}; let jacobi_pre = JacobiPreconLinOp::new(a_test.as_ref()); let (res_x, err, iters) = gmres(a_test.as_ref(), b.as_ref(), x0.as_ref(), 10, 1e-8, Some(&jacobi_pre)).unwrap(); diff --git a/src/lib.rs b/src/lib.rs index 110af7a..018c499 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -34,6 +34,9 @@ use faer::prelude::*; use faer::sparse::*; use faer::mat; +use faer::linop::LinOp; +use faer::Parallelism; +use faer::dyn_stack::PodStack; use reborrow::*; use faer::utils::simd::Write; use num_traits::Float; @@ -63,36 +66,78 @@ impl fmt::Display for GmresError } } -pub trait LinOp - where - T: faer::RealField + Float -{ - fn apply_linop_to_vec(&self, target: MatMut); -} +// pub trait LinOp +// where +// T: faer::RealField + Float +// { +// fn apply_linop_to_vec(&self, target: MatRef) -> Mat; +// } -#[derive(Clone)] +#[derive(Clone,Debug)] pub struct JacobiPreconLinOp<'a, T> where T: faer::RealField + Float { m: SparseColMatRef<'a, usize, T>, } + impl <'a, T> LinOp for JacobiPreconLinOp<'a, T> where T: faer::RealField + Float + faer::SimpleEntity { - fn apply_linop_to_vec(&self, mut target: MatMut) { + fn apply_req( + &self, + rhs_ncols: usize, + parallelism: Parallelism, + ) -> Result { + let _ = parallelism; + let _ = rhs_ncols; + Ok(faer::dyn_stack::StackReq::empty()) + } + + fn nrows(&self) -> usize { + self.m.nrows() + } + + fn ncols(&self) -> usize { + self.m.ncols() + } + + fn apply( + &self, + out: MatMut, + rhs: MatRef, + parallelism: Parallelism, + stack: PodStack<'_>, + ) + { + // unused + _ = parallelism; + _ = stack; + + let mut out = out; let eps = T::from(1e-12).unwrap(); let one_c = T::from(1.0).unwrap(); - let zero_c = T::from(0.0).unwrap(); - for i in 0..target.nrows() + for i in 0..rhs.nrows() { - let v = target.read(i, 0); - target.write(i, 0, + let v = rhs.read(i, 0); + out.write(i, 0, v * (one_c / (*self.m.as_ref().get(i, i).unwrap_or(&one_c) + eps) )); } } + + fn conj_apply( + &self, + out: MatMut<'_, T>, + rhs: MatRef<'_, T>, + parallelism: Parallelism, + stack: PodStack<'_>, + ) { + // Not implented error! + panic!("Not Implemented"); + } } + impl <'a, T> JacobiPreconLinOp <'a, T> where T: faer::RealField + Float @@ -104,7 +149,6 @@ impl <'a, T> JacobiPreconLinOp <'a, T> } } - /// Calculate the givens rotation matrix fn givens_rotation(v1: T, v2: T) -> (T, T) where @@ -165,7 +209,12 @@ fn arnoldi<'a, T>( // Apply left preconditioner if supplied match m { - Some(m) => m.apply_linop_to_vec(qv.as_mut()), + Some(m) => { + let mut _dummy_podstack: [u8;1] = [0u8;1]; + let mut lp_out = faer::Mat::zeros(qv.nrows(), qv.ncols()); + m.apply(lp_out.as_mut(), qv.as_ref(), faer::get_global_parallelism(), PodStack::new(&mut _dummy_podstack)); + qv = lp_out; + }, _ => {} } @@ -197,8 +246,14 @@ pub fn gmres<'a, T>( { // compute initial residual let mut r = b - a * x.as_ref(); + match &m { - Some(m) => (&m).apply_linop_to_vec(r.as_mut()), + Some(m) => { + let mut _dummy_podstack: [u8;1] = [0u8;1]; + let mut lp_out = faer::Mat::zeros(r.nrows(), r.ncols()); + (&m).apply(lp_out.as_mut(), r.as_ref(), faer::get_global_parallelism(), PodStack::new(&mut _dummy_podstack)); + r = lp_out; + }, _ => {} }