Merge pull request #9 from wgurecky/col_row_prod

Use faer LinOp and use faer Row Col product in arnoldi

Cargo.toml

@@ -17,4 +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 = "*"

examples/ex_1.rs

@@ -44,12 +44,12 @@ fn main() {
     }
 
     // init guess
-    let x0 = faer::Mat::zeros(216, 1);
+    let mut x0 = faer::Mat::zeros(216, 1);
 
     // solve the system
     let jacobi_pre = JacobiPreconLinOp::new(a_test.as_ref());
-    let (res_x, err, iters) = gmres(a_test.as_ref(), rhs.as_ref(), x0.as_ref(), 500, 1e-8, Some(&jacobi_pre)).unwrap();
-    println!("Result x: {:?}", res_x);
+    let (err, iters) = gmres(a_test.as_ref(), rhs.as_ref(), x0.as_mut(), 500, 1e-8, Some(&jacobi_pre)).unwrap();
+    println!("Result x: {:?}", x0);
     println!("Error x: {:?}", err);
     println!("Iters : {:?}", iters);
 

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::<usize>().unwrap()-1;
+        let idx_j = tmp_vec_line[1].parse::<usize>().unwrap()-1;
+        let val = tmp_vec_line[tmp_vec_line.len()-1].parse::<f64>().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::<f64>().unwrap();
+        b_rhs.push(val);
+    }
+
+    // create sparse mat
+    let a_test = SparseColMat::<usize, f64>::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);
+
+}

readme.md

@@ -2,10 +2,11 @@ About
 =====
 
 [![Crates.io](https://img.shields.io/crates/v/faer_gmres)](https://crates.io/crates/faer_gmres)
+![status](https://github.com/wgurecky/faer-gmres/actions/workflows/rust.yml/badge.svg)
 
 GMRES in rust using [faer](https://github.com/sarah-ek/faer-rs).
 
-Solves linear systems of the form: Ax=b, where A is sparse.  Depends on faer for sparse matrix implementation and sparse QR solver.
+Solves linear systems of the form: Ax=b, where A is sparse.  Depends on faer for sparse matrix implementation.
 
 Use
 ===
@@ -35,15 +36,16 @@ Example use:
         ];
 
     // init sol guess
-    let x0 = faer::mat![
+    // Note: x is modified in-place, the result is stored in x
+    let mut x = faer::mat![
         [0.0],
         [0.0],
         [0.0],
         ];
 
     // the final None arg means do not apply left preconditioning
-    let (res_x, err, iters) = gmres(a_test.as_ref(), b.as_ref(), x0.as_ref(), 10, 1e-8, None).unwrap();
-    println!("Result x: {:?}", res_x);
+    let (err, iters) = gmres(a_test.as_ref(), b.as_ref(), x.as_mut(), 10, 1e-8, None).unwrap();
+    println!("Result x: {:?}", x);
     println!("Error x: {:?}", err);
     println!("Iters : {:?}", iters);
 
@@ -52,9 +54,9 @@ 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();
+    let (err, iters) = gmres(a_test.as_ref(), b.as_ref(), x.as_mut(), 10, 1e-8, Some(&jacobi_pre)).unwrap();
 
 ## Restarted GMRES:
 
@@ -63,8 +65,8 @@ A restarted GMRES routine is provided:
     use faer_gmres::restarted_gmres;
     let max_inner = 30;
     let max_outer = 50;
-    let (res_x, err, iters) = restarted_gmres(
-        a_test.as_ref(), b.as_ref(), x0.as_ref(), max_inner, max_outer, 1e-8, None).unwrap();
+    let (err, iters) = restarted_gmres(
+        a_test.as_ref(), b.as_ref(), x.as_mut(), max_inner, max_outer, 1e-8, None).unwrap();
 
 This will repeatedly call the inner GMRES routine, using the previous outer iteration's solution as the inital guess for the next outer solve.  The current implementation of restarted GMRES in this package can reduce the memory requirements needed, but slow convergence.