[Important] Remove Generic SparseMatrix and make complex sparse solvers take SparseCooMatrix instead. Thus, eliminating the memoization problem

Author: cpmech

russell_ode/src/radau5.rs

@@ -2,8 +2,9 @@ use crate::StrError;
 use crate::{OdeSolverTrait, Params, System, Workspace};
 use russell_lab::math::SQRT_6;
 use russell_lab::{complex_vec_zip, cpx, format_fortran, vec_copy, Complex64, ComplexVector, Vector};
-use russell_sparse::{numerical_jacobian, ComplexCscMatrix, CooMatrix, CscMatrix};
-use russell_sparse::{ComplexLinSolver, ComplexSparseMatrix, Genie, LinSolver};
+use russell_sparse::{
+    numerical_jacobian, ComplexCooMatrix, ComplexCscMatrix, ComplexLinSolver, CooMatrix, CscMatrix, Genie, LinSolver,
+};
 use std::thread;
 
 /// Implements the Radau5 method (Radau IIA) (implicit, order 5, embedded) for ODEs and DAEs
@@ -41,7 +42,7 @@ pub(crate) struct Radau5<'a, A> {
     kk_real: CooMatrix,
 
     /// Coefficient matrix (for real system). K_comp = (α + βi) M - J
-    kk_comp: ComplexSparseMatrix,
+    kk_comp: ComplexCooMatrix,
 
     /// Linear solver (for real system)
     solver_real: LinSolver<'a>,
@@ -148,7 +149,7 @@ impl<'a, A> Radau5<'a, A> {
             mass,
             jj: CooMatrix::new(ndim, ndim, jac_nnz, sym).unwrap(),
             kk_real: CooMatrix::new(ndim, ndim, nnz, sym).unwrap(),
-            kk_comp: ComplexSparseMatrix::new_coo(ndim, ndim, nnz, sym).unwrap(),
+            kk_comp: ComplexCooMatrix::new(ndim, ndim, nnz, sym).unwrap(),
             solver_real: LinSolver::new(params.newton.genie).unwrap(),
             solver_comp: ComplexLinSolver::new(params.newton.genie).unwrap(),
             reuse_jacobian: false,
@@ -195,7 +196,7 @@ impl<'a, A> Radau5<'a, A> {
         // auxiliary
         let jj = &mut self.jj; // J = df/dy
         let kk_real = &mut self.kk_real; // K_real = γ M - J
-        let kk_comp = self.kk_comp.get_coo_mut().unwrap(); // K_comp = (α + βi) M - J
+        let kk_comp = &mut self.kk_comp; // K_comp = (α + βi) M - J
 
         // Jacobian matrix
         if self.reuse_jacobian {
@@ -261,7 +262,7 @@ impl<'a, A> Radau5<'a, A> {
             .factorize(&self.kk_real, self.params.newton.lin_sol_params)?;
         self.solver_comp
             .actual
-            .factorize(&mut self.kk_comp, self.params.newton.lin_sol_params)
+            .factorize(&self.kk_comp, self.params.newton.lin_sol_params)
     }
 
     /// Factorizes the real and complex systems concurrently
@@ -276,7 +277,7 @@ impl<'a, A> Radau5<'a, A> {
             let handle_comp = scope.spawn(|| {
                 self.solver_comp
                     .actual
-                    .factorize(&mut self.kk_comp, self.params.newton.lin_sol_params)
+                    .factorize(&self.kk_comp, self.params.newton.lin_sol_params)
                     .unwrap();
             });
             let err_real = handle_real.join();
@@ -296,9 +297,7 @@ impl<'a, A> Radau5<'a, A> {
     /// Solves the real and complex linear systems
     fn solve_lin_sys(&mut self) -> Result<(), StrError> {
         self.solver_real.actual.solve(&mut self.dw0, &self.v0, false)?;
-        self.solver_comp
-            .actual
-            .solve(&mut self.dw12, &self.kk_comp, &self.v12, false)?;
+        self.solver_comp.actual.solve(&mut self.dw12, &self.v12, false)?;
         Ok(())
     }
 
@@ -309,10 +308,7 @@ impl<'a, A> Radau5<'a, A> {
                 self.solver_real.actual.solve(&mut self.dw0, &self.v0, false).unwrap();
             });
             let handle_comp = scope.spawn(|| {
-                self.solver_comp
-                    .actual
-                    .solve(&mut self.dw12, &self.kk_comp, &self.v12, false)
-                    .unwrap();
+                self.solver_comp.actual.solve(&mut self.dw12, &self.v12, false).unwrap();
             });
             let err_real = handle_real.join();
             let err_comp = handle_comp.join();

russell_sparse/examples/complex_system.rs

@@ -42,7 +42,7 @@ fn solve(genie: Genie) -> Result<(), StrError> {
     let nnz = 16; // number of non-zero values, including duplicates
 
     // input matrix in Complex Triplet format
-    let mut coo = ComplexSparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?;
+    let mut coo = ComplexCooMatrix::new(ndim, ndim, nnz, Sym::No)?;
 
     // first column
     coo.put(0, 0, cpx!(19.73, 0.00))?;
@@ -88,8 +88,8 @@ fn solve(genie: Genie) -> Result<(), StrError> {
 
     // call factorize and solve
     let mut x = ComplexVector::new(ndim);
-    solver.actual.factorize(&mut coo, Some(params))?;
-    solver.actual.solve(&mut x, &coo, &b, false)?;
+    solver.actual.factorize(&coo, Some(params))?;
+    solver.actual.solve(&mut x, &b, false)?;
     println!("x =\n{}", x);
 
     // check

russell_sparse/src/aliases.rs

@@ -1,4 +1,4 @@
-use crate::{NumCooMatrix, NumCscMatrix, NumCsrMatrix, NumSparseMatrix};
+use crate::{NumCooMatrix, NumCscMatrix, NumCsrMatrix};
 use russell_lab::Complex64;
 
 /// Defines an alias to NumCooMatrix with f64
@@ -10,9 +10,6 @@ pub type CscMatrix = NumCscMatrix<f64>;
 /// Defines an alias to NumCsrMatrix with f64
 pub type CsrMatrix = NumCsrMatrix<f64>;
 
-/// Defines an alias to NumSparseMatrix with f64
-pub type SparseMatrix = NumSparseMatrix<f64>;
-
 /// Defines an alias to NumCooMatrix with Complex64
 pub type ComplexCooMatrix = NumCooMatrix<Complex64>;
 
@@ -21,6 +18,3 @@ pub type ComplexCscMatrix = NumCscMatrix<Complex64>;
 
 /// Defines an alias to NumCsrMatrix with Complex64
 pub type ComplexCsrMatrix = NumCsrMatrix<Complex64>;
-
-/// Defines an alias to NumSparseMatrix with Complex64
-pub type ComplexSparseMatrix = NumSparseMatrix<Complex64>;

russell_sparse/src/bin/mem_check.rs

@@ -73,9 +73,8 @@ fn test_complex_solver(genie: Genie) {
         Genie::Mumps => Samples::complex_symmetric_3x3_lower().0,
         Genie::Umfpack => Samples::complex_symmetric_3x3_full().0,
     };
-    let mut mat = ComplexSparseMatrix::from_coo(coo);
 
-    match solver.actual.factorize(&mut mat, None) {
+    match solver.actual.factorize(&coo, None) {
         Err(e) => {
             println!("FAIL(factorize): {}", e);
             return;
@@ -86,15 +85,15 @@ fn test_complex_solver(genie: Genie) {
     let mut x = ComplexVector::new(3);
     let rhs = ComplexVector::from(&[cpx!(-3.0, 3.0), cpx!(2.0, -2.0), cpx!(9.0, 7.0)]);
 
-    match solver.actual.solve(&mut x, &mut mat, &rhs, false) {
+    match solver.actual.solve(&mut x, &rhs, false) {
         Err(e) => {
             println!("FAIL(solve): {}", e);
             return;
         }
         _ => (),
     }
 
-    match solver.actual.solve(&mut x, &mat, &rhs, false) {
+    match solver.actual.solve(&mut x, &rhs, false) {
         Err(e) => {
             println!("FAIL(solve again): {}", e);
             return;

russell_sparse/src/bin/solve_matrix_market.rs

@@ -173,11 +173,8 @@ fn main() -> Result<(), StrError> {
             csc.write_matrix_market("/tmp/russell_sparse/solve_matrix_market_complex.smat", true)?;
         }
 
-        // save the COO matrix as a generic SparseMatrix
-        let mut mat = ComplexSparseMatrix::from_coo(coo);
-
         // save information about the matrix
-        let (nrow, ncol, nnz, sym) = mat.get_info();
+        let (nrow, ncol, nnz, sym) = coo.get_info();
         stats.set_matrix_name_from_path(&opt.matrix_market_file);
         stats.matrix.nrow = nrow;
         stats.matrix.ncol = ncol;
@@ -189,18 +186,18 @@ fn main() -> Result<(), StrError> {
         let mut solver = ComplexLinSolver::new(genie)?;
 
         // call factorize
-        solver.actual.factorize(&mut mat, Some(params))?;
+        solver.actual.factorize(&coo, Some(params))?;
 
         // allocate vectors
         let mut x = ComplexVector::new(nrow);
         let rhs = ComplexVector::filled(nrow, cpx!(1.0, 1.0));
 
         // solve linear system
-        solver.actual.solve(&mut x, &mat, &rhs, opt.verbose)?;
+        solver.actual.solve(&mut x, &rhs, opt.verbose)?;
 
         // verify the solution
         sw.reset();
-        stats.verify = VerifyLinSys::from_complex(&mat, &x, &rhs)?;
+        stats.verify = VerifyLinSys::from_complex(&coo, &x, &rhs)?;
         stats.time_nanoseconds.verify = sw.stop();
 
         // update stats

russell_sparse/src/complex_lin_solver.rs

@@ -1,9 +1,8 @@
 #[cfg(feature = "with_mumps")]
 use super::ComplexSolverMUMPS;
 
-use super::{ComplexSparseMatrix, Genie, LinSolParams, StatsLinSol};
+use super::{ComplexCooMatrix, ComplexSolverKLU, ComplexSolverUMFPACK, Genie, LinSolParams, StatsLinSol};
 use crate::StrError;
-use crate::{ComplexSolverKLU, ComplexSolverUMFPACK};
 use russell_lab::ComplexVector;
 
 /// Defines a unified interface for complex linear system solvers
@@ -24,7 +23,7 @@ pub trait ComplexLinSolTrait: Send {
     ///    kept the same for the next calls.
     /// 3. If the structure of the matrix needs to be changed, the solver must
     ///    be "dropped" and a new solver allocated.
-    fn factorize(&mut self, mat: &mut ComplexSparseMatrix, params: Option<LinSolParams>) -> Result<(), StrError>;
+    fn factorize(&mut self, mat: &ComplexCooMatrix, params: Option<LinSolParams>) -> Result<(), StrError>;
 
     /// Computes the solution of the linear system
     ///
@@ -46,13 +45,7 @@ pub trait ComplexLinSolTrait: Send {
     /// * `verbose` -- shows messages
     ///
     /// **Warning:** the matrix must be same one used in `factorize`.
-    fn solve(
-        &mut self,
-        x: &mut ComplexVector,
-        mat: &ComplexSparseMatrix,
-        rhs: &ComplexVector,
-        verbose: bool,
-    ) -> Result<(), StrError>;
+    fn solve(&mut self, x: &mut ComplexVector, rhs: &ComplexVector, verbose: bool) -> Result<(), StrError>;
 
     /// Updates the stats structure (should be called after solve)
     fn update_stats(&self, stats: &mut StatsLinSol);
@@ -125,14 +118,14 @@ impl<'a> ComplexLinSolver<'a> {
     pub fn compute(
         genie: Genie,
         x: &mut ComplexVector,
-        mat: &mut ComplexSparseMatrix,
+        mat: &ComplexCooMatrix,
         rhs: &ComplexVector,
         params: Option<LinSolParams>,
     ) -> Result<Self, StrError> {
         let mut solver = ComplexLinSolver::new(genie)?;
         solver.actual.factorize(mat, params)?;
         let verbose = if let Some(p) = params { p.verbose } else { false };
-        solver.actual.solve(x, mat, rhs, verbose)?;
+        solver.actual.solve(x, rhs, verbose)?;
         Ok(solver)
     }
 }
@@ -142,7 +135,7 @@ impl<'a> ComplexLinSolver<'a> {
 #[cfg(test)]
 mod tests {
     use super::ComplexLinSolver;
-    use crate::{ComplexSparseMatrix, Genie, Samples};
+    use crate::{Genie, Samples};
     use russell_lab::{complex_vec_approx_eq, cpx, Complex64, ComplexVector};
 
     #[cfg(feature = "with_mumps")]
@@ -151,10 +144,9 @@ mod tests {
     #[test]
     fn complex_lin_solver_compute_works_klu() {
         let (coo, _, _, _) = Samples::complex_symmetric_3x3_full();
-        let mut mat = ComplexSparseMatrix::from_coo(coo);
         let mut x = ComplexVector::new(3);
         let rhs = ComplexVector::from(&[cpx!(-3.0, 3.0), cpx!(2.0, -2.0), cpx!(9.0, 7.0)]);
-        ComplexLinSolver::compute(Genie::Klu, &mut x, &mut mat, &rhs, None).unwrap();
+        ComplexLinSolver::compute(Genie::Klu, &mut x, &coo, &rhs, None).unwrap();
         let x_correct = &[cpx!(1.0, 1.0), cpx!(2.0, -2.0), cpx!(3.0, 3.0)];
         complex_vec_approx_eq(&x, x_correct, 1e-15);
     }
@@ -164,21 +156,19 @@ mod tests {
     #[cfg(feature = "with_mumps")]
     fn complex_lin_solver_compute_works_mumps() {
         let (coo, _, _, _) = Samples::complex_symmetric_3x3_lower();
-        let mut mat = ComplexSparseMatrix::from_coo(coo);
         let mut x = ComplexVector::new(3);
         let rhs = ComplexVector::from(&[cpx!(-3.0, 3.0), cpx!(2.0, -2.0), cpx!(9.0, 7.0)]);
-        ComplexLinSolver::compute(Genie::Mumps, &mut x, &mut mat, &rhs, None).unwrap();
+        ComplexLinSolver::compute(Genie::Mumps, &mut x, &coo, &rhs, None).unwrap();
         let x_correct = &[cpx!(1.0, 1.0), cpx!(2.0, -2.0), cpx!(3.0, 3.0)];
         complex_vec_approx_eq(&x, x_correct, 1e-15);
     }
 
     #[test]
     fn complex_lin_solver_compute_works_umfpack() {
         let (coo, _, _, _) = Samples::complex_symmetric_3x3_full();
-        let mut mat = ComplexSparseMatrix::from_coo(coo);
         let mut x = ComplexVector::new(3);
         let rhs = ComplexVector::from(&[cpx!(-3.0, 3.0), cpx!(2.0, -2.0), cpx!(9.0, 7.0)]);
-        ComplexLinSolver::compute(Genie::Umfpack, &mut x, &mut mat, &rhs, None).unwrap();
+        ComplexLinSolver::compute(Genie::Umfpack, &mut x, &coo, &rhs, None).unwrap();
         let x_correct = &[cpx!(1.0, 1.0), cpx!(2.0, -2.0), cpx!(3.0, 3.0)];
         complex_vec_approx_eq(&x, x_correct, 1e-15);
     }