From b9a8254e9ca99c00c2a635b36c8a8f9022025221 Mon Sep 17 00:00:00 2001 From: Dorival Pedroso Date: Sat, 9 Mar 2024 11:42:35 +1000 Subject: [PATCH 1/3] [sparse] Rename Symmetry to Sym --- russell_ode/src/samples.rs | 8 +- russell_ode/src/system.rs | 10 +- russell_sparse/src/complex_coo_matrix.rs | 6 +- russell_sparse/src/complex_solver_mumps.rs | 16 +-- russell_sparse/src/complex_solver_umfpack.rs | 10 +- russell_sparse/src/coo_matrix.rs | 54 ++++----- russell_sparse/src/csc_matrix.rs | 24 ++-- russell_sparse/src/csr_matrix.rs | 24 ++-- russell_sparse/src/enums.rs | 108 +++++++++--------- russell_sparse/src/read_matrix_market.rs | 18 +-- russell_sparse/src/samples.rs | 20 ++-- russell_sparse/src/solver_mumps.rs | 12 +- russell_sparse/src/solver_umfpack.rs | 12 +- russell_sparse/src/sparse_matrix.rs | 40 +++---- russell_sparse/src/write_matrix_market.rs | 6 +- .../tests/test_complex_coo_matrix.rs | 2 +- 16 files changed, 182 insertions(+), 188 deletions(-) diff --git a/russell_ode/src/samples.rs b/russell_ode/src/samples.rs index 085cb23f..09344475 100644 --- a/russell_ode/src/samples.rs +++ b/russell_ode/src/samples.rs @@ -2,7 +2,7 @@ use crate::StrError; use crate::{HasJacobian, NoArgs, System}; use russell_lab::math::PI; use russell_lab::Vector; -use russell_sparse::{CooMatrix, Symmetry}; +use russell_sparse::{CooMatrix, Sym}; /// Holds data corresponding to a sample ODE problem pub struct SampleData { @@ -164,7 +164,7 @@ impl Samples { ) { // selected symmetry option (for both Mass and Jacobian matrices) let symmetry = if lower_triangle { - Some(Symmetry::new_general_lower()) + Some(Sym::new_general_lower()) } else { None }; @@ -789,7 +789,7 @@ mod tests { use super::{NoArgs, Samples}; use crate::StrError; use russell_lab::{deriv_central5, mat_approx_eq, vec_approx_eq, Matrix, Vector}; - use russell_sparse::{CooMatrix, Symmetry}; + use russell_sparse::{CooMatrix, Sym}; fn numerical_jacobian(ndim: usize, x0: f64, y0: Vector, function: F, multiplier: f64) -> Matrix where @@ -1045,6 +1045,6 @@ mod tests { println!("{}", mass.as_dense()); let ndim = system.ndim; let nnz_mass = 5 + 4; - assert_eq!(mass.get_info(), (ndim, ndim, nnz_mass, Symmetry::No)); + assert_eq!(mass.get_info(), (ndim, ndim, nnz_mass, Sym::No)); } } diff --git a/russell_ode/src/system.rs b/russell_ode/src/system.rs index 2edd5ed4..078e8cbd 100644 --- a/russell_ode/src/system.rs +++ b/russell_ode/src/system.rs @@ -1,6 +1,6 @@ use crate::{HasJacobian, StrError}; use russell_lab::Vector; -use russell_sparse::{CooMatrix, Symmetry}; +use russell_sparse::{CooMatrix, Sym}; use std::marker::PhantomData; /// Indicates that the system functions do not require extra arguments @@ -72,7 +72,7 @@ where pub(crate) jac_nnz: usize, /// Symmetry properties of the Jacobian matrix - pub(crate) jac_symmetry: Symmetry, + pub(crate) jac_symmetry: Sym, /// Holds the mass matrix pub(crate) mass_matrix: Option, @@ -106,7 +106,7 @@ where jacobian: J, has_ana_jacobian: HasJacobian, jac_nnz: Option, - jac_symmetry: Option, + jac_symmetry: Option, ) -> Self { let jac_available = match has_ana_jacobian { HasJacobian::Yes => true, @@ -118,7 +118,7 @@ where jacobian, jac_available, jac_nnz: if let Some(n) = jac_nnz { n } else { ndim * ndim }, - jac_symmetry: if let Some(s) = jac_symmetry { s } else { Symmetry::No }, + jac_symmetry: if let Some(s) = jac_symmetry { s } else { Sym::No }, mass_matrix: None, phantom: PhantomData, } @@ -132,7 +132,7 @@ where /// /// * `max_nnz` -- Max number of non-zero values pub fn init_mass_matrix(&mut self, max_nnz: usize) -> Result<(), StrError> { - let sym = if self.jac_symmetry == Symmetry::No { + let sym = if self.jac_symmetry == Sym::No { None } else { Some(self.jac_symmetry) diff --git a/russell_sparse/src/complex_coo_matrix.rs b/russell_sparse/src/complex_coo_matrix.rs index 629ae12a..34ea0482 100644 --- a/russell_sparse/src/complex_coo_matrix.rs +++ b/russell_sparse/src/complex_coo_matrix.rs @@ -84,13 +84,13 @@ impl ComplexCooMatrix { #[cfg(test)] mod tests { - use crate::{ComplexCooMatrix, CooMatrix, Storage, Symmetry}; + use crate::{ComplexCooMatrix, CooMatrix, Storage, Sym}; use num_complex::Complex64; use russell_lab::cpx; #[test] fn assign_capture_errors() { - let sym = Some(Symmetry::General(Storage::Full)); + let sym = Some(Sym::General(Storage::Full)); let nnz_a = 1; let nnz_b = 2; // wrong: must be ≤ nnz_a let mut a_1x2 = ComplexCooMatrix::new(1, 2, nnz_a, None).unwrap(); @@ -148,7 +148,7 @@ mod tests { #[test] fn augment_capture_errors() { - let sym = Some(Symmetry::General(Storage::Full)); + let sym = Some(Sym::General(Storage::Full)); let nnz_a = 1; let nnz_b = 1; let mut a_1x2 = ComplexCooMatrix::new(1, 2, nnz_a /* + nnz_b */, None).unwrap(); diff --git a/russell_sparse/src/complex_solver_mumps.rs b/russell_sparse/src/complex_solver_mumps.rs index dd8df0eb..229e3d41 100644 --- a/russell_sparse/src/complex_solver_mumps.rs +++ b/russell_sparse/src/complex_solver_mumps.rs @@ -1,5 +1,5 @@ use super::{handle_mumps_error_code, mumps_ordering, mumps_scaling}; -use super::{ComplexLinSolTrait, ComplexSparseMatrix, LinSolParams, StatsLinSol, Symmetry}; +use super::{ComplexLinSolTrait, ComplexSparseMatrix, LinSolParams, StatsLinSol, Sym}; use super::{ MUMPS_ORDERING_AMD, MUMPS_ORDERING_AMF, MUMPS_ORDERING_AUTO, MUMPS_ORDERING_METIS, MUMPS_ORDERING_PORD, MUMPS_ORDERING_QAMD, MUMPS_ORDERING_SCOTCH, MUMPS_SCALING_AUTO, MUMPS_SCALING_COLUMN, MUMPS_SCALING_DIAGONAL, @@ -81,7 +81,7 @@ pub struct ComplexSolverMUMPS { factorized: bool, /// Holds the symmetry type used in the initialize - initialized_symmetry: Symmetry, + initialized_symmetry: Sym, /// Holds the matrix dimension saved in initialize initialized_ndim: usize, @@ -161,7 +161,7 @@ impl ComplexSolverMUMPS { solver, initialized: false, factorized: false, - initialized_symmetry: Symmetry::No, + initialized_symmetry: Sym::No, initialized_ndim: 0, initialized_nnz: 0, effective_ordering: -1, @@ -455,7 +455,7 @@ impl ComplexLinSolTrait for ComplexSolverMUMPS { #[cfg(test)] mod tests { use super::*; - use crate::{ComplexCooMatrix, ComplexSparseMatrix, LinSolParams, Ordering, Samples, Scaling, Storage, Symmetry}; + use crate::{ComplexCooMatrix, ComplexSparseMatrix, LinSolParams, Ordering, Samples, Scaling, Storage, Sym}; use num_complex::Complex64; use russell_lab::{complex_approx_eq, complex_vec_approx_eq, cpx, ComplexVector}; use serial_test::serial; @@ -492,7 +492,7 @@ mod tests { solver.factorize(&mut mat, None).err(), Some("the COO matrix must have at least one non-zero value") ); - let mut coo = ComplexCooMatrix::new(1, 1, 1, Some(Symmetry::General(Storage::Full))).unwrap(); + let mut coo = ComplexCooMatrix::new(1, 1, 1, Some(Sym::General(Storage::Full))).unwrap(); coo.put(0, 0, cpx!(4.0, 4.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( @@ -508,7 +508,7 @@ mod tests { // ... factorize once => OK solver.factorize(&mut mat, None).unwrap(); // ... change matrix (symmetry) - let mut coo = ComplexCooMatrix::new(2, 2, 2, Some(Symmetry::General(Storage::Full))).unwrap(); + let mut coo = ComplexCooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); coo.put(1, 1, cpx!(2.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); @@ -576,7 +576,7 @@ mod tests { ); // wrong symmetry let rhs = ComplexVector::new(2); - let mut coo_wrong = ComplexCooMatrix::new(2, 2, 2, Some(Symmetry::General(Storage::Full))).unwrap(); + let mut coo_wrong = ComplexCooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); coo_wrong.put(0, 0, cpx!(123.0, 1.0)).unwrap(); coo_wrong.put(1, 1, cpx!(456.0, 2.0)).unwrap(); let mut mat_wrong = ComplexSparseMatrix::from_coo(coo_wrong); @@ -658,7 +658,7 @@ mod tests { complex_vec_approx_eq(x_again.as_data(), x_correct, 1e-14); // solve with positive-definite matrix works - let sym = Some(Symmetry::PositiveDefinite(Storage::Lower)); + let sym = Some(Sym::PositiveDefinite(Storage::Lower)); let nrow = 5; let ncol = 5; let mut coo_pd_lower = ComplexCooMatrix::new(nrow, ncol, 9, sym).unwrap(); diff --git a/russell_sparse/src/complex_solver_umfpack.rs b/russell_sparse/src/complex_solver_umfpack.rs index 535b8d49..1730b737 100644 --- a/russell_sparse/src/complex_solver_umfpack.rs +++ b/russell_sparse/src/complex_solver_umfpack.rs @@ -1,5 +1,5 @@ use super::{handle_umfpack_error_code, umfpack_ordering, umfpack_scaling}; -use super::{ComplexLinSolTrait, ComplexSparseMatrix, LinSolParams, StatsLinSol, Symmetry}; +use super::{ComplexLinSolTrait, ComplexSparseMatrix, LinSolParams, StatsLinSol, Sym}; use super::{ UMFPACK_ORDERING_AMD, UMFPACK_ORDERING_BEST, UMFPACK_ORDERING_CHOLMOD, UMFPACK_ORDERING_METIS, UMFPACK_ORDERING_NONE, UMFPACK_SCALE_MAX, UMFPACK_SCALE_NONE, UMFPACK_SCALE_SUM, UMFPACK_STRATEGY_AUTO, @@ -83,7 +83,7 @@ pub struct ComplexSolverUMFPACK { factorized: bool, /// Holds the symmetry type used in initialize - initialized_symmetry: Symmetry, + initialized_symmetry: Sym, /// Holds the matrix dimension saved in initialize initialized_ndim: usize, @@ -152,7 +152,7 @@ impl ComplexSolverUMFPACK { solver, initialized: false, factorized: false, - initialized_symmetry: Symmetry::No, + initialized_symmetry: Sym::No, initialized_ndim: 0, initialized_nnz: 0, effective_strategy: -1, @@ -456,7 +456,7 @@ mod tests { // ... factorize once => OK solver.factorize(&mut mat, None).unwrap(); // ... change matrix (symmetry) - let mut coo = ComplexCooMatrix::new(2, 2, 2, Some(Symmetry::General(Storage::Full))).unwrap(); + let mut coo = ComplexCooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); coo.put(1, 1, cpx!(2.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); @@ -549,7 +549,7 @@ mod tests { ); // wrong symmetry let rhs = ComplexVector::new(2); - let mut coo_wrong = ComplexCooMatrix::new(2, 2, 2, Some(Symmetry::General(Storage::Full))).unwrap(); + let mut coo_wrong = ComplexCooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); coo_wrong.put(0, 0, cpx!(123.0, 1.0)).unwrap(); coo_wrong.put(1, 1, cpx!(456.0, 2.0)).unwrap(); let mut mat_wrong = ComplexSparseMatrix::from_coo(coo_wrong); diff --git a/russell_sparse/src/coo_matrix.rs b/russell_sparse/src/coo_matrix.rs index 1fc6c702..bf352cf6 100644 --- a/russell_sparse/src/coo_matrix.rs +++ b/russell_sparse/src/coo_matrix.rs @@ -1,4 +1,4 @@ -use super::{Storage, Symmetry}; +use super::{Storage, Sym}; use crate::to_i32; use crate::StrError; use num_traits::{Num, NumCast}; @@ -23,7 +23,7 @@ where T: AddAssign + MulAssign + Num + NumCast + Copy + DeserializeOwned + Serialize, { /// Defines the symmetry and storage: lower-triangular, upper-triangular, full-matrix - pub(crate) symmetry: Symmetry, + pub(crate) symmetry: Sym, /// Holds the number of rows (must fit i32) pub(crate) nrow: usize, @@ -167,7 +167,7 @@ where /// Ok(()) /// } /// ``` - pub fn new(nrow: usize, ncol: usize, max_nnz: usize, symmetry: Option) -> Result { + pub fn new(nrow: usize, ncol: usize, max_nnz: usize, symmetry: Option) -> Result { if nrow < 1 { return Err("nrow must be ≥ 1"); } @@ -178,7 +178,7 @@ where return Err("max_nnz must be ≥ 1"); } Ok(NumCooMatrix { - symmetry: if let Some(v) = symmetry { v } else { Symmetry::No }, + symmetry: if let Some(v) = symmetry { v } else { Sym::No }, nrow, ncol, nnz: 0, @@ -242,7 +242,7 @@ where row_indices: Vec, col_indices: Vec, values: Vec, - symmetry: Option, + symmetry: Option, ) -> Result { if nrow < 1 { return Err("nrow must be ≥ 1"); @@ -271,7 +271,7 @@ where } } Ok(NumCooMatrix { - symmetry: if let Some(v) = symmetry { v } else { Symmetry::No }, + symmetry: if let Some(v) = symmetry { v } else { Sym::No }, nrow, ncol, nnz, @@ -325,7 +325,7 @@ where if self.nnz >= self.max_nnz { return Err("COO matrix: max number of items has been reached"); } - if self.symmetry != Symmetry::No { + if self.symmetry != Sym::No { if self.symmetry.lower() { if j > i { return Err("COO matrix: j > i is incorrect for lower triangular storage"); @@ -634,18 +634,18 @@ where /// Ok(()) /// } /// ``` - pub fn get_info(&self) -> (usize, usize, usize, Symmetry) { + pub fn get_info(&self) -> (usize, usize, usize, Sym) { (self.nrow, self.ncol, self.nnz, self.symmetry) } /// Returns the storage corresponding to the symmetry type (if any) pub fn get_storage(&self) -> Storage { - Symmetry::storage(self.symmetry) + Sym::storage(self.symmetry) } /// Returns whether the symmetry flag corresponds to a symmetric matrix or not pub fn get_symmetric(&self) -> bool { - self.symmetry != Symmetry::No + self.symmetry != Sym::No } /// Get an access to the row indices @@ -692,7 +692,7 @@ where #[cfg(test)] mod tests { use super::NumCooMatrix; - use crate::{Samples, Storage, Symmetry}; + use crate::{Samples, Storage, Sym}; use num_complex::Complex64; use russell_lab::{complex_vec_approx_eq, cpx, vec_approx_eq, ComplexVector, NumMatrix, NumVector}; @@ -709,7 +709,7 @@ mod tests { #[test] fn new_works() { let coo = NumCooMatrix::::new(1, 1, 3, None).unwrap(); - assert_eq!(coo.symmetry, Symmetry::No); + assert_eq!(coo.symmetry, Sym::No); assert_eq!(coo.nrow, 1); assert_eq!(coo.ncol, 1); assert_eq!(coo.nnz, 0); @@ -736,7 +736,7 @@ mod tests { #[test] fn from_works() { let coo = NumCooMatrix::::from(1, 1, vec![0], vec![0], vec![123.0], None).unwrap(); - assert_eq!(coo.symmetry, Symmetry::No); + assert_eq!(coo.symmetry, Sym::No); assert_eq!(coo.nrow, 1); assert_eq!(coo.ncol, 1); assert_eq!(coo.nnz, 1); @@ -744,9 +744,9 @@ mod tests { assert_eq!(coo.indices_i, &[0]); assert_eq!(coo.indices_j, &[0]); assert_eq!(coo.values, &[123.0]); - let sym = Some(Symmetry::new_general_full()); + let sym = Some(Sym::new_general_full()); let coo = NumCooMatrix::::from(1, 1, vec![0], vec![0], vec![123.0], sym).unwrap(); - assert_eq!(coo.symmetry, Symmetry::General(Storage::Full)); + assert_eq!(coo.symmetry, Sym::General(Storage::Full)); } #[test] @@ -756,7 +756,7 @@ mod tests { assert_eq!(nrow, 1); assert_eq!(ncol, 2); assert_eq!(nnz, 0); - assert_eq!(symmetry, Symmetry::No); + assert_eq!(symmetry, Sym::No); } #[test] @@ -775,13 +775,13 @@ mod tests { coo.put(0, 0, 0).err(), Some("COO matrix: max number of items has been reached") ); - let sym = Some(Symmetry::General(Storage::Lower)); + let sym = Some(Sym::General(Storage::Lower)); let mut coo = NumCooMatrix::::new(2, 2, 4, sym).unwrap(); assert_eq!( coo.put(0, 1, 0).err(), Some("COO matrix: j > i is incorrect for lower triangular storage") ); - let sym = Some(Symmetry::General(Storage::Upper)); + let sym = Some(Sym::General(Storage::Upper)); let mut coo = NumCooMatrix::::new(2, 2, 4, sym).unwrap(); assert_eq!( coo.put(1, 0, 0).err(), @@ -898,7 +898,7 @@ mod tests { #[test] fn to_dense_symmetric_lower_works() { - let sym = Some(Symmetry::General(Storage::Lower)); + let sym = Some(Sym::General(Storage::Lower)); let mut coo = NumCooMatrix::::new(3, 3, 4, sym).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(1, 0, 2.0).unwrap(); @@ -916,7 +916,7 @@ mod tests { #[test] fn to_dense_symmetric_upper_works() { - let sym = Some(Symmetry::General(Storage::Upper)); + let sym = Some(Sym::General(Storage::Upper)); let mut coo = NumCooMatrix::::new(3, 3, 4, sym).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(0, 1, 2.0).unwrap(); @@ -986,7 +986,7 @@ mod tests { // 3 1 1 7 // 2 1 5 1 8 let (nrow, ncol, nnz) = (5, 5, 15); - let sym = Some(Symmetry::General(Storage::Lower)); + let sym = Some(Sym::General(Storage::Lower)); let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, sym).unwrap(); coo.put(0, 0, 2.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); @@ -1022,7 +1022,7 @@ mod tests { // 3 1 1 7 1 // 2 1 5 1 8 let (nrow, ncol, nnz) = (5, 5, 25); - let sym = Some(Symmetry::General(Storage::Full)); + let sym = Some(Sym::General(Storage::Full)); let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, sym).unwrap(); coo.put(0, 0, 2.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); @@ -1066,7 +1066,7 @@ mod tests { // -1 2 -1 => -1 2 // -1 2 -1 2 let (nrow, ncol, nnz) = (3, 3, 5); - let sym = Some(Symmetry::PositiveDefinite(Storage::Lower)); + let sym = Some(Sym::PositiveDefinite(Storage::Lower)); let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, sym).unwrap(); coo.put(0, 0, 2.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); @@ -1104,7 +1104,7 @@ mod tests { #[test] fn assign_capture_errors() { - let sym = Some(Symmetry::General(Storage::Full)); + let sym = Some(Sym::General(Storage::Full)); let nnz_a = 1; let nnz_b = 2; // wrong: must be ≤ nnz_a let mut a_1x2 = NumCooMatrix::::new(1, 2, nnz_a, None).unwrap(); @@ -1156,7 +1156,7 @@ mod tests { #[test] fn augment_capture_errors() { - let sym = Some(Symmetry::General(Storage::Full)); + let sym = Some(Sym::General(Storage::Full)); let nnz_a = 1; let nnz_b = 1; let mut a_1x2 = NumCooMatrix::::new(1, 2, nnz_a /* + nnz_b */, None).unwrap(); @@ -1209,14 +1209,14 @@ mod tests { #[test] fn getters_are_correct() { let (coo, _, _, _) = Samples::rectangular_1x2(false, false); - assert_eq!(coo.get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(coo.get_info(), (1, 2, 2, Sym::No)); assert_eq!(coo.get_storage(), Storage::Full); assert_eq!(coo.get_symmetric(), false); assert_eq!(coo.get_row_indices(), &[0, 0]); assert_eq!(coo.get_col_indices(), &[0, 1]); assert_eq!(coo.get_values(), &[10.0, 20.0]); - let sym = Some(Symmetry::new_general_full()); + let sym = Some(Sym::new_general_full()); let coo = NumCooMatrix::::new(2, 2, 2, sym).unwrap(); assert_eq!(coo.get_symmetric(), true); diff --git a/russell_sparse/src/csc_matrix.rs b/russell_sparse/src/csc_matrix.rs index f2b1e849..cc9d7204 100644 --- a/russell_sparse/src/csc_matrix.rs +++ b/russell_sparse/src/csc_matrix.rs @@ -1,4 +1,4 @@ -use super::{to_i32, NumCooMatrix, NumCsrMatrix, Symmetry}; +use super::{to_i32, NumCooMatrix, NumCsrMatrix, Sym}; use crate::StrError; use num_traits::{Num, NumCast}; use russell_lab::{NumMatrix, NumVector}; @@ -55,7 +55,7 @@ where T: AddAssign + MulAssign + Num + NumCast + Copy + DeserializeOwned + Serialize, { /// Defines the symmetry and storage: lower-triangular, upper-triangular, full-matrix - pub(crate) symmetry: Symmetry, + pub(crate) symmetry: Sym, /// Holds the number of rows (must fit i32) pub(crate) nrow: usize, @@ -195,7 +195,7 @@ where col_pointers: Vec, row_indices: Vec, values: Vec, - symmetry: Option, + symmetry: Option, ) -> Result { if nrow < 1 { return Err("nrow must be ≥ 1"); @@ -242,7 +242,7 @@ where } } Ok(NumCscMatrix { - symmetry: if let Some(v) = symmetry { v } else { Symmetry::No }, + symmetry: if let Some(v) = symmetry { v } else { Sym::No }, nrow, ncol, col_pointers, @@ -780,7 +780,7 @@ where /// Ok(()) /// } /// ``` - pub fn get_info(&self) -> (usize, usize, usize, Symmetry) { + pub fn get_info(&self) -> (usize, usize, usize, Sym) { ( self.nrow, self.ncol, @@ -839,7 +839,7 @@ where #[cfg(test)] mod tests { use super::NumCscMatrix; - use crate::{CooMatrix, Samples, Storage, Symmetry}; + use crate::{CooMatrix, Samples, Storage, Sym}; use num_complex::Complex64; use russell_lab::{complex_vec_approx_eq, cpx, vec_approx_eq, ComplexVector, Matrix, Vector}; @@ -904,7 +904,7 @@ mod tests { fn new_works() { let (_, csc_correct, _, _) = Samples::rectangular_1x2(false, false); let csc = NumCscMatrix::::new(1, 2, vec![0, 1, 2], vec![0, 0], vec![10.0, 20.0], None).unwrap(); - assert_eq!(csc.symmetry, Symmetry::No); + assert_eq!(csc.symmetry, Sym::No); assert_eq!(csc.nrow, 1); assert_eq!(csc.ncol, 2); assert_eq!(&csc.col_pointers, &csc_correct.col_pointers); @@ -1011,8 +1011,8 @@ mod tests { fn update_from_coo_captures_errors() { let (coo, _, _, _) = Samples::rectangular_1x2(false, false, ); let mut csc = NumCscMatrix::::from_coo(&coo).unwrap(); - let yes = Symmetry::General(Storage::Lower); - let no = Symmetry::No; + let yes = Sym::General(Storage::Lower); + let no = Sym::No; assert_eq!(csc.update_from_coo(&CooMatrix { symmetry: yes, nrow: 1, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.symmetry must be equal to csc.symmetry")); assert_eq!(csc.update_from_coo(&CooMatrix { symmetry: no, nrow: 2, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.nrow must be equal to csc.nrow")); assert_eq!(csc.update_from_coo(&CooMatrix { symmetry: no, nrow: 1, ncol: 1, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.ncol must be equal to csc.ncol")); @@ -1240,20 +1240,20 @@ mod tests { #[test] fn getters_are_correct() { let (_, csc, _, _) = Samples::rectangular_1x2(false, false); - assert_eq!(csc.get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(csc.get_info(), (1, 2, 2, Sym::No)); assert_eq!(csc.get_col_pointers(), &[0, 1, 2]); assert_eq!(csc.get_row_indices(), &[0, 0]); assert_eq!(csc.get_values(), &[10.0, 20.0]); // with duplicates let (coo, _, _, _) = Samples::rectangular_1x2(false, false); let csc = NumCscMatrix::::from_coo(&coo).unwrap(); - assert_eq!(csc.get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(csc.get_info(), (1, 2, 2, Sym::No)); assert_eq!(csc.get_col_pointers(), &[0, 1, 2]); assert_eq!(csc.get_row_indices(), &[0, 0]); assert_eq!(csc.get_values(), &[10.0, 20.0]); // mutable let mut csc = NumCscMatrix:: { - symmetry: Symmetry::No, + symmetry: Sym::No, nrow: 1, ncol: 2, values: vec![10.0, 20.0], diff --git a/russell_sparse/src/csr_matrix.rs b/russell_sparse/src/csr_matrix.rs index af0e5f42..f247710f 100644 --- a/russell_sparse/src/csr_matrix.rs +++ b/russell_sparse/src/csr_matrix.rs @@ -1,4 +1,4 @@ -use super::{to_i32, NumCooMatrix, NumCscMatrix, Symmetry}; +use super::{to_i32, NumCooMatrix, NumCscMatrix, Sym}; use crate::StrError; use num_traits::{Num, NumCast}; use russell_lab::{NumMatrix, NumVector}; @@ -55,7 +55,7 @@ where T: AddAssign + MulAssign + Num + NumCast + Copy + DeserializeOwned + Serialize, { /// Defines the symmetry and storage: lower-triangular, upper-triangular, full-matrix - pub(crate) symmetry: Symmetry, + pub(crate) symmetry: Sym, /// Holds the number of rows (must fit i32) pub(crate) nrow: usize, @@ -191,7 +191,7 @@ where row_pointers: Vec, col_indices: Vec, values: Vec, - symmetry: Option, + symmetry: Option, ) -> Result { if nrow < 1 { return Err("nrow must be ≥ 1"); @@ -238,7 +238,7 @@ where } } Ok(NumCsrMatrix { - symmetry: if let Some(v) = symmetry { v } else { Symmetry::No }, + symmetry: if let Some(v) = symmetry { v } else { Sym::No }, nrow, ncol, row_pointers, @@ -754,7 +754,7 @@ where /// Ok(()) /// } /// ``` - pub fn get_info(&self) -> (usize, usize, usize, Symmetry) { + pub fn get_info(&self) -> (usize, usize, usize, Sym) { ( self.nrow, self.ncol, @@ -813,7 +813,7 @@ where #[cfg(test)] mod tests { use super::NumCsrMatrix; - use crate::{CooMatrix, Samples, Storage, Symmetry}; + use crate::{CooMatrix, Samples, Storage, Sym}; use num_complex::Complex64; use russell_lab::{complex_vec_approx_eq, cpx, vec_approx_eq, ComplexVector, Matrix, Vector}; @@ -877,7 +877,7 @@ mod tests { fn new_works() { let (_, _, csr_correct, _) = Samples::rectangular_1x2(false, false); let csr = NumCsrMatrix::::new(1, 2, vec![0, 2], vec![0, 1], vec![10.0, 20.0], None).unwrap(); - assert_eq!(csr.symmetry, Symmetry::No); + assert_eq!(csr.symmetry, Sym::No); assert_eq!(csr.nrow, 1); assert_eq!(csr.ncol, 2); assert_eq!(&csr.row_pointers, &csr_correct.row_pointers); @@ -999,8 +999,8 @@ mod tests { fn update_from_coo_captures_errors() { let (coo, _, _, _) = Samples::rectangular_1x2(false, false); let mut csr = NumCsrMatrix::::from_coo(&coo).unwrap(); - let yes = Symmetry::General(Storage::Lower); - let no = Symmetry::No; + let yes = Sym::General(Storage::Lower); + let no = Sym::No; assert_eq!(csr.update_from_coo(&CooMatrix { symmetry: yes, nrow: 1, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.symmetry must be equal to csr.symmetry")); assert_eq!(csr.update_from_coo(&CooMatrix { symmetry: no, nrow: 2, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.nrow must be equal to csr.nrow")); assert_eq!(csr.update_from_coo(&CooMatrix { symmetry: no, nrow: 1, ncol: 1, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.ncol must be equal to csr.ncol")); @@ -1228,20 +1228,20 @@ mod tests { #[test] fn getters_are_correct() { let (_, _, csr, _) = Samples::rectangular_1x2(false, false); - assert_eq!(csr.get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(csr.get_info(), (1, 2, 2, Sym::No)); assert_eq!(csr.get_row_pointers(), &[0, 2]); assert_eq!(csr.get_col_indices(), &[0, 1]); assert_eq!(csr.get_values(), &[10.0, 20.0]); // with duplicates let (coo, _, _, _) = Samples::rectangular_1x2(false, false); let csr = NumCsrMatrix::::from_coo(&coo).unwrap(); - assert_eq!(csr.get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(csr.get_info(), (1, 2, 2, Sym::No)); assert_eq!(csr.get_row_pointers(), &[0, 2]); assert_eq!(csr.get_col_indices(), &[0, 1]); assert_eq!(csr.get_values(), &[10.0, 20.0]); // mutable let mut csr = NumCsrMatrix:: { - symmetry: Symmetry::No, + symmetry: Sym::No, nrow: 1, ncol: 2, values: vec![10.0, 20.0], diff --git a/russell_sparse/src/enums.rs b/russell_sparse/src/enums.rs index 9c82de93..c8d702fd 100644 --- a/russell_sparse/src/enums.rs +++ b/russell_sparse/src/enums.rs @@ -30,7 +30,7 @@ pub enum Storage { /// Specifies the type of matrix symmetry #[derive(Clone, Copy, Debug, Eq, PartialEq, Deserialize, Serialize)] -pub enum Symmetry { +pub enum Sym { /// Unknown symmetry (possibly unsymmetric) No, @@ -176,53 +176,53 @@ impl Genie { } /// Returns the solver's required symmetry/storage configuration - pub fn symmetry(&self, symmetric: bool, positive_definite: bool) -> Symmetry { + pub fn symmetry(&self, symmetric: bool, positive_definite: bool) -> Sym { let storage = self.storage(); if symmetric || positive_definite { if positive_definite { - Symmetry::PositiveDefinite(storage) + Sym::PositiveDefinite(storage) } else { - Symmetry::General(storage) + Sym::General(storage) } } else { - Symmetry::No + Sym::No } } } -impl Symmetry { +impl Sym { /// Returns a new general symmetry flag with lower storage pub fn new_general_lower() -> Self { - Symmetry::General(Storage::Lower) + Sym::General(Storage::Lower) } /// Returns a new general symmetry flag with upper storage pub fn new_general_upper() -> Self { - Symmetry::General(Storage::Upper) + Sym::General(Storage::Upper) } /// Returns a new general symmetry flag with full storage pub fn new_general_full() -> Self { - Symmetry::General(Storage::Full) + Sym::General(Storage::Full) } /// Returns a new positive-definite symmetry flag with lower storage pub fn new_pos_def_lower() -> Self { - Symmetry::PositiveDefinite(Storage::Lower) + Sym::PositiveDefinite(Storage::Lower) } /// Returns a new positive-definite symmetry flag with upper storage pub fn new_pos_def_upper() -> Self { - Symmetry::PositiveDefinite(Storage::Upper) + Sym::PositiveDefinite(Storage::Upper) } /// Returns a new positive-definite symmetry flag with full storage pub fn new_pos_def_full() -> Self { - Symmetry::PositiveDefinite(Storage::Full) + Sym::PositiveDefinite(Storage::Full) } /// Returns which type of storage is used, if symmetric - pub fn storage(symmetry: Symmetry) -> Storage { + pub fn storage(symmetry: Sym) -> Storage { if symmetry.lower() { Storage::Lower } else if symmetry.upper() { @@ -235,27 +235,27 @@ impl Symmetry { /// Returns true if the storage is triangular (lower or upper) pub fn triangular(&self) -> bool { match self { - Symmetry::No => false, - Symmetry::General(storage) => *storage != Storage::Full, - Symmetry::PositiveDefinite(storage) => *storage != Storage::Full, + Sym::No => false, + Sym::General(storage) => *storage != Storage::Full, + Sym::PositiveDefinite(storage) => *storage != Storage::Full, } } /// Returns true if the storage is lower triangular pub fn lower(&self) -> bool { match self { - Symmetry::No => false, - Symmetry::General(storage) => *storage == Storage::Lower, - Symmetry::PositiveDefinite(storage) => *storage == Storage::Lower, + Sym::No => false, + Sym::General(storage) => *storage == Storage::Lower, + Sym::PositiveDefinite(storage) => *storage == Storage::Lower, } } /// Returns true if the storage is upper triangular pub fn upper(&self) -> bool { match self { - Symmetry::No => false, - Symmetry::General(storage) => *storage == Storage::Upper, - Symmetry::PositiveDefinite(storage) => *storage == Storage::Upper, + Sym::No => false, + Sym::General(storage) => *storage == Storage::Upper, + Sym::PositiveDefinite(storage) => *storage == Storage::Upper, } } @@ -274,8 +274,8 @@ impl Symmetry { /// * `positive_definite` -- 1 if true, 0 otherwise pub fn status(&self, must_be_lower: bool, must_be_upper: bool) -> Result<(i32, i32), StrError> { match self { - Symmetry::No => Ok((0, 0)), - Symmetry::General(storage) => { + Sym::No => Ok((0, 0)), + Sym::General(storage) => { if must_be_lower && *storage != Storage::Lower { return Err("if the matrix is general symmetric, the required storage is lower triangular"); } @@ -284,7 +284,7 @@ impl Symmetry { } Ok((1, 0)) } - Symmetry::PositiveDefinite(storage) => { + Sym::PositiveDefinite(storage) => { if must_be_lower && *storage != Storage::Lower { return Err("if the matrix is positive-definite, the required storage is lower triangular"); } @@ -387,14 +387,14 @@ mod tests { let from_json: Storage = serde_json::from_str(&json).unwrap(); assert_eq!(from_json, storage); - let symmetry = Symmetry::PositiveDefinite(Storage::Lower); + let symmetry = Sym::PositiveDefinite(Storage::Lower); let copy = symmetry; let clone = symmetry.clone(); assert_eq!(format!("{:?}", symmetry), "PositiveDefinite(Lower)"); - assert_eq!(copy, Symmetry::PositiveDefinite(Storage::Lower)); - assert_eq!(clone, Symmetry::PositiveDefinite(Storage::Lower)); + assert_eq!(copy, Sym::PositiveDefinite(Storage::Lower)); + assert_eq!(clone, Sym::PositiveDefinite(Storage::Lower)); let json = serde_json::to_string(&symmetry).unwrap(); - let from_json: Symmetry = serde_json::from_str(&json).unwrap(); + let from_json: Sym = serde_json::from_str(&json).unwrap(); assert_eq!(from_json, symmetry); let handling = MMsymOption::LeaveAsLower; @@ -492,16 +492,16 @@ mod tests { let l = Storage::Lower; let f = Storage::Full; - let gl = Symmetry::General(l); - let gf = Symmetry::General(f); + let gl = Sym::General(l); + let gf = Sym::General(f); - let pl = Symmetry::PositiveDefinite(l); - let pf = Symmetry::PositiveDefinite(f); + let pl = Sym::PositiveDefinite(l); + let pf = Sym::PositiveDefinite(f); let genie = Genie::Mumps; assert_eq!(genie.to_string(), "mumps"); assert_eq!(genie.storage(), l); - assert_eq!(genie.symmetry(false, false), Symmetry::No); + assert_eq!(genie.symmetry(false, false), Sym::No); assert_eq!(genie.symmetry(true, false), gl); assert_eq!(genie.symmetry(false, true), pl); assert_eq!(genie.symmetry(true, true), pl); @@ -509,7 +509,7 @@ mod tests { let genie = Genie::Umfpack; assert_eq!(genie.to_string(), "umfpack"); assert_eq!(genie.storage(), f); - assert_eq!(genie.symmetry(false, false), Symmetry::No); + assert_eq!(genie.symmetry(false, false), Sym::No); assert_eq!(genie.symmetry(true, false), gf); assert_eq!(genie.symmetry(false, true), pf); assert_eq!(genie.symmetry(true, true), pf); @@ -521,20 +521,20 @@ mod tests { let u = Storage::Upper; let f = Storage::Full; - let gl = Symmetry::General(l); - let gu = Symmetry::General(u); - let gf = Symmetry::General(f); + let gl = Sym::General(l); + let gu = Sym::General(u); + let gf = Sym::General(f); - let pl = Symmetry::PositiveDefinite(l); - let pu = Symmetry::PositiveDefinite(u); - let pf = Symmetry::PositiveDefinite(f); + let pl = Sym::PositiveDefinite(l); + let pu = Sym::PositiveDefinite(u); + let pf = Sym::PositiveDefinite(f); - assert_eq!(Symmetry::new_general_lower(), gl); - assert_eq!(Symmetry::new_general_upper(), gu); - assert_eq!(Symmetry::new_general_full(), gf); - assert_eq!(Symmetry::storage(gl), Storage::Lower); - assert_eq!(Symmetry::storage(gu), Storage::Upper); - assert_eq!(Symmetry::storage(gf), Storage::Full); + assert_eq!(Sym::new_general_lower(), gl); + assert_eq!(Sym::new_general_upper(), gu); + assert_eq!(Sym::new_general_full(), gf); + assert_eq!(Sym::storage(gl), Storage::Lower); + assert_eq!(Sym::storage(gu), Storage::Upper); + assert_eq!(Sym::storage(gf), Storage::Full); assert_eq!(gl.triangular(), true); assert_eq!(gu.triangular(), true); assert_eq!(gf.triangular(), false); @@ -559,12 +559,12 @@ mod tests { Err("if the matrix is general symmetric, the required storage is lower triangular") ); - assert_eq!(Symmetry::new_pos_def_lower(), pl); - assert_eq!(Symmetry::new_pos_def_upper(), pu); - assert_eq!(Symmetry::new_pos_def_full(), pf); - assert_eq!(Symmetry::storage(pl), Storage::Lower); - assert_eq!(Symmetry::storage(pu), Storage::Upper); - assert_eq!(Symmetry::storage(pf), Storage::Full); + assert_eq!(Sym::new_pos_def_lower(), pl); + assert_eq!(Sym::new_pos_def_upper(), pu); + assert_eq!(Sym::new_pos_def_full(), pf); + assert_eq!(Sym::storage(pl), Storage::Lower); + assert_eq!(Sym::storage(pu), Storage::Upper); + assert_eq!(Sym::storage(pf), Storage::Full); assert_eq!(pl.triangular(), true); assert_eq!(pu.triangular(), true); assert_eq!(pf.triangular(), false); @@ -589,6 +589,6 @@ mod tests { Err("if the matrix is positive-definite, the required storage is lower triangular") ); - assert_eq!(Symmetry::storage(Symmetry::No), Storage::Full); + assert_eq!(Sym::storage(Sym::No), Storage::Full); } } diff --git a/russell_sparse/src/read_matrix_market.rs b/russell_sparse/src/read_matrix_market.rs index 9eb99d19..b7eebfef 100644 --- a/russell_sparse/src/read_matrix_market.rs +++ b/russell_sparse/src/read_matrix_market.rs @@ -1,4 +1,4 @@ -use super::{CooMatrix, MMsymOption, Storage, Symmetry}; +use super::{CooMatrix, MMsymOption, Storage, Sym}; use crate::StrError; use std::ffi::OsStr; use std::fs::File; @@ -354,9 +354,9 @@ where return Err("MatrixMarket data is invalid: the number of rows must be equal the number of columns for symmetric matrices"); } match symmetric_handling { - MMsymOption::LeaveAsLower => Some(Symmetry::General(Storage::Lower)), - MMsymOption::SwapToUpper => Some(Symmetry::General(Storage::Upper)), - MMsymOption::MakeItFull => Some(Symmetry::General(Storage::Full)), + MMsymOption::LeaveAsLower => Some(Sym::General(Storage::Lower)), + MMsymOption::SwapToUpper => Some(Sym::General(Storage::Upper)), + MMsymOption::MakeItFull => Some(Sym::General(Storage::Full)), } } else { None @@ -414,7 +414,7 @@ where #[cfg(test)] mod tests { use super::{read_matrix_market, MatrixMarketData}; - use crate::{MMsymOption, Storage, Symmetry}; + use crate::{MMsymOption, Storage, Sym}; use russell_lab::Matrix; #[test] @@ -564,7 +564,7 @@ mod tests { let h = MMsymOption::LeaveAsLower; let filepath = "./data/matrix_market/ok_general.mtx".to_string(); let coo = read_matrix_market(&filepath, h).unwrap(); - assert_eq!(coo.symmetry, Symmetry::No); + assert_eq!(coo.symmetry, Sym::No); assert_eq!((coo.nrow, coo.ncol, coo.nnz, coo.max_nnz), (5, 5, 12, 12)); assert_eq!(coo.indices_i, &[0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4]); assert_eq!(coo.indices_j, &[0, 0, 1, 1, 1, 2, 2, 2, 2, 3, 4, 4]); @@ -579,7 +579,7 @@ mod tests { let h = MMsymOption::LeaveAsLower; let filepath = "./data/matrix_market/ok_symmetric.mtx".to_string(); let coo = read_matrix_market(&filepath, h).unwrap(); - assert_eq!(coo.symmetry, Symmetry::General(Storage::Lower)); + assert_eq!(coo.symmetry, Sym::General(Storage::Lower)); assert_eq!((coo.nrow, coo.ncol, coo.nnz, coo.max_nnz), (5, 5, 15, 15)); assert_eq!(coo.indices_i, &[0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4]); assert_eq!(coo.indices_j, &[0, 1, 2, 3, 4, 0, 0, 0, 0, 1, 1, 1, 2, 2, 3]); @@ -594,7 +594,7 @@ mod tests { let h = MMsymOption::SwapToUpper; let filepath = "./data/matrix_market/ok_symmetric.mtx".to_string(); let coo = read_matrix_market(&filepath, h).unwrap(); - assert_eq!(coo.symmetry, Symmetry::General(Storage::Upper)); + assert_eq!(coo.symmetry, Sym::General(Storage::Upper)); assert_eq!((coo.nrow, coo.ncol, coo.nnz, coo.max_nnz), (5, 5, 15, 15)); assert_eq!(coo.indices_i, &[0, 1, 2, 3, 4, 0, 0, 0, 0, 1, 1, 1, 2, 2, 3]); assert_eq!(coo.indices_j, &[0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4]); @@ -609,7 +609,7 @@ mod tests { let h = MMsymOption::MakeItFull; let filepath = "./data/matrix_market/ok_symmetric_small.mtx".to_string(); let coo = read_matrix_market(&filepath, h).unwrap(); - assert_eq!(coo.symmetry, Symmetry::General(Storage::Full)); + assert_eq!(coo.symmetry, Sym::General(Storage::Full)); assert_eq!((coo.nrow, coo.ncol, coo.nnz, coo.max_nnz), (5, 5, 11, 14)); assert_eq!(coo.indices_i, &[0, 1, 0, 2, 1, 3, 2, 3, 4, 1, 4, 0, 0, 0]); assert_eq!(coo.indices_j, &[0, 0, 1, 1, 2, 2, 3, 3, 1, 4, 4, 0, 0, 0]); diff --git a/russell_sparse/src/samples.rs b/russell_sparse/src/samples.rs index ecf0e952..ad23f452 100644 --- a/russell_sparse/src/samples.rs +++ b/russell_sparse/src/samples.rs @@ -1,5 +1,5 @@ use crate::{ComplexCooMatrix, ComplexCscMatrix, ComplexCsrMatrix}; -use crate::{CooMatrix, CscMatrix, CsrMatrix, Storage, Symmetry}; +use crate::{CooMatrix, CscMatrix, CsrMatrix, Storage, Sym}; use num_complex::Complex64; use russell_lab::cpx; @@ -72,7 +72,7 @@ impl Samples { /// ``` pub fn positive_definite_3x3() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { let (nrow, ncol, nnz) = (3, 3, 6); - let sym = Some(Symmetry::PositiveDefinite(Storage::Lower)); + let sym = Some(Sym::PositiveDefinite(Storage::Lower)); let mut coo = CooMatrix::new(nrow, ncol, nnz, sym).unwrap(); coo.put(1, 0, -0.5).unwrap(); // duplicate coo.put(0, 0, 2.0).unwrap(); @@ -118,7 +118,7 @@ impl Samples { /// ``` pub fn complex_symmetric_3x3_lower() -> (ComplexCooMatrix, ComplexCscMatrix, ComplexCsrMatrix, Complex64) { let (nrow, ncol, nnz) = (3, 3, 6); - let sym = Some(Symmetry::General(Storage::Lower)); + let sym = Some(Sym::General(Storage::Lower)); let mut coo = ComplexCooMatrix::new(nrow, ncol, nnz, sym).unwrap(); coo.put(1, 0, cpx!(-0.5, -0.5)).unwrap(); // duplicate coo.put(0, 0, cpx!(2.0, 1.0)).unwrap(); @@ -166,7 +166,7 @@ impl Samples { /// ``` pub fn complex_symmetric_3x3_full() -> (ComplexCooMatrix, ComplexCscMatrix, ComplexCsrMatrix, Complex64) { let (nrow, ncol, nnz) = (3, 3, 8); - let sym = Some(Symmetry::General(Storage::Full)); + let sym = Some(Sym::General(Storage::Full)); let mut coo = ComplexCooMatrix::new(nrow, ncol, nnz, sym).unwrap(); coo.put(1, 0, cpx!(-0.5, -0.5)).unwrap(); // duplicate coo.put(0, 0, cpx!(2.0, 1.0)).unwrap(); @@ -218,7 +218,7 @@ impl Samples { /// ``` pub fn lower_symmetric_5x5() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { let (nrow, ncol, nnz) = (5, 5, 18); - let sym = Some(Symmetry::PositiveDefinite(Storage::Lower)); + let sym = Some(Sym::PositiveDefinite(Storage::Lower)); let mut coo = CooMatrix::new(nrow, ncol, nnz, sym).unwrap(); coo.put(1, 1, 2.0).unwrap(); coo.put(4, 2, 2.5).unwrap(); // duplicate @@ -655,7 +655,7 @@ impl Samples { /// x_correct = vec![-979.0 / 3.0, 983.0, 1961.0 / 12.0, 398.0, 123.0 / 2.0]; /// ``` pub fn mkl_positive_definite_5x5_lower() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = Some(Symmetry::PositiveDefinite(Storage::Lower)); + let sym = Some(Sym::PositiveDefinite(Storage::Lower)); let nrow = 5; let ncol = 5; let mut coo = CooMatrix::new(nrow, ncol, 9, sym).unwrap(); @@ -729,7 +729,7 @@ impl Samples { /// x_correct = vec![-979.0 / 3.0, 983.0, 1961.0 / 12.0, 398.0, 123.0 / 2.0]; /// ``` pub fn mkl_positive_definite_5x5_upper() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = Some(Symmetry::PositiveDefinite(Storage::Upper)); + let sym = Some(Sym::PositiveDefinite(Storage::Upper)); let nrow = 5; let ncol = 5; let mut coo = CooMatrix::new(nrow, ncol, 9, sym).unwrap(); @@ -806,7 +806,7 @@ impl Samples { shuffle_coo_entries: bool, duplicate_coo_entries: bool, ) -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = Some(Symmetry::General(Storage::Lower)); + let sym = Some(Sym::General(Storage::Lower)); let nrow = 5; let ncol = 5; let max_nnz = 13; @@ -928,7 +928,7 @@ impl Samples { shuffle_coo_entries: bool, duplicate_coo_entries: bool, ) -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = Some(Symmetry::General(Storage::Upper)); + let sym = Some(Sym::General(Storage::Upper)); let nrow = 5; let ncol = 5; let max_nnz = 15; @@ -1046,7 +1046,7 @@ impl Samples { /// x_correct = vec![-979.0 / 3.0, 983.0, 1961.0 / 12.0, 398.0, 123.0 / 2.0]; /// ``` pub fn mkl_symmetric_5x5_full() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = Some(Symmetry::General(Storage::Full)); + let sym = Some(Sym::General(Storage::Full)); let nrow = 5; let ncol = 5; let max_nnz = 13; diff --git a/russell_sparse/src/solver_mumps.rs b/russell_sparse/src/solver_mumps.rs index b0fcdca0..63980402 100644 --- a/russell_sparse/src/solver_mumps.rs +++ b/russell_sparse/src/solver_mumps.rs @@ -1,4 +1,4 @@ -use super::{LinSolParams, LinSolTrait, Ordering, Scaling, SparseMatrix, StatsLinSol, Symmetry}; +use super::{LinSolParams, LinSolTrait, Ordering, Scaling, SparseMatrix, StatsLinSol, Sym}; use crate::auxiliary_and_constants::*; use crate::StrError; use russell_lab::{using_intel_mkl, vec_copy, Stopwatch, Vector}; @@ -124,7 +124,7 @@ pub struct SolverMUMPS { factorized: bool, /// Holds the symmetry type used in the initialize - initialized_symmetry: Symmetry, + initialized_symmetry: Sym, /// Holds the matrix dimension saved in initialize initialized_ndim: usize, @@ -199,7 +199,7 @@ impl SolverMUMPS { solver, initialized: false, factorized: false, - initialized_symmetry: Symmetry::No, + initialized_symmetry: Sym::No, initialized_ndim: 0, initialized_nnz: 0, effective_ordering: -1, @@ -616,7 +616,7 @@ pub(crate) fn handle_mumps_error_code(err: i32) -> StrError { #[cfg(test)] mod tests { use super::*; - use crate::{CooMatrix, LinSolParams, LinSolTrait, Ordering, Samples, Scaling, SparseMatrix, Storage, Symmetry}; + use crate::{CooMatrix, LinSolParams, LinSolTrait, Ordering, Samples, Scaling, SparseMatrix, Storage, Sym}; use russell_lab::{approx_eq, vec_approx_eq, Vector}; use serial_test::serial; @@ -666,7 +666,7 @@ mod tests { // ... factorize once => OK solver.factorize(&mut mat, None).unwrap(); // ... change matrix (symmetry) - let mut coo = CooMatrix::new(2, 2, 2, Some(Symmetry::General(Storage::Full))).unwrap(); + let mut coo = CooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); @@ -734,7 +734,7 @@ mod tests { ); // wrong symmetry let rhs = Vector::new(2); - let mut coo_wrong = CooMatrix::new(2, 2, 2, Some(Symmetry::General(Storage::Full))).unwrap(); + let mut coo_wrong = CooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); coo_wrong.put(0, 0, 123.0).unwrap(); coo_wrong.put(1, 1, 456.0).unwrap(); let mut mat_wrong = SparseMatrix::from_coo(coo_wrong); diff --git a/russell_sparse/src/solver_umfpack.rs b/russell_sparse/src/solver_umfpack.rs index eab1dd46..194c4cdc 100644 --- a/russell_sparse/src/solver_umfpack.rs +++ b/russell_sparse/src/solver_umfpack.rs @@ -1,4 +1,4 @@ -use super::{LinSolParams, LinSolTrait, Ordering, Scaling, SparseMatrix, StatsLinSol, Symmetry}; +use super::{LinSolParams, LinSolTrait, Ordering, Scaling, SparseMatrix, StatsLinSol, Sym}; use crate::auxiliary_and_constants::*; use crate::StrError; use russell_lab::{Stopwatch, Vector}; @@ -75,7 +75,7 @@ pub struct SolverUMFPACK { factorized: bool, /// Holds the symmetry type used in initialize - initialized_symmetry: Symmetry, + initialized_symmetry: Sym, /// Holds the matrix dimension saved in initialize initialized_ndim: usize, @@ -139,7 +139,7 @@ impl SolverUMFPACK { solver, initialized: false, factorized: false, - initialized_symmetry: Symmetry::No, + initialized_symmetry: Sym::No, initialized_ndim: 0, initialized_nnz: 0, effective_strategy: -1, @@ -465,7 +465,7 @@ pub(crate) fn handle_umfpack_error_code(err: i32) -> StrError { #[cfg(test)] mod tests { use super::*; - use crate::{CooMatrix, LinSolParams, LinSolTrait, Ordering, Samples, Scaling, SparseMatrix, Storage, Symmetry}; + use crate::{CooMatrix, LinSolParams, LinSolTrait, Ordering, Samples, Scaling, SparseMatrix, Storage, Sym}; use russell_lab::{approx_eq, vec_approx_eq, Vector}; #[test] @@ -510,7 +510,7 @@ mod tests { // ... factorize once => OK solver.factorize(&mut mat, None).unwrap(); // ... change matrix (symmetry) - let mut coo = CooMatrix::new(2, 2, 2, Some(Symmetry::General(Storage::Full))).unwrap(); + let mut coo = CooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); @@ -601,7 +601,7 @@ mod tests { ); // wrong symmetry let rhs = Vector::new(2); - let mut coo_wrong = CooMatrix::new(2, 2, 2, Some(Symmetry::General(Storage::Full))).unwrap(); + let mut coo_wrong = CooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); coo_wrong.put(0, 0, 123.0).unwrap(); coo_wrong.put(1, 1, 456.0).unwrap(); let mut mat_wrong = SparseMatrix::from_coo(coo_wrong); diff --git a/russell_sparse/src/sparse_matrix.rs b/russell_sparse/src/sparse_matrix.rs index c84ff200..4b9821fa 100644 --- a/russell_sparse/src/sparse_matrix.rs +++ b/russell_sparse/src/sparse_matrix.rs @@ -1,4 +1,4 @@ -use super::{NumCooMatrix, NumCscMatrix, NumCsrMatrix, Symmetry}; +use super::{NumCooMatrix, NumCscMatrix, NumCsrMatrix, Sym}; use crate::StrError; use num_traits::{Num, NumCast}; use russell_lab::{NumMatrix, NumVector}; @@ -59,7 +59,7 @@ where /// * `max_nnz` -- (≥ 1) Maximum number of entries ≥ nnz (number of non-zeros), /// including entries with repeated indices. (must be fit i32) /// * `symmetry` -- Defines the symmetry/storage, if any - pub fn new_coo(nrow: usize, ncol: usize, max_nnz: usize, symmetry: Option) -> Result { + pub fn new_coo(nrow: usize, ncol: usize, max_nnz: usize, symmetry: Option) -> Result { Ok(NumSparseMatrix { coo: Some(NumCooMatrix::new(nrow, ncol, max_nnz, symmetry)?), csc: None, @@ -98,7 +98,7 @@ where col_pointers: Vec, row_indices: Vec, values: Vec, - symmetry: Option, + symmetry: Option, ) -> Result { Ok(NumSparseMatrix { coo: None, @@ -145,7 +145,7 @@ where row_pointers: Vec, col_indices: Vec, values: Vec, - symmetry: Option, + symmetry: Option, ) -> Result { Ok(NumSparseMatrix { coo: None, @@ -193,7 +193,7 @@ where /// Returns `(nrow, ncol, nnz, symmetry)` /// /// **Priority**: CSC -> CSR -> COO - pub fn get_info(&self) -> (usize, usize, usize, Symmetry) { + pub fn get_info(&self) -> (usize, usize, usize, Sym) { match &self.csc { Some(csc) => csc.get_info(), None => match &self.csr { @@ -437,7 +437,7 @@ where #[cfg(test)] mod tests { use super::NumSparseMatrix; - use crate::{Samples, Symmetry}; + use crate::{Samples, Sym}; use russell_lab::{vec_approx_eq, Matrix, Vector}; #[test] @@ -471,22 +471,22 @@ mod tests { let mut wrong = Vector::new(2); // COO let coo_mat = NumSparseMatrix::::from_coo(coo); - assert_eq!(coo_mat.get_info(), (1, 2, 2, Symmetry::No)); - assert_eq!(coo_mat.get_coo().unwrap().get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(coo_mat.get_info(), (1, 2, 2, Sym::No)); + assert_eq!(coo_mat.get_coo().unwrap().get_info(), (1, 2, 2, Sym::No)); assert_eq!(coo_mat.get_csc().err(), Some("CSC matrix is not available")); assert_eq!(coo_mat.get_csr().err(), Some("CSR matrix is not available")); assert_eq!(coo_mat.get_values(), &[10.0, 20.0]); // CSC let csc_mat = NumSparseMatrix::::from_csc(csc); - assert_eq!(csc_mat.get_info(), (1, 2, 2, Symmetry::No)); - assert_eq!(csc_mat.get_csc().unwrap().get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(csc_mat.get_info(), (1, 2, 2, Sym::No)); + assert_eq!(csc_mat.get_csc().unwrap().get_info(), (1, 2, 2, Sym::No)); assert_eq!(csc_mat.get_coo().err(), Some("COO matrix is not available")); assert_eq!(csc_mat.get_csr().err(), Some("CSR matrix is not available")); assert_eq!(csc_mat.get_values(), &[10.0, 20.0]); // CSR let csr_mat = NumSparseMatrix::::from_csr(csr); - assert_eq!(csr_mat.get_info(), (1, 2, 2, Symmetry::No)); - assert_eq!(csr_mat.get_csr().unwrap().get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(csr_mat.get_info(), (1, 2, 2, Sym::No)); + assert_eq!(csr_mat.get_csr().unwrap().get_info(), (1, 2, 2, Sym::No)); assert_eq!(csr_mat.get_csc().err(), Some("CSC matrix is not available")); assert_eq!(csr_mat.get_coo().err(), Some("COO matrix is not available")); assert_eq!(csr_mat.get_values(), &[10.0, 20.0]); @@ -522,7 +522,7 @@ mod tests { wrong.put(0, 0, 3.0).unwrap(); // COO let mut coo_mat = NumSparseMatrix::::from_coo(coo); - assert_eq!(coo_mat.get_coo_mut().unwrap().get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(coo_mat.get_coo_mut().unwrap().get_info(), (1, 2, 2, Sym::No)); assert_eq!(coo_mat.get_csc_mut().err(), Some("CSC matrix is not available")); assert_eq!(coo_mat.get_csr_mut().err(), Some("CSR matrix is not available")); let mut empty = NumSparseMatrix::::new_coo(1, 1, 1, None).unwrap(); @@ -530,13 +530,10 @@ mod tests { assert_eq!(empty.get_csr_or_from_coo().err(), Some("COO to CSR requires nnz > 0")); // CSC let mut csc_mat = NumSparseMatrix::::from_csc(csc); - assert_eq!(csc_mat.get_csc_mut().unwrap().get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(csc_mat.get_csc_mut().unwrap().get_info(), (1, 2, 2, Sym::No)); assert_eq!(csc_mat.get_coo_mut().err(), Some("COO matrix is not available")); assert_eq!(csc_mat.get_csr_mut().err(), Some("CSR matrix is not available")); - assert_eq!( - csc_mat.get_csc_or_from_coo().unwrap().get_info(), - (1, 2, 2, Symmetry::No) - ); + assert_eq!(csc_mat.get_csc_or_from_coo().unwrap().get_info(), (1, 2, 2, Sym::No)); assert_eq!( csc_mat.get_csr_or_from_coo().err(), Some("CSR is not available and COO matrix is not available to convert to CSR") @@ -559,13 +556,10 @@ mod tests { ); // CSR let mut csr_mat = NumSparseMatrix::::from_csr(csr); - assert_eq!(csr_mat.get_csr_mut().unwrap().get_info(), (1, 2, 2, Symmetry::No)); + assert_eq!(csr_mat.get_csr_mut().unwrap().get_info(), (1, 2, 2, Sym::No)); assert_eq!(csr_mat.get_csc_mut().err(), Some("CSC matrix is not available")); assert_eq!(csr_mat.get_coo_mut().err(), Some("COO matrix is not available")); - assert_eq!( - csr_mat.get_csr_or_from_coo().unwrap().get_info(), - (1, 2, 2, Symmetry::No) - ); + assert_eq!(csr_mat.get_csr_or_from_coo().unwrap().get_info(), (1, 2, 2, Sym::No)); assert_eq!( csr_mat.get_csc_or_from_coo().err(), Some("CSC is not available and COO matrix is not available to convert to CSC") diff --git a/russell_sparse/src/write_matrix_market.rs b/russell_sparse/src/write_matrix_market.rs index cacd2580..7448417c 100644 --- a/russell_sparse/src/write_matrix_market.rs +++ b/russell_sparse/src/write_matrix_market.rs @@ -1,4 +1,4 @@ -use super::{CscMatrix, CsrMatrix, Symmetry}; +use super::{CscMatrix, CsrMatrix, Sym}; use crate::StrError; use std::ffi::OsStr; use std::fmt::Write; @@ -35,7 +35,7 @@ where // write header if !vismatrix { - if symmetry == Symmetry::No { + if symmetry == Sym::No { write!(&mut buffer, "%%MatrixMarket matrix coordinate real general\n").unwrap(); } else { write!(&mut buffer, "%%MatrixMarket matrix coordinate real symmetric\n").unwrap(); @@ -103,7 +103,7 @@ where // write header if !vismatrix { - if symmetry == Symmetry::No { + if symmetry == Sym::No { write!(&mut buffer, "%%MatrixMarket matrix coordinate real general\n").unwrap(); } else { write!(&mut buffer, "%%MatrixMarket matrix coordinate real symmetric\n").unwrap(); diff --git a/russell_sparse/tests/test_complex_coo_matrix.rs b/russell_sparse/tests/test_complex_coo_matrix.rs index b10abf12..6943e2f1 100644 --- a/russell_sparse/tests/test_complex_coo_matrix.rs +++ b/russell_sparse/tests/test_complex_coo_matrix.rs @@ -5,7 +5,7 @@ use russell_sparse::StrError; #[test] fn test_complex_coo_matrix() -> Result<(), StrError> { - let sym = Some(Symmetry::new_general_lower()); + let sym = Some(Sym::new_general_lower()); let mut coo = ComplexCooMatrix::new(3, 3, 4, sym)?; coo.put(0, 0, cpx!(1.0, 0.1))?; coo.put(1, 0, cpx!(2.0, 0.2))?; From 36b3ce11fe8ff0934792db2fccc36d274123ff03 Mon Sep 17 00:00:00 2001 From: Dorival Pedroso Date: Sat, 9 Mar 2024 11:54:25 +1000 Subject: [PATCH 2/3] [sparse] Rename sym flag --- russell_sparse/examples/doc_coo_from_arrays.rs | 4 ++-- russell_sparse/examples/doc_csc_from_arrays.rs | 4 ++-- russell_sparse/examples/doc_csr_from_arrays.rs | 4 ++-- russell_sparse/src/bin/solve_matrix_market.rs | 4 ++-- russell_sparse/src/complex_lin_solver.rs | 2 +- russell_sparse/src/complex_solver_mumps.rs | 6 +++--- russell_sparse/src/complex_solver_umfpack.rs | 6 +++--- russell_sparse/src/coo_matrix.rs | 18 +++++++++--------- russell_sparse/src/csc_matrix.rs | 13 ++++++------- russell_sparse/src/csr_matrix.rs | 4 ++-- russell_sparse/src/lin_solver.rs | 2 +- russell_sparse/src/read_matrix_market.rs | 8 ++++---- russell_sparse/src/solver_mumps.rs | 6 +++--- russell_sparse/src/solver_umfpack.rs | 6 +++--- russell_sparse/src/sparse_matrix.rs | 14 +++++++------- russell_sparse/src/write_matrix_market.rs | 8 ++++---- 16 files changed, 54 insertions(+), 55 deletions(-) diff --git a/russell_sparse/examples/doc_coo_from_arrays.rs b/russell_sparse/examples/doc_coo_from_arrays.rs index 7a468207..cbc86c65 100644 --- a/russell_sparse/examples/doc_coo_from_arrays.rs +++ b/russell_sparse/examples/doc_coo_from_arrays.rs @@ -15,8 +15,8 @@ fn main() -> Result<(), StrError> { let values = vec![ 1.0, /*dup*/ 1.0, 3.0, 3.0, -1.0, 4.0, 4.0, -3.0, 1.0, 2.0, 2.0, 6.0, 1.0, ]; - let symmetry = None; - let coo = CooMatrix::from(nrow, ncol, row_indices, col_indices, values, symmetry)?; + let sym = None; + let coo = CooMatrix::from(nrow, ncol, row_indices, col_indices, values, sym)?; // covert to dense let a = coo.as_dense(); diff --git a/russell_sparse/examples/doc_csc_from_arrays.rs b/russell_sparse/examples/doc_csc_from_arrays.rs index 8909495c..47c4c486 100644 --- a/russell_sparse/examples/doc_csc_from_arrays.rs +++ b/russell_sparse/examples/doc_csc_from_arrays.rs @@ -29,8 +29,8 @@ fn main() -> Result<(), StrError> { 6.0, 1.0, // j = 4, count = 10, 11, // 12 ]; - let symmetry = None; - let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, symmetry)?; + let sym = None; + let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, sym)?; // covert to dense let a = csc.as_dense(); diff --git a/russell_sparse/examples/doc_csr_from_arrays.rs b/russell_sparse/examples/doc_csr_from_arrays.rs index 8ffb5ce4..0dc5cba5 100644 --- a/russell_sparse/examples/doc_csr_from_arrays.rs +++ b/russell_sparse/examples/doc_csr_from_arrays.rs @@ -29,8 +29,8 @@ fn main() -> Result<(), StrError> { 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 // count = 12 ]; - let symmetry = None; - let csr = CsrMatrix::new(nrow, ncol, row_pointers, col_indices, values, symmetry)?; + let sym = None; + let csr = CsrMatrix::new(nrow, ncol, row_pointers, col_indices, values, sym)?; // covert to dense let a = csr.as_dense(); diff --git a/russell_sparse/src/bin/solve_matrix_market.rs b/russell_sparse/src/bin/solve_matrix_market.rs index 59b54c4d..f878bef3 100644 --- a/russell_sparse/src/bin/solve_matrix_market.rs +++ b/russell_sparse/src/bin/solve_matrix_market.rs @@ -105,12 +105,12 @@ fn main() -> Result<(), StrError> { let mut mat = SparseMatrix::from_coo(coo); // save information about the matrix - let (nrow, ncol, nnz, symmetry) = mat.get_info(); + let (nrow, ncol, nnz, sym) = mat.get_info(); stats.set_matrix_name_from_path(&opt.matrix_market_file); stats.matrix.nrow = nrow; stats.matrix.ncol = ncol; stats.matrix.nnz = nnz; - stats.matrix.symmetry = format!("{:?}", symmetry); + stats.matrix.symmetry = format!("{:?}", sym); // allocate and configure the solver let mut solver = LinSolver::new(genie)?; diff --git a/russell_sparse/src/complex_lin_solver.rs b/russell_sparse/src/complex_lin_solver.rs index 9a60fe3d..3f98f317 100644 --- a/russell_sparse/src/complex_lin_solver.rs +++ b/russell_sparse/src/complex_lin_solver.rs @@ -13,7 +13,7 @@ pub trait ComplexLinSolTrait { /// /// # Notes /// - /// 1. The structure of the matrix (nrow, ncol, nnz, symmetry) must be + /// 1. The structure of the matrix (nrow, ncol, nnz, sym) must be /// exactly the same among multiple calls to `factorize`. The values may differ /// from call to call, nonetheless. /// 2. The first call to `factorize` will define the structure which must be diff --git a/russell_sparse/src/complex_solver_mumps.rs b/russell_sparse/src/complex_solver_mumps.rs index 229e3d41..cef226ea 100644 --- a/russell_sparse/src/complex_solver_mumps.rs +++ b/russell_sparse/src/complex_solver_mumps.rs @@ -195,7 +195,7 @@ impl ComplexLinSolTrait for ComplexSolverMUMPS { /// /// # Notes /// - /// 1. The structure of the matrix (nrow, ncol, nnz, symmetry) must be + /// 1. The structure of the matrix (nrow, ncol, nnz, sym) must be /// exactly the same among multiple calls to `factorize`. The values may differ /// from call to call, nonetheless. /// 2. The first call to `factorize` will define the structure which must be @@ -362,8 +362,8 @@ impl ComplexLinSolTrait for ComplexSolverMUMPS { let coo = mat.get_coo()?; // check already factorized data - let (nrow, ncol, nnz, symmetry) = coo.get_info(); - if symmetry != self.initialized_symmetry { + let (nrow, ncol, nnz, sym) = coo.get_info(); + if sym != self.initialized_symmetry { return Err("solve must use the same matrix (symmetry differs)"); } if nrow != self.initialized_ndim || ncol != self.initialized_ndim { diff --git a/russell_sparse/src/complex_solver_umfpack.rs b/russell_sparse/src/complex_solver_umfpack.rs index 1730b737..cd1f1bb7 100644 --- a/russell_sparse/src/complex_solver_umfpack.rs +++ b/russell_sparse/src/complex_solver_umfpack.rs @@ -183,7 +183,7 @@ impl ComplexLinSolTrait for ComplexSolverUMFPACK { /// /// # Notes /// - /// 1. The structure of the matrix (nrow, ncol, nnz, symmetry) must be + /// 1. The structure of the matrix (nrow, ncol, nnz, sym) must be /// exactly the same among multiple calls to `factorize`. The values may differ /// from call to call, nonetheless. /// 2. The first call to `factorize` will define the structure which must be @@ -325,8 +325,8 @@ impl ComplexLinSolTrait for ComplexSolverUMFPACK { let csc = mat.get_csc()?; // check already factorized data - let (nrow, ncol, nnz, symmetry) = csc.get_info(); - if symmetry != self.initialized_symmetry { + let (nrow, ncol, nnz, sym) = csc.get_info(); + if sym != self.initialized_symmetry { return Err("solve must use the same matrix (symmetry differs)"); } if nrow != self.initialized_ndim || ncol != self.initialized_ndim { diff --git a/russell_sparse/src/coo_matrix.rs b/russell_sparse/src/coo_matrix.rs index bf352cf6..e07c8262 100644 --- a/russell_sparse/src/coo_matrix.rs +++ b/russell_sparse/src/coo_matrix.rs @@ -220,8 +220,8 @@ where /// let values = vec![ /// 1.0, /*dup*/ 1.0, 3.0, 3.0, -1.0, 4.0, 4.0, -3.0, 1.0, 2.0, 2.0, 6.0, 1.0, /// ]; - /// let symmetry = None; - /// let coo = CooMatrix::from(nrow, ncol, row_indices, col_indices, values, symmetry)?; + /// let sym = None; + /// let coo = CooMatrix::from(nrow, ncol, row_indices, col_indices, values, sym)?; /// /// // covert to dense /// let a = coo.as_dense(); @@ -371,11 +371,11 @@ where /// │ 0 0 3 │\n\ /// └ ┘"; /// coo.reset(); - /// let (nrow, ncol, nnz, symmetry) = coo.get_info(); + /// let (nrow, ncol, nnz, sym) = coo.get_info(); /// assert_eq!(nrow, 3); /// assert_eq!(ncol, 3); /// assert_eq!(nnz, 0); - /// assert_eq!(symmetry, Symmetry::No); + /// assert_eq!(sym, Sym::No); /// Ok(()) /// } /// ``` @@ -616,7 +616,7 @@ where /// Returns information about the dimensions and symmetry type /// - /// Returns `(nrow, ncol, nnz, symmetry)` + /// Returns `(nrow, ncol, nnz, sym)` /// /// # Example /// @@ -626,11 +626,11 @@ where /// /// fn main() -> Result<(), StrError> { /// let coo = CooMatrix::new(1, 2, 3, None)?; - /// let (nrow, ncol, nnz, symmetry) = coo.get_info(); + /// let (nrow, ncol, nnz, sym) = coo.get_info(); /// assert_eq!(nrow, 1); /// assert_eq!(ncol, 2); /// assert_eq!(nnz, 0); - /// assert_eq!(symmetry, Symmetry::No); + /// assert_eq!(sym, Sym::No); /// Ok(()) /// } /// ``` @@ -752,11 +752,11 @@ mod tests { #[test] fn get_info_works() { let coo = NumCooMatrix::::new(1, 2, 10, None).unwrap(); - let (nrow, ncol, nnz, symmetry) = coo.get_info(); + let (nrow, ncol, nnz, sym) = coo.get_info(); assert_eq!(nrow, 1); assert_eq!(ncol, 2); assert_eq!(nnz, 0); - assert_eq!(symmetry, Sym::No); + assert_eq!(sym, Sym::No); } #[test] diff --git a/russell_sparse/src/csc_matrix.rs b/russell_sparse/src/csc_matrix.rs index cc9d7204..a80264f1 100644 --- a/russell_sparse/src/csc_matrix.rs +++ b/russell_sparse/src/csc_matrix.rs @@ -173,8 +173,8 @@ where /// 6.0, 1.0, // j = 4, count = 10, 11, /// // 12 /// ]; - /// let symmetry = None; - /// let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, symmetry)?; + /// let sym = None; + /// let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, sym)?; /// /// // covert to dense /// let a = csc.as_dense(); @@ -611,9 +611,8 @@ where /// 6.0, 1.0, // j = 4, count = 10, 11, /// // 12 /// ]; - /// let symmetry = None; - /// let csc = CscMatrix::new(nrow, ncol, - /// col_pointers, row_indices, values, symmetry)?; + /// let sym = None; + /// let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, sym)?; /// /// // covert to dense /// let a = csc.as_dense(); @@ -767,11 +766,11 @@ where /// let values = vec![10.0, 20.0]; /// let csc = CscMatrix::new(1, 2, /// col_pointers, row_indices, values, None)?; - /// let (nrow, ncol, nnz, symmetry) = csc.get_info(); + /// let (nrow, ncol, nnz, sym) = csc.get_info(); /// assert_eq!(nrow, 1); /// assert_eq!(ncol, 2); /// assert_eq!(nnz, 2); - /// assert_eq!(symmetry, Symmetry::No); + /// assert_eq!(sym, Sym::No); /// let a = csc.as_dense(); /// let correct = "┌ ┐\n\ /// │ 10 20 │\n\ diff --git a/russell_sparse/src/csr_matrix.rs b/russell_sparse/src/csr_matrix.rs index f247710f..805b654f 100644 --- a/russell_sparse/src/csr_matrix.rs +++ b/russell_sparse/src/csr_matrix.rs @@ -741,11 +741,11 @@ where /// let values = vec![10.0, 20.0]; /// let csr = CsrMatrix::new(1, 2, /// row_pointers, col_indices, values, None)?; - /// let (nrow, ncol, nnz, symmetry) = csr.get_info(); + /// let (nrow, ncol, nnz, sym) = csr.get_info(); /// assert_eq!(nrow, 1); /// assert_eq!(ncol, 2); /// assert_eq!(nnz, 2); - /// assert_eq!(symmetry, Symmetry::No); + /// assert_eq!(sym, Sym::No); /// let a = csr.as_dense(); /// let correct = "┌ ┐\n\ /// │ 10 20 │\n\ diff --git a/russell_sparse/src/lin_solver.rs b/russell_sparse/src/lin_solver.rs index 69e14238..60a66fac 100644 --- a/russell_sparse/src/lin_solver.rs +++ b/russell_sparse/src/lin_solver.rs @@ -13,7 +13,7 @@ pub trait LinSolTrait { /// /// # Notes /// - /// 1. The structure of the matrix (nrow, ncol, nnz, symmetry) must be + /// 1. The structure of the matrix (nrow, ncol, nnz, sym) must be /// exactly the same among multiple calls to `factorize`. The values may differ /// from call to call, nonetheless. /// 2. The first call to `factorize` will define the structure which must be diff --git a/russell_sparse/src/read_matrix_market.rs b/russell_sparse/src/read_matrix_market.rs index b7eebfef..4659f33a 100644 --- a/russell_sparse/src/read_matrix_market.rs +++ b/russell_sparse/src/read_matrix_market.rs @@ -265,11 +265,11 @@ impl MatrixMarketData { /// fn main() -> Result<(), StrError> { /// let name = "./data/matrix_market/ok_simple_general.mtx"; /// let coo = read_matrix_market(name, MMsymOption::LeaveAsLower)?; -/// let (nrow, ncol, nnz, symmetry) = coo.get_info(); +/// let (nrow, ncol, nnz, sym) = coo.get_info(); /// assert_eq!(nrow, 3); /// assert_eq!(ncol, 3); /// assert_eq!(nnz, 5); -/// assert_eq!(symmetry, Symmetry::No); +/// assert_eq!(sym, Sym::No); /// let a = coo.as_dense(); /// let correct = "┌ ┐\n\ /// │ 1 2 0 │\n\ @@ -304,11 +304,11 @@ impl MatrixMarketData { /// fn main() -> Result<(), StrError> { /// let name = "./data/matrix_market/ok_simple_symmetric.mtx"; /// let coo = read_matrix_market(name, MMsymOption::LeaveAsLower)?; -/// let (nrow, ncol, nnz, symmetry) = coo.get_info(); +/// let (nrow, ncol, nnz, sym) = coo.get_info(); /// assert_eq!(nrow, 3); /// assert_eq!(ncol, 3); /// assert_eq!(nnz, 4); -/// assert_eq!(symmetry, Symmetry::General(Storage::Lower)); +/// assert_eq!(sym, Sym::General(Storage::Lower)); /// let a = coo.as_dense(); /// let correct = "┌ ┐\n\ /// │ 1 2 0 │\n\ diff --git a/russell_sparse/src/solver_mumps.rs b/russell_sparse/src/solver_mumps.rs index 63980402..d92df2e6 100644 --- a/russell_sparse/src/solver_mumps.rs +++ b/russell_sparse/src/solver_mumps.rs @@ -232,7 +232,7 @@ impl LinSolTrait for SolverMUMPS { /// /// # Notes /// - /// 1. The structure of the matrix (nrow, ncol, nnz, symmetry) must be + /// 1. The structure of the matrix (nrow, ncol, nnz, sym) must be /// exactly the same among multiple calls to `factorize`. The values may differ /// from call to call, nonetheless. /// 2. The first call to `factorize` will define the structure which must be @@ -392,8 +392,8 @@ impl LinSolTrait for SolverMUMPS { let coo = mat.get_coo()?; // check already factorized data - let (nrow, ncol, nnz, symmetry) = coo.get_info(); - if symmetry != self.initialized_symmetry { + let (nrow, ncol, nnz, sym) = coo.get_info(); + if sym != self.initialized_symmetry { return Err("solve must use the same matrix (symmetry differs)"); } if nrow != self.initialized_ndim || ncol != self.initialized_ndim { diff --git a/russell_sparse/src/solver_umfpack.rs b/russell_sparse/src/solver_umfpack.rs index 194c4cdc..a003244d 100644 --- a/russell_sparse/src/solver_umfpack.rs +++ b/russell_sparse/src/solver_umfpack.rs @@ -169,7 +169,7 @@ impl LinSolTrait for SolverUMFPACK { /// /// # Notes /// - /// 1. The structure of the matrix (nrow, ncol, nnz, symmetry) must be + /// 1. The structure of the matrix (nrow, ncol, nnz, sym) must be /// exactly the same among multiple calls to `factorize`. The values may differ /// from call to call, nonetheless. /// 2. The first call to `factorize` will define the structure which must be @@ -304,8 +304,8 @@ impl LinSolTrait for SolverUMFPACK { let csc = mat.get_csc()?; // check already factorized data - let (nrow, ncol, nnz, symmetry) = csc.get_info(); - if symmetry != self.initialized_symmetry { + let (nrow, ncol, nnz, sym) = csc.get_info(); + if sym != self.initialized_symmetry { return Err("solve must use the same matrix (symmetry differs)"); } if nrow != self.initialized_ndim || ncol != self.initialized_ndim { diff --git a/russell_sparse/src/sparse_matrix.rs b/russell_sparse/src/sparse_matrix.rs index 4b9821fa..e25fe539 100644 --- a/russell_sparse/src/sparse_matrix.rs +++ b/russell_sparse/src/sparse_matrix.rs @@ -58,10 +58,10 @@ where /// * `ncol` -- (≥ 1) Is the number of columns of the sparse matrix (must be fit i32) /// * `max_nnz` -- (≥ 1) Maximum number of entries ≥ nnz (number of non-zeros), /// including entries with repeated indices. (must be fit i32) - /// * `symmetry` -- Defines the symmetry/storage, if any - pub fn new_coo(nrow: usize, ncol: usize, max_nnz: usize, symmetry: Option) -> Result { + /// * `sym` -- Defines the symmetry/storage, if any + pub fn new_coo(nrow: usize, ncol: usize, max_nnz: usize, sym: Option) -> Result { Ok(NumSparseMatrix { - coo: Some(NumCooMatrix::new(nrow, ncol, max_nnz, symmetry)?), + coo: Some(NumCooMatrix::new(nrow, ncol, max_nnz, sym)?), csc: None, csr: None, }) @@ -190,7 +190,7 @@ where /// Returns information about the dimensions and symmetry type /// - /// Returns `(nrow, ncol, nnz, symmetry)` + /// Returns `(nrow, ncol, nnz, sym)` /// /// **Priority**: CSC -> CSR -> COO pub fn get_info(&self) -> (usize, usize, usize, Sym) { @@ -652,7 +652,7 @@ mod tests { #[test] fn derive_methods_work() { let (coo, _, _, _) = Samples::tiny_1x1(); - let (nrow, ncol, nnz, symmetry) = coo.get_info(); + let (nrow, ncol, nnz, sym) = coo.get_info(); let mat = NumSparseMatrix::::from_coo(coo); let mut clone = mat.clone(); clone.get_coo_mut().unwrap().values[0] *= 2.0; @@ -665,8 +665,8 @@ mod tests { r#"{"coo":{"symmetry":"No","nrow":1,"ncol":1,"nnz":1,"max_nnz":1,"indices_i":[0],"indices_j":[0],"values":[123.0]},"csc":null,"csr":null}"# ); let from_json: NumSparseMatrix = serde_json::from_str(&json).unwrap(); - let (json_nrow, json_ncol, json_nnz, json_symmetry) = from_json.get_coo().unwrap().get_info(); - assert_eq!(json_symmetry, symmetry); + let (json_nrow, json_ncol, json_nnz, json_sym) = from_json.get_coo().unwrap().get_info(); + assert_eq!(json_sym, sym); assert_eq!(json_nrow, nrow); assert_eq!(json_ncol, ncol); assert_eq!(json_nnz, nnz); diff --git a/russell_sparse/src/write_matrix_market.rs b/russell_sparse/src/write_matrix_market.rs index 7448417c..36c98a07 100644 --- a/russell_sparse/src/write_matrix_market.rs +++ b/russell_sparse/src/write_matrix_market.rs @@ -31,11 +31,11 @@ where let d = if vismatrix { 0 } else { 1 }; // info - let (nrow, ncol, nnz, symmetry) = mat.get_info(); + let (nrow, ncol, nnz, sym) = mat.get_info(); // write header if !vismatrix { - if symmetry == Sym::No { + if sym == Sym::No { write!(&mut buffer, "%%MatrixMarket matrix coordinate real general\n").unwrap(); } else { write!(&mut buffer, "%%MatrixMarket matrix coordinate real symmetric\n").unwrap(); @@ -99,11 +99,11 @@ where let d = if vismatrix { 0 } else { 1 }; // info - let (nrow, ncol, nnz, symmetry) = mat.get_info(); + let (nrow, ncol, nnz, sym) = mat.get_info(); // write header if !vismatrix { - if symmetry == Sym::No { + if sym == Sym::No { write!(&mut buffer, "%%MatrixMarket matrix coordinate real general\n").unwrap(); } else { write!(&mut buffer, "%%MatrixMarket matrix coordinate real symmetric\n").unwrap(); From 88d63dbea5dbe6c6492db6dfb9e31891bb54ceac Mon Sep 17 00:00:00 2001 From: Dorival Pedroso Date: Sat, 9 Mar 2024 18:49:51 +1000 Subject: [PATCH 3/3] [IMPORTANT] Simplify the Sym enum --- README.md | 2 +- russell_ode/src/euler_backward.rs | 3 +- russell_ode/src/radau5.rs | 162 ++++------ russell_ode/src/samples.rs | 55 ++-- russell_ode/src/system.rs | 27 +- russell_sparse/README.md | 2 +- russell_sparse/examples/complex_system.rs | 2 +- .../examples/doc_coo_from_arrays.rs | 3 +- .../examples/doc_coo_new_put_reset.rs | 2 +- .../examples/doc_csc_from_arrays.rs | 3 +- russell_sparse/examples/doc_csc_from_coo.rs | 2 +- .../examples/doc_csr_from_arrays.rs | 3 +- russell_sparse/examples/doc_csr_from_coo.rs | 2 +- .../examples/doc_lin_solver_compute.rs | 2 +- .../examples/doc_lin_solver_umfpack_tiny.rs | 2 +- .../examples/doc_umfpack_quickstart_coo.rs | 2 +- .../examples/doc_umfpack_quickstart_csc.rs | 2 +- russell_sparse/examples/doc_umfpack_tiny.rs | 2 +- russell_sparse/examples/mumps_solve_small.rs | 2 +- .../examples/nonlinear_system_4eqs.rs | 2 +- russell_sparse/src/bin/mem_check.rs | 2 +- russell_sparse/src/bin/solve_matrix_market.rs | 4 +- russell_sparse/src/complex_coo_matrix.rs | 44 ++- russell_sparse/src/complex_solver_mumps.rs | 86 +++--- russell_sparse/src/complex_solver_umfpack.rs | 70 +++-- russell_sparse/src/coo_matrix.rs | 205 ++++++------- russell_sparse/src/csc_matrix.rs | 82 +++-- russell_sparse/src/csr_matrix.rs | 83 +++-- russell_sparse/src/enums.rs | 286 +++--------------- russell_sparse/src/lib.rs | 6 +- russell_sparse/src/lin_sol_params.rs | 4 + russell_sparse/src/lin_solver.rs | 2 +- russell_sparse/src/read_matrix_market.rs | 54 ++-- russell_sparse/src/samples.rs | 50 +-- russell_sparse/src/solver_mumps.rs | 108 ++++--- russell_sparse/src/solver_umfpack.rs | 94 +++--- russell_sparse/src/sparse_matrix.rs | 45 +-- russell_sparse/src/verify_lin_sys.rs | 10 +- russell_sparse/src/write_matrix_market.rs | 4 +- .../tests/test_complex_coo_matrix.rs | 3 +- russell_sparse/tests/test_complex_mumps.rs | 2 +- russell_sparse/tests/test_complex_umfpack.rs | 2 +- russell_sparse/tests/test_mumps.rs | 2 +- russell_sparse/tests/test_nonlinear_system.rs | 4 +- russell_sparse/tests/test_umfpack.rs | 2 +- 45 files changed, 642 insertions(+), 894 deletions(-) diff --git a/README.md b/README.md index bfe3248b..a8afa0c9 100644 --- a/README.md +++ b/README.md @@ -294,7 +294,7 @@ fn main() -> Result<(), StrError> { // . -1 -3 2 . // . . 1 . . // . 4 2 . 1 - let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?; + let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(1, 0, 3.0)?; diff --git a/russell_ode/src/euler_backward.rs b/russell_ode/src/euler_backward.rs index 0a6eb880..d0653a35 100644 --- a/russell_ode/src/euler_backward.rs +++ b/russell_ode/src/euler_backward.rs @@ -50,7 +50,6 @@ where system.jac_nnz }; let nnz = jac_nnz + ndim; // +ndim corresponds to the diagonal I matrix - let symmetry = Some(system.jac_symmetry); EulerBackward { params, system, @@ -58,7 +57,7 @@ where w: Vector::new(ndim), r: Vector::new(ndim), dy: Vector::new(ndim), - kk: SparseMatrix::new_coo(ndim, ndim, nnz, symmetry).unwrap(), + kk: SparseMatrix::new_coo(ndim, ndim, nnz, system.jac_sym).unwrap(), solver: LinSolver::new(params.newton.genie).unwrap(), } } diff --git a/russell_ode/src/radau5.rs b/russell_ode/src/radau5.rs index 71ba11cc..1b05cb37 100644 --- a/russell_ode/src/radau5.rs +++ b/russell_ode/src/radau5.rs @@ -110,7 +110,6 @@ where /// Allocates a new instance pub fn new(params: Params, system: &'a System) -> Self { let ndim = system.ndim; - let symmetry = Some(system.jac_symmetry); let mass_nnz = match system.mass_matrix.as_ref() { Some(mass) => mass.get_info().2, None => ndim, @@ -125,9 +124,9 @@ where Radau5 { params, system, - jj: SparseMatrix::new_coo(ndim, ndim, jac_nnz, symmetry).unwrap(), - kk_real: SparseMatrix::new_coo(ndim, ndim, nnz, symmetry).unwrap(), - kk_comp: ComplexSparseMatrix::new_coo(ndim, ndim, nnz, symmetry).unwrap(), + jj: SparseMatrix::new_coo(ndim, ndim, jac_nnz, system.jac_sym).unwrap(), + kk_real: SparseMatrix::new_coo(ndim, ndim, nnz, system.jac_sym).unwrap(), + kk_comp: ComplexSparseMatrix::new_coo(ndim, ndim, nnz, system.jac_sym).unwrap(), solver_real: LinSolver::new(params.newton.genie).unwrap(), solver_comp: ComplexLinSolver::new(params.newton.genie).unwrap(), reuse_jacobian: false, @@ -847,109 +846,62 @@ mod tests { assert_eq!(work.bench.n_function, n_fcn_correct); } - #[test] - fn radau5_works_mass_matrix() { - // problem - let (system, data, mut args) = Samples::simple_system_with_mass_matrix(false); - let yfx = data.y_analytical.unwrap(); - let ndim = system.ndim; - - // allocate structs - let params = Params::new(Method::Radau5); - let mut solver = Radau5::new(params, &system); - let mut work = Workspace::new(Method::Radau5); - - // message - println!("{:>4}{:>10}{:>10}{:>10}", "step", "err_y0", "err_y1", "err_y2"); - - // numerical approximation - let h = 0.1; - let mut x = data.x0; - let mut y = data.y0.clone(); - let mut y_ana = Vector::new(ndim); - for n in 0..5 { - // call step - solver.step(&mut work, x, &y, h, &mut args).unwrap(); - - // important: update n_accepted (must precede `accept`) - work.bench.n_accepted += 1; - - // call accept - solver.accept(&mut work, &mut x, &mut y, h, &mut args).unwrap(); - - // important: save previous stepsize and relative error (must succeed `accept`) - work.h_prev = h; - work.rel_error_prev = f64::max(params.step.rel_error_prev_min, work.rel_error); - - // check the results - yfx(&mut y_ana, x); - let err_y0 = f64::abs(y[0] - y_ana[0]); - let err_y1 = f64::abs(y[1] - y_ana[1]); - let err_y2 = f64::abs(y[2] - y_ana[2]); - println!( - "{:>4}{}{}{}", - n, - format_scientific(err_y0, 10, 2), - format_scientific(err_y1, 10, 2), - format_scientific(err_y2, 10, 2) - ); - assert!(err_y0 < 1e-9); - assert!(err_y1 < 1e-9); - assert!(err_y2 < 1e-8); - } - } - #[test] #[serial] - fn radau5_works_mass_matrix_symmetric_mumps() { - // problem - let (system, data, mut args) = Samples::simple_system_with_mass_matrix(true); - let yfx = data.y_analytical.unwrap(); - let ndim = system.ndim; - - // allocate structs - let mut params = Params::new(Method::Radau5); - params.newton.genie = Genie::Mumps; - let mut solver = Radau5::new(params, &system); - let mut work = Workspace::new(Method::Radau5); - - // message - println!("{:>4}{:>10}{:>10}{:>10}", "step", "err_y0", "err_y1", "err_y2"); - - // numerical approximation - let h = 0.1; - let mut x = data.x0; - let mut y = data.y0.clone(); - let mut y_ana = Vector::new(ndim); - for n in 0..5 { - // call step - solver.step(&mut work, x, &y, h, &mut args).unwrap(); - - // important: update n_accepted (must precede `accept`) - work.bench.n_accepted += 1; - - // call accept - solver.accept(&mut work, &mut x, &mut y, h, &mut args).unwrap(); - - // important: save previous stepsize and relative error (must succeed `accept`) - work.h_prev = h; - work.rel_error_prev = f64::max(params.step.rel_error_prev_min, work.rel_error); - - // check the results - yfx(&mut y_ana, x); - let err_y0 = f64::abs(y[0] - y_ana[0]); - let err_y1 = f64::abs(y[1] - y_ana[1]); - let err_y2 = f64::abs(y[2] - y_ana[2]); - println!( - "{:>4}{}{}{}", - n, - format_scientific(err_y0, 10, 2), - format_scientific(err_y1, 10, 2), - format_scientific(err_y2, 10, 2) - ); - assert!(err_y0 < 1e-9); - assert!(err_y1 < 1e-9); - assert!(err_y2 < 1e-8); + fn radau5_works_mass_matrix() { + for symmetric in [true, false] { + for genie in [Genie::Umfpack, Genie::Mumps] { + // problem + let (system, data, mut args) = Samples::simple_system_with_mass_matrix(symmetric, genie); + let yfx = data.y_analytical.unwrap(); + let ndim = system.ndim; + + // allocate structs + let mut params = Params::new(Method::Radau5); + params.newton.genie = genie; + let mut solver = Radau5::new(params, &system); + let mut work = Workspace::new(Method::Radau5); + + // message + println!("\nsymmetric = {:?} --- {:?}", symmetric, genie); + println!("{:>4}{:>10}{:>10}{:>10}", "step", "err_y0", "err_y1", "err_y2"); + + // numerical approximation + let h = 0.1; + let mut x = data.x0; + let mut y = data.y0.clone(); + let mut y_ana = Vector::new(ndim); + for n in 0..4 { + // call step + solver.step(&mut work, x, &y, h, &mut args).unwrap(); + + // important: update n_accepted (must precede `accept`) + work.bench.n_accepted += 1; + + // call accept + solver.accept(&mut work, &mut x, &mut y, h, &mut args).unwrap(); + + // important: save previous stepsize and relative error (must succeed `accept`) + work.h_prev = h; + work.rel_error_prev = f64::max(params.step.rel_error_prev_min, work.rel_error); + + // check the results + yfx(&mut y_ana, x); + let err_y0 = f64::abs(y[0] - y_ana[0]); + let err_y1 = f64::abs(y[1] - y_ana[1]); + let err_y2 = f64::abs(y[2] - y_ana[2]); + println!( + "{:>4}{}{}{}", + n, + format_scientific(err_y0, 10, 2), + format_scientific(err_y1, 10, 2), + format_scientific(err_y2, 10, 2) + ); + assert!(err_y0 < 1e-9); + assert!(err_y1 < 1e-9); + assert!(err_y2 < 1e-8); + } + } } } } diff --git a/russell_ode/src/samples.rs b/russell_ode/src/samples.rs index 09344475..be96ced5 100644 --- a/russell_ode/src/samples.rs +++ b/russell_ode/src/samples.rs @@ -2,7 +2,7 @@ use crate::StrError; use crate::{HasJacobian, NoArgs, System}; use russell_lab::math::PI; use russell_lab::Vector; -use russell_sparse::{CooMatrix, Sym}; +use russell_sparse::{CooMatrix, Genie, Sym}; /// Holds data corresponding to a sample ODE problem pub struct SampleData { @@ -139,7 +139,9 @@ impl Samples { /// /// # Input /// - /// * `lower_triangle` -- Considers the symmetry of the Jacobian and Mass matrices, and generates a lower triangular representation + /// * `symmetric` -- considers the symmetry of the Jacobian and Mass matrices + /// * `genie` -- if symmetric, this information is required to decide on the lower-triangle/full + /// representation which is consistent with the linear solver employed /// /// # Output /// @@ -152,7 +154,8 @@ impl Samples { /// * `data: SampleData` -- holds the initial values /// * `args: NoArgs` -- is a placeholder variable with the arguments to F and J pub fn simple_system_with_mass_matrix( - lower_triangle: bool, + symmetric: bool, + genie: Genie, ) -> ( System< impl Fn(&mut Vector, f64, &Vector, &mut NoArgs) -> Result<(), StrError>, @@ -162,12 +165,9 @@ impl Samples { SampleData, NoArgs, ) { - // selected symmetry option (for both Mass and Jacobian matrices) - let symmetry = if lower_triangle { - Some(Sym::new_general_lower()) - } else { - None - }; + // selected symmetric option (for both Mass and Jacobian matrices) + let sym = genie.symmetry(symmetric); + let triangular = sym.triangular(); // initial values let x0 = 0.0; @@ -176,7 +176,7 @@ impl Samples { // ODE system let ndim = 3; - let jac_nnz = if lower_triangle { 3 } else { 4 }; + let jac_nnz = if triangular { 3 } else { 4 }; let mut system = System::new( ndim, move |f: &mut Vector, x: f64, y: &Vector, _args: &mut NoArgs| { @@ -188,7 +188,7 @@ impl Samples { move |jj: &mut CooMatrix, _x: f64, _y: &Vector, m: f64, _args: &mut NoArgs| { jj.reset(); jj.put(0, 0, m * (-1.0)).unwrap(); - if !lower_triangle { + if !triangular { jj.put(0, 1, m * (1.0)).unwrap(); } jj.put(1, 0, m * (1.0)).unwrap(); @@ -197,14 +197,14 @@ impl Samples { }, HasJacobian::Yes, Some(jac_nnz), - symmetry, + if sym == Sym::No { None } else { Some(sym) }, ); // mass matrix - let mass_nnz = if lower_triangle { 4 } else { 5 }; + let mass_nnz = if triangular { 4 } else { 5 }; system.init_mass_matrix(mass_nnz).unwrap(); system.mass_put(0, 0, 1.0).unwrap(); - if !lower_triangle { + if !triangular { system.mass_put(0, 1, 1.0).unwrap(); } system.mass_put(1, 0, 1.0).unwrap(); @@ -843,8 +843,7 @@ mod tests { } // compute the analytical Jacobian matrix - let symmetry = Some(system.jac_symmetry); - let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap(); + let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap(); // compute the numerical Jacobian matrix @@ -871,8 +870,7 @@ mod tests { } // compute the analytical Jacobian matrix - let symmetry = Some(system.jac_symmetry); - let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap(); + let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap(); // compute the numerical Jacobian matrix @@ -899,8 +897,7 @@ mod tests { } // compute the analytical Jacobian matrix - let symmetry = Some(system.jac_symmetry); - let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap(); + let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap(); // compute the numerical Jacobian matrix @@ -927,8 +924,7 @@ mod tests { } // compute the analytical Jacobian matrix - let symmetry = Some(system.jac_symmetry); - let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap(); + let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap(); // compute the numerical Jacobian matrix @@ -947,8 +943,7 @@ mod tests { let (system, data, mut args) = Samples::robertson(); // compute the analytical Jacobian matrix - let symmetry = Some(system.jac_symmetry); - let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap(); + let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap(); // compute the numerical Jacobian matrix @@ -967,8 +962,7 @@ mod tests { let (system, data, mut args) = Samples::van_der_pol(None, false); // compute the analytical Jacobian matrix - let symmetry = Some(system.jac_symmetry); - let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap(); + let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap(); // compute the numerical Jacobian matrix @@ -987,8 +981,7 @@ mod tests { let (system, data, mut args) = Samples::van_der_pol(None, true); // compute the analytical Jacobian matrix - let symmetry = Some(system.jac_symmetry); - let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap(); + let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap(); // compute the numerical Jacobian matrix @@ -1007,8 +1000,7 @@ mod tests { let (system, data, mut args) = Samples::arenstorf(); // compute the analytical Jacobian matrix - let symmetry = Some(system.jac_symmetry); - let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap(); + let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap(); // compute the numerical Jacobian matrix @@ -1027,8 +1019,7 @@ mod tests { let (system, data, mut args) = Samples::amplifier1t(); // compute the analytical Jacobian matrix - let symmetry = Some(system.jac_symmetry); - let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap(); + let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, system.jac_sym).unwrap(); (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap(); // compute the numerical Jacobian matrix diff --git a/russell_ode/src/system.rs b/russell_ode/src/system.rs index 078e8cbd..96e119cd 100644 --- a/russell_ode/src/system.rs +++ b/russell_ode/src/system.rs @@ -71,8 +71,8 @@ where /// Number of non-zeros in the Jacobian matrix pub(crate) jac_nnz: usize, - /// Symmetry properties of the Jacobian matrix - pub(crate) jac_symmetry: Sym, + /// Symmetric flag for the Jacobian and mass matrices + pub(crate) jac_sym: Sym, /// Holds the mass matrix pub(crate) mass_matrix: Option, @@ -94,8 +94,8 @@ where /// * `function` -- implements the function: `dy/dx = f(x, y)` /// * `jacobian` -- implements the Jacobian: `J = df/dy` /// * `has_jacobian` -- indicates that the analytical Jacobian is available (input by `jacobian`) - /// * `jac_nnz` -- the number of non-zeros in the Jacobian; use None to indicate a full matrix (i.e., nnz = ndim * ndim) - /// * `jac_symmetry` -- specifies the type of symmetry representation for the Jacobian matrix + /// * `jac_nnz` -- the number of non-zeros in the Jacobian; use None to indicate a dense matrix (i.e., nnz = ndim * ndim) + /// * `jac_sym` -- specifies the symmetric flag for the Jacobian and mass matrices /// /// # Generics /// @@ -106,7 +106,7 @@ where jacobian: J, has_ana_jacobian: HasJacobian, jac_nnz: Option, - jac_symmetry: Option, + jac_sym: Option, ) -> Self { let jac_available = match has_ana_jacobian { HasJacobian::Yes => true, @@ -117,8 +117,8 @@ where function, jacobian, jac_available, - jac_nnz: if let Some(n) = jac_nnz { n } else { ndim * ndim }, - jac_symmetry: if let Some(s) = jac_symmetry { s } else { Sym::No }, + jac_nnz: if let Some(nnz) = jac_nnz { nnz } else { ndim * ndim }, + jac_sym: if let Some(sym) = jac_sym { sym } else { Sym::No }, mass_matrix: None, phantom: PhantomData, } @@ -132,12 +132,7 @@ where /// /// * `max_nnz` -- Max number of non-zero values pub fn init_mass_matrix(&mut self, max_nnz: usize) -> Result<(), StrError> { - let sym = if self.jac_symmetry == Sym::No { - None - } else { - Some(self.jac_symmetry) - }; - self.mass_matrix = Some(CooMatrix::new(self.ndim, self.ndim, max_nnz, sym).unwrap()); + self.mass_matrix = Some(CooMatrix::new(self.ndim, self.ndim, max_nnz, self.jac_sym).unwrap()); Ok(()) } @@ -219,7 +214,7 @@ mod tests { use super::{no_jacobian, System}; use crate::HasJacobian; use russell_lab::Vector; - use russell_sparse::CooMatrix; + use russell_sparse::{CooMatrix, Sym}; #[test] fn ode_system_most_none_works() { @@ -251,7 +246,7 @@ mod tests { let mut k = Vector::new(2); (ode.function)(&mut k, x, &y, &mut args).unwrap(); // call jacobian function - let mut jj = CooMatrix::new(2, 2, 2, None).unwrap(); + let mut jj = CooMatrix::new(2, 2, 2, Sym::No).unwrap(); let m = 1.0; assert_eq!( (ode.jacobian)(&mut jj, x, &y, m, &mut args), @@ -310,7 +305,7 @@ mod tests { let mut k = Vector::new(2); (ode.function)(&mut k, x, &y, &mut args).unwrap(); // call jacobian function - let mut jj = CooMatrix::new(2, 2, 2, None).unwrap(); + let mut jj = CooMatrix::new(2, 2, 2, Sym::No).unwrap(); let m = 1.0; (ode.jacobian)(&mut jj, x, &y, m, &mut args).unwrap(); // check diff --git a/russell_sparse/README.md b/russell_sparse/README.md index b691a95a..ff485dec 100644 --- a/russell_sparse/README.md +++ b/russell_sparse/README.md @@ -75,7 +75,7 @@ fn main() -> Result<(), StrError> { let mut umfpack = SolverUMFPACK::new()?; // allocate the coefficient matrix - let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?; + let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; coo.put(0, 0, 0.2)?; coo.put(0, 1, 0.2)?; coo.put(1, 0, 0.5)?; diff --git a/russell_sparse/examples/complex_system.rs b/russell_sparse/examples/complex_system.rs index 79385a0f..fbbc8053 100644 --- a/russell_sparse/examples/complex_system.rs +++ b/russell_sparse/examples/complex_system.rs @@ -43,7 +43,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, None)?; + let mut coo = ComplexSparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; // first column coo.put(0, 0, cpx!(19.73, 0.00))?; diff --git a/russell_sparse/examples/doc_coo_from_arrays.rs b/russell_sparse/examples/doc_coo_from_arrays.rs index cbc86c65..225264b6 100644 --- a/russell_sparse/examples/doc_coo_from_arrays.rs +++ b/russell_sparse/examples/doc_coo_from_arrays.rs @@ -15,8 +15,7 @@ fn main() -> Result<(), StrError> { let values = vec![ 1.0, /*dup*/ 1.0, 3.0, 3.0, -1.0, 4.0, 4.0, -3.0, 1.0, 2.0, 2.0, 6.0, 1.0, ]; - let sym = None; - let coo = CooMatrix::from(nrow, ncol, row_indices, col_indices, values, sym)?; + let coo = CooMatrix::from(nrow, ncol, row_indices, col_indices, values, Sym::No)?; // covert to dense let a = coo.as_dense(); diff --git a/russell_sparse/examples/doc_coo_new_put_reset.rs b/russell_sparse/examples/doc_coo_new_put_reset.rs index 272582bf..db97cbbe 100644 --- a/russell_sparse/examples/doc_coo_new_put_reset.rs +++ b/russell_sparse/examples/doc_coo_new_put_reset.rs @@ -8,7 +8,7 @@ fn main() -> Result<(), StrError> { // . -1 -3 2 . // . . 1 . . // . 4 2 . 1 - let mut coo = CooMatrix::new(5, 5, 13, None)?; + let mut coo = CooMatrix::new(5, 5, 13, Sym::No)?; coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(1, 0, 3.0)?; diff --git a/russell_sparse/examples/doc_csc_from_arrays.rs b/russell_sparse/examples/doc_csc_from_arrays.rs index 47c4c486..e7d80cb5 100644 --- a/russell_sparse/examples/doc_csc_from_arrays.rs +++ b/russell_sparse/examples/doc_csc_from_arrays.rs @@ -29,8 +29,7 @@ fn main() -> Result<(), StrError> { 6.0, 1.0, // j = 4, count = 10, 11, // 12 ]; - let sym = None; - let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, sym)?; + let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, Sym::No)?; // covert to dense let a = csc.as_dense(); diff --git a/russell_sparse/examples/doc_csc_from_coo.rs b/russell_sparse/examples/doc_csc_from_coo.rs index 6911072c..1b709542 100644 --- a/russell_sparse/examples/doc_csc_from_coo.rs +++ b/russell_sparse/examples/doc_csc_from_coo.rs @@ -9,7 +9,7 @@ fn main() -> Result<(), StrError> { // . . 1 . . // . 4 2 . 1 let (nrow, ncol, nnz) = (5, 5, 13); - let mut coo = CooMatrix::new(nrow, ncol, nnz, None)?; + let mut coo = CooMatrix::new(nrow, ncol, nnz, Sym::No)?; coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(1, 0, 3.0)?; diff --git a/russell_sparse/examples/doc_csr_from_arrays.rs b/russell_sparse/examples/doc_csr_from_arrays.rs index 0dc5cba5..3b5d084e 100644 --- a/russell_sparse/examples/doc_csr_from_arrays.rs +++ b/russell_sparse/examples/doc_csr_from_arrays.rs @@ -29,8 +29,7 @@ fn main() -> Result<(), StrError> { 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 // count = 12 ]; - let sym = None; - let csr = CsrMatrix::new(nrow, ncol, row_pointers, col_indices, values, sym)?; + let csr = CsrMatrix::new(nrow, ncol, row_pointers, col_indices, values, Sym::No)?; // covert to dense let a = csr.as_dense(); diff --git a/russell_sparse/examples/doc_csr_from_coo.rs b/russell_sparse/examples/doc_csr_from_coo.rs index deaac306..95a47b35 100644 --- a/russell_sparse/examples/doc_csr_from_coo.rs +++ b/russell_sparse/examples/doc_csr_from_coo.rs @@ -9,7 +9,7 @@ fn main() -> Result<(), StrError> { // . . 1 . . // . 4 2 . 1 let (nrow, ncol, nnz) = (5, 5, 13); - let mut coo = CooMatrix::new(nrow, ncol, nnz, None)?; + let mut coo = CooMatrix::new(nrow, ncol, nnz, Sym::No)?; coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(1, 0, 3.0)?; diff --git a/russell_sparse/examples/doc_lin_solver_compute.rs b/russell_sparse/examples/doc_lin_solver_compute.rs index 0e558a9a..43204ee6 100644 --- a/russell_sparse/examples/doc_lin_solver_compute.rs +++ b/russell_sparse/examples/doc_lin_solver_compute.rs @@ -8,7 +8,7 @@ fn main() -> Result<(), StrError> { let nnz = 5; // number of non-zero values // allocate the coefficient matrix - let mut mat = SparseMatrix::new_coo(ndim, ndim, nnz, None)?; + let mut mat = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; mat.put(0, 0, 0.2)?; mat.put(0, 1, 0.2)?; mat.put(1, 0, 0.5)?; diff --git a/russell_sparse/examples/doc_lin_solver_umfpack_tiny.rs b/russell_sparse/examples/doc_lin_solver_umfpack_tiny.rs index e9ecbdc8..5a36dc49 100644 --- a/russell_sparse/examples/doc_lin_solver_umfpack_tiny.rs +++ b/russell_sparse/examples/doc_lin_solver_umfpack_tiny.rs @@ -11,7 +11,7 @@ fn main() -> Result<(), StrError> { let mut solver = LinSolver::new(Genie::Umfpack)?; // allocate the coefficient matrix - let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?; + let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; coo.put(0, 0, 0.2)?; coo.put(0, 1, 0.2)?; coo.put(1, 0, 0.5)?; diff --git a/russell_sparse/examples/doc_umfpack_quickstart_coo.rs b/russell_sparse/examples/doc_umfpack_quickstart_coo.rs index d262b5f4..de1639f7 100644 --- a/russell_sparse/examples/doc_umfpack_quickstart_coo.rs +++ b/russell_sparse/examples/doc_umfpack_quickstart_coo.rs @@ -16,7 +16,7 @@ fn main() -> Result<(), StrError> { // . -1 -3 2 . // . . 1 . . // . 4 2 . 1 - let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?; + let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(1, 0, 3.0)?; diff --git a/russell_sparse/examples/doc_umfpack_quickstart_csc.rs b/russell_sparse/examples/doc_umfpack_quickstart_csc.rs index 73678ba0..07dd9aa7 100644 --- a/russell_sparse/examples/doc_umfpack_quickstart_csc.rs +++ b/russell_sparse/examples/doc_umfpack_quickstart_csc.rs @@ -9,7 +9,7 @@ fn main() -> Result<(), StrError> { let ap = vec![0, 2, 5, 9, 10, 12]; let ai = vec![0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4]; let ax = vec![2.0, 3.0, 3.0, -1.0, 4.0, 4.0, -3.0, 1.0, 2.0, 2.0, 6.0, 1.0]; - let mut csc = SparseMatrix::new_csc(n, n, ap, ai, ax, None)?; + let mut csc = SparseMatrix::new_csc(n, n, ap, ai, ax, Sym::No)?; // print the coefficient matrix let a = csc.as_dense(); diff --git a/russell_sparse/examples/doc_umfpack_tiny.rs b/russell_sparse/examples/doc_umfpack_tiny.rs index 9271822e..220fdbe5 100644 --- a/russell_sparse/examples/doc_umfpack_tiny.rs +++ b/russell_sparse/examples/doc_umfpack_tiny.rs @@ -11,7 +11,7 @@ fn main() -> Result<(), StrError> { let mut umfpack = SolverUMFPACK::new()?; // allocate the coefficient matrix - let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?; + let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; coo.put(0, 0, 0.2)?; coo.put(0, 1, 0.2)?; coo.put(1, 0, 0.5)?; diff --git a/russell_sparse/examples/mumps_solve_small.rs b/russell_sparse/examples/mumps_solve_small.rs index 90635492..6ef0d5e2 100644 --- a/russell_sparse/examples/mumps_solve_small.rs +++ b/russell_sparse/examples/mumps_solve_small.rs @@ -16,7 +16,7 @@ fn main() -> Result<(), StrError> { // . -1 -3 2 . // . . 1 . . // . 4 2 . 1 - let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?; + let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate coo.put(1, 0, 3.0)?; diff --git a/russell_sparse/examples/nonlinear_system_4eqs.rs b/russell_sparse/examples/nonlinear_system_4eqs.rs index 49c26e15..9610e60e 100644 --- a/russell_sparse/examples/nonlinear_system_4eqs.rs +++ b/russell_sparse/examples/nonlinear_system_4eqs.rs @@ -29,7 +29,7 @@ fn main() -> Result<(), StrError> { // allocate Jacobian matrix (J) as SparseMatrix let (neq, nnz) = (4, 16); - let mut jj = SparseMatrix::new_coo(neq, neq, nnz, None).unwrap(); + let mut jj = SparseMatrix::new_coo(neq, neq, nnz, Sym::No).unwrap(); // allocate residual (rr), vector of unknowns (uu), and minus-uu (mdu) let mut rr = Vector::new(neq); diff --git a/russell_sparse/src/bin/mem_check.rs b/russell_sparse/src/bin/mem_check.rs index 249522c5..4eac13ef 100644 --- a/russell_sparse/src/bin/mem_check.rs +++ b/russell_sparse/src/bin/mem_check.rs @@ -124,7 +124,7 @@ fn test_solver_singular(genie: Genie) { } }; - let mut coo_singular = match SparseMatrix::new_coo(ndim, ndim, nnz, None) { + let mut coo_singular = match SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No) { Ok(v) => v, Err(e) => { println!("FAIL(new CooMatrix): {}", e); diff --git a/russell_sparse/src/bin/solve_matrix_market.rs b/russell_sparse/src/bin/solve_matrix_market.rs index f878bef3..ed3e19a0 100644 --- a/russell_sparse/src/bin/solve_matrix_market.rs +++ b/russell_sparse/src/bin/solve_matrix_market.rs @@ -71,8 +71,8 @@ fn main() -> Result<(), StrError> { // select the symmetric handling option let handling = match genie { - Genie::Mumps => MMsymOption::LeaveAsLower, - Genie::Umfpack => MMsymOption::MakeItFull, + Genie::Mumps => MMsym::LeaveAsLower, + Genie::Umfpack => MMsym::MakeItFull, }; // configuration parameters diff --git a/russell_sparse/src/complex_coo_matrix.rs b/russell_sparse/src/complex_coo_matrix.rs index 34ea0482..7f17d716 100644 --- a/russell_sparse/src/complex_coo_matrix.rs +++ b/russell_sparse/src/complex_coo_matrix.rs @@ -30,8 +30,8 @@ impl ComplexCooMatrix { if other.ncol != self.ncol { return Err("matrices must have the same ncol"); } - if other.symmetry != self.symmetry { - return Err("matrices must have the same symmetry"); + if other.symmetric != self.symmetric { + return Err("matrices must have the same symmetric flag"); } self.reset(); for p in 0..other.nnz { @@ -68,8 +68,8 @@ impl ComplexCooMatrix { if other.ncol != self.ncol { return Err("matrices must have the same ncol"); } - if other.symmetry != self.symmetry { - return Err("matrices must have the same symmetry"); + if other.symmetric != self.symmetric { + return Err("matrices must have the same symmetric flag"); } for p in 0..other.nnz { let i = other.indices_i[p] as usize; @@ -84,20 +84,19 @@ impl ComplexCooMatrix { #[cfg(test)] mod tests { - use crate::{ComplexCooMatrix, CooMatrix, Storage, Sym}; + use crate::{ComplexCooMatrix, CooMatrix, Sym}; use num_complex::Complex64; use russell_lab::cpx; #[test] fn assign_capture_errors() { - let sym = Some(Sym::General(Storage::Full)); let nnz_a = 1; let nnz_b = 2; // wrong: must be ≤ nnz_a - let mut a_1x2 = ComplexCooMatrix::new(1, 2, nnz_a, None).unwrap(); - let b_2x1 = CooMatrix::new(2, 1, nnz_b, None).unwrap(); - let b_1x3 = CooMatrix::new(1, 3, nnz_b, None).unwrap(); - let b_1x2_sym = CooMatrix::new(1, 2, nnz_b, sym).unwrap(); - let mut b_1x2 = CooMatrix::new(1, 2, nnz_b, None).unwrap(); + let mut a_1x2 = ComplexCooMatrix::new(1, 2, nnz_a, Sym::No).unwrap(); + let b_2x1 = CooMatrix::new(2, 1, nnz_b, Sym::No).unwrap(); + let b_1x3 = CooMatrix::new(1, 3, nnz_b, Sym::No).unwrap(); + let b_1x2_sym = CooMatrix::new(1, 2, nnz_b, Sym::YesFull).unwrap(); + let mut b_1x2 = CooMatrix::new(1, 2, nnz_b, Sym::No).unwrap(); a_1x2.put(0, 0, cpx!(123.0, 321.0)).unwrap(); b_1x2.put(0, 0, 456.0).unwrap(); b_1x2.put(0, 1, 654.0).unwrap(); @@ -111,7 +110,7 @@ mod tests { ); assert_eq!( a_1x2.assign_real(2.0, 3.0, &b_1x2_sym).err(), - Some("matrices must have the same symmetry") + Some("matrices must have the same symmetric flag") ); assert_eq!( a_1x2.assign_real(2.0, 3.0, &b_1x2).err(), @@ -122,8 +121,8 @@ mod tests { #[test] fn assign_works() { let nnz = 2; - let mut a = ComplexCooMatrix::new(3, 2, nnz, None).unwrap(); - let mut b = CooMatrix::new(3, 2, nnz, None).unwrap(); + let mut a = ComplexCooMatrix::new(3, 2, nnz, Sym::No).unwrap(); + let mut b = CooMatrix::new(3, 2, nnz, Sym::No).unwrap(); a.put(2, 1, cpx!(1000.0, 2000.0)).unwrap(); b.put(0, 0, 10.0).unwrap(); b.put(2, 1, 20.0).unwrap(); @@ -148,14 +147,13 @@ mod tests { #[test] fn augment_capture_errors() { - let sym = Some(Sym::General(Storage::Full)); let nnz_a = 1; let nnz_b = 1; - let mut a_1x2 = ComplexCooMatrix::new(1, 2, nnz_a /* + nnz_b */, None).unwrap(); - let b_2x1 = CooMatrix::new(2, 1, nnz_b, None).unwrap(); - let b_1x3 = CooMatrix::new(1, 3, nnz_b, None).unwrap(); - let b_1x2_sym = CooMatrix::new(1, 2, nnz_b, sym).unwrap(); - let mut b_1x2 = CooMatrix::new(1, 2, nnz_b, None).unwrap(); + let mut a_1x2 = ComplexCooMatrix::new(1, 2, nnz_a /* + nnz_b */, Sym::No).unwrap(); + let b_2x1 = CooMatrix::new(2, 1, nnz_b, Sym::No).unwrap(); + let b_1x3 = CooMatrix::new(1, 3, nnz_b, Sym::No).unwrap(); + let b_1x2_sym = CooMatrix::new(1, 2, nnz_b, Sym::YesFull).unwrap(); + let mut b_1x2 = CooMatrix::new(1, 2, nnz_b, Sym::No).unwrap(); a_1x2.put(0, 0, cpx!(123.0, 321.0)).unwrap(); b_1x2.put(0, 0, 456.0).unwrap(); assert_eq!( @@ -168,7 +166,7 @@ mod tests { ); assert_eq!( a_1x2.augment_real(2.0, 3.0, &b_1x2_sym).err(), - Some("matrices must have the same symmetry") + Some("matrices must have the same symmetric flag") ); assert_eq!( a_1x2.augment_real(2.0, 3.0, &b_1x2).err(), @@ -180,8 +178,8 @@ mod tests { fn augment_works() { let nnz_a = 1; let nnz_b = 2; - let mut a = ComplexCooMatrix::new(3, 2, nnz_a + nnz_b, None).unwrap(); - let mut b = CooMatrix::new(3, 2, nnz_b, None).unwrap(); + let mut a = ComplexCooMatrix::new(3, 2, nnz_a + nnz_b, Sym::No).unwrap(); + let mut b = CooMatrix::new(3, 2, nnz_b, Sym::No).unwrap(); a.put(2, 1, cpx!(1000.0, 2000.0)).unwrap(); b.put(0, 0, 10.0).unwrap(); b.put(2, 1, 20.0).unwrap(); diff --git a/russell_sparse/src/complex_solver_mumps.rs b/russell_sparse/src/complex_solver_mumps.rs index cef226ea..75eb347f 100644 --- a/russell_sparse/src/complex_solver_mumps.rs +++ b/russell_sparse/src/complex_solver_mumps.rs @@ -80,8 +80,8 @@ pub struct ComplexSolverMUMPS { /// Indicates whether the sparse matrix has been factorized or not factorized: bool, - /// Holds the symmetry type used in the initialize - initialized_symmetry: Sym, + /// Holds the symmetric flag saved in initialize + initialized_sym: Sym, /// Holds the matrix dimension saved in initialize initialized_ndim: usize, @@ -161,7 +161,7 @@ impl ComplexSolverMUMPS { solver, initialized: false, factorized: false, - initialized_symmetry: Sym::No, + initialized_sym: Sym::No, initialized_ndim: 0, initialized_nnz: 0, effective_ordering: -1, @@ -190,7 +190,7 @@ impl ComplexLinSolTrait for ComplexSolverMUMPS { /// /// * `mat` -- the coefficient matrix A (one-base **COO** only, not CSC and not CSR). /// Also, the matrix must be square (`nrow = ncol`) and, if symmetric, - /// the symmetry/storage must [crate::Storage::Lower]. + /// the symmetric flag must be [Sym::YesLower] /// * `params` -- configuration parameters; None => use default /// /// # Notes @@ -202,24 +202,16 @@ impl ComplexLinSolTrait for ComplexSolverMUMPS { /// 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. - /// 4. For symmetric matrices, `MUMPS` requires that the symmetry/storage be [crate::Storage::Lower]. + /// 4. For symmetric matrices, `MUMPS` requires [Sym::YesLower]. /// 5. The COO matrix must be one-based. fn factorize(&mut self, mat: &mut ComplexSparseMatrix, params: Option) -> Result<(), StrError> { // get COO matrix let coo = mat.get_coo()?; - // check the COO matrix - if coo.nrow != coo.ncol { - return Err("the COO matrix must be square"); - } - if coo.nnz < 1 { - return Err("the COO matrix must have at least one non-zero value"); - } - - // check already initialized data + // check if self.initialized { - if coo.symmetry != self.initialized_symmetry { - return Err("subsequent factorizations must use the same matrix (symmetry differs)"); + if coo.symmetric != self.initialized_sym { + return Err("subsequent factorizations must use the same matrix (symmetric differs)"); } if coo.nrow != self.initialized_ndim { return Err("subsequent factorizations must use the same matrix (ndim differs)"); @@ -228,7 +220,16 @@ impl ComplexLinSolTrait for ComplexSolverMUMPS { return Err("subsequent factorizations must use the same matrix (nnz differs)"); } } else { - self.initialized_symmetry = coo.symmetry; + if coo.nrow != coo.ncol { + return Err("the COO matrix must be square"); + } + if coo.nnz < 1 { + return Err("the COO matrix must have at least one non-zero value"); + } + if coo.symmetric == Sym::YesFull || coo.symmetric == Sym::YesUpper { + return Err("MUMPS requires Sym::YesLower for symmetric matrices"); + } + self.initialized_sym = coo.symmetric; self.initialized_ndim = coo.nrow; self.initialized_nnz = coo.nnz; self.fortran_indices_i = vec![0; coo.nnz]; @@ -267,10 +268,9 @@ impl ComplexLinSolTrait for ComplexSolverMUMPS { let compute_determinant = if par.compute_determinant { 1 } else { 0 }; let verbose = if par.verbose { 1 } else { 0 }; - // extract the symmetry flags and check the storage type - let (general_symmetric, positive_definite) = coo.symmetry.status(true, false)?; - // matrix config + let general_symmetric = if coo.symmetric == Sym::YesLower { 1 } else { 0 }; + let positive_definite = if par.positive_definite { 1 } else { 0 }; let ndim = to_i32(coo.nrow); let nnz = to_i32(coo.nnz); @@ -341,7 +341,7 @@ impl ComplexLinSolTrait for ComplexSolverMUMPS { /// /// # Input /// - /// * `mat` -- the coefficient matrix A; must be square and, if symmetric, [crate::Storage::Lower]. + /// * `mat` -- the coefficient matrix A; must be square and, if symmetric, [Sym::YesLower]. /// * `rhs` -- the right-hand side vector with know values an dimension equal to mat.nrow /// * `verbose` -- shows messages /// @@ -363,8 +363,8 @@ impl ComplexLinSolTrait for ComplexSolverMUMPS { // check already factorized data let (nrow, ncol, nnz, sym) = coo.get_info(); - if sym != self.initialized_symmetry { - return Err("solve must use the same matrix (symmetry differs)"); + if sym != self.initialized_sym { + return Err("solve must use the same matrix (symmetric differs)"); } if nrow != self.initialized_ndim || ncol != self.initialized_ndim { return Err("solve must use the same matrix (ndim differs)"); @@ -455,7 +455,7 @@ impl ComplexLinSolTrait for ComplexSolverMUMPS { #[cfg(test)] mod tests { use super::*; - use crate::{ComplexCooMatrix, ComplexSparseMatrix, LinSolParams, Ordering, Samples, Scaling, Storage, Sym}; + use crate::{ComplexCooMatrix, ComplexSparseMatrix, LinSolParams, Ordering, Samples, Scaling, Sym}; use num_complex::Complex64; use russell_lab::{complex_approx_eq, complex_vec_approx_eq, cpx, ComplexVector}; use serial_test::serial; @@ -479,45 +479,45 @@ mod tests { ); // check COO matrix - let mut coo = ComplexCooMatrix::new(2, 1, 1, None).unwrap(); + let mut coo = ComplexCooMatrix::new(2, 1, 1, Sym::No).unwrap(); coo.put(0, 0, cpx!(4.0, 4.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), Some("the COO matrix must be square") ); - let coo = ComplexCooMatrix::new(1, 1, 1, None).unwrap(); + let coo = ComplexCooMatrix::new(1, 1, 1, Sym::No).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), Some("the COO matrix must have at least one non-zero value") ); - let mut coo = ComplexCooMatrix::new(1, 1, 1, Some(Sym::General(Storage::Full))).unwrap(); + let mut coo = ComplexCooMatrix::new(1, 1, 1, Sym::YesFull).unwrap(); coo.put(0, 0, cpx!(4.0, 4.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), - Some("if the matrix is general symmetric, the required storage is lower triangular") + Some("MUMPS requires Sym::YesLower for symmetric matrices") ); // check already factorized data - let mut coo = ComplexCooMatrix::new(2, 2, 2, None).unwrap(); + let mut coo = ComplexCooMatrix::new(2, 2, 2, Sym::No).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); coo.put(1, 1, cpx!(2.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); // ... factorize once => OK solver.factorize(&mut mat, None).unwrap(); - // ... change matrix (symmetry) - let mut coo = ComplexCooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); + // ... change matrix (symmetric) + let mut coo = ComplexCooMatrix::new(2, 2, 2, Sym::YesFull).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); coo.put(1, 1, cpx!(2.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), - Some("subsequent factorizations must use the same matrix (symmetry differs)") + Some("subsequent factorizations must use the same matrix (symmetric differs)") ); // ... change matrix (ndim) - let mut coo = ComplexCooMatrix::new(1, 1, 1, None).unwrap(); + let mut coo = ComplexCooMatrix::new(1, 1, 1, Sym::No).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( @@ -525,7 +525,7 @@ mod tests { Some("subsequent factorizations must use the same matrix (ndim differs)") ); // ... change matrix (nnz) - let mut coo = ComplexCooMatrix::new(2, 2, 1, None).unwrap(); + let mut coo = ComplexCooMatrix::new(2, 2, 1, Sym::No).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( @@ -537,7 +537,7 @@ mod tests { #[test] #[serial] fn factorize_fails_on_singular_matrix() { - let mut mat_singular = ComplexSparseMatrix::new_coo(3, 3, 2, None).unwrap(); + let mut mat_singular = ComplexSparseMatrix::new_coo(3, 3, 2, Sym::No).unwrap(); mat_singular.put(0, 0, cpx!(1.0, 0.0)).unwrap(); mat_singular.put(1, 1, cpx!(1.0, 0.0)).unwrap(); let mut solver = ComplexSolverMUMPS::new().unwrap(); @@ -550,7 +550,7 @@ mod tests { #[test] #[serial] fn solve_handles_errors() { - let mut coo = ComplexCooMatrix::new(2, 2, 2, None).unwrap(); + let mut coo = ComplexCooMatrix::new(2, 2, 2, Sym::No).unwrap(); coo.put(0, 0, cpx!(123.0, 1.0)).unwrap(); coo.put(1, 1, cpx!(456.0, 2.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); @@ -574,19 +574,19 @@ mod tests { solver.solve(&mut x, &mut mat, &rhs, false), Err("the dimension of the right-hand side vector is incorrect") ); - // wrong symmetry + // wrong symmetric let rhs = ComplexVector::new(2); - let mut coo_wrong = ComplexCooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); + let mut coo_wrong = ComplexCooMatrix::new(2, 2, 2, Sym::YesFull).unwrap(); coo_wrong.put(0, 0, cpx!(123.0, 1.0)).unwrap(); coo_wrong.put(1, 1, cpx!(456.0, 2.0)).unwrap(); let mut mat_wrong = ComplexSparseMatrix::from_coo(coo_wrong); mat_wrong.get_csc_or_from_coo().unwrap(); // make sure to convert to CSC (because we're not calling factorize on this wrong matrix) assert_eq!( solver.solve(&mut x, &mut mat_wrong, &rhs, false), - Err("solve must use the same matrix (symmetry differs)") + Err("solve must use the same matrix (symmetric differs)") ); // wrong ndim - let mut coo_wrong = ComplexCooMatrix::new(1, 1, 1, None).unwrap(); + let mut coo_wrong = ComplexCooMatrix::new(1, 1, 1, Sym::No).unwrap(); coo_wrong.put(0, 0, cpx!(123.0, 1.0)).unwrap(); let mut mat_wrong = ComplexSparseMatrix::from_coo(coo_wrong); mat_wrong.get_csc_or_from_coo().unwrap(); // make sure to convert to CSC (because we're not calling factorize on this wrong matrix) @@ -595,7 +595,7 @@ mod tests { Err("solve must use the same matrix (ndim differs)") ); // wrong nnz - let mut coo_wrong = ComplexCooMatrix::new(2, 2, 3, None).unwrap(); + let mut coo_wrong = ComplexCooMatrix::new(2, 2, 3, Sym::No).unwrap(); coo_wrong.put(0, 0, cpx!(123.0, 1.0)).unwrap(); coo_wrong.put(1, 1, cpx!(456.0, 2.0)).unwrap(); coo_wrong.put(0, 1, cpx!(100.0, 1.0)).unwrap(); @@ -658,10 +658,10 @@ mod tests { complex_vec_approx_eq(x_again.as_data(), x_correct, 1e-14); // solve with positive-definite matrix works - let sym = Some(Sym::PositiveDefinite(Storage::Lower)); + params.positive_definite = true; let nrow = 5; let ncol = 5; - let mut coo_pd_lower = ComplexCooMatrix::new(nrow, ncol, 9, sym).unwrap(); + let mut coo_pd_lower = ComplexCooMatrix::new(nrow, ncol, 9, Sym::YesLower).unwrap(); coo_pd_lower.put(0, 0, cpx!(9.0, 0.0)).unwrap(); coo_pd_lower.put(1, 1, cpx!(0.5, 0.0)).unwrap(); coo_pd_lower.put(2, 2, cpx!(12.0, 0.0)).unwrap(); diff --git a/russell_sparse/src/complex_solver_umfpack.rs b/russell_sparse/src/complex_solver_umfpack.rs index cd1f1bb7..402a506b 100644 --- a/russell_sparse/src/complex_solver_umfpack.rs +++ b/russell_sparse/src/complex_solver_umfpack.rs @@ -82,8 +82,8 @@ pub struct ComplexSolverUMFPACK { /// Indicates whether the sparse matrix has been factorized or not factorized: bool, - /// Holds the symmetry type used in initialize - initialized_symmetry: Sym, + /// Holds the symmetric flag saved in initialize + initialized_sym: Sym, /// Holds the matrix dimension saved in initialize initialized_ndim: usize, @@ -152,7 +152,7 @@ impl ComplexSolverUMFPACK { solver, initialized: false, factorized: false, - initialized_symmetry: Sym::No, + initialized_sym: Sym::No, initialized_ndim: 0, initialized_nnz: 0, effective_strategy: -1, @@ -178,7 +178,7 @@ impl ComplexLinSolTrait for ComplexSolverUMFPACK { /// /// * `mat` -- the coefficient matrix A (**COO** or **CSC**, but not CSR). /// Also, the matrix must be square (`nrow = ncol`) and, if symmetric, - /// the symmetry/storage must [crate::Storage::Full]. + /// the symmetric flag must be [Sym::YesFull] /// * `params` -- configuration parameters; None => use default /// /// # Notes @@ -190,24 +190,16 @@ impl ComplexLinSolTrait for ComplexSolverUMFPACK { /// 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. - /// 4. For symmetric matrices, `UMFPACK` requires that the symmetry/storage be [crate::Storage::Full]. + /// 4. For symmetric matrices, `UMFPACK` requires [Sym::YesFull]. fn factorize(&mut self, mat: &mut ComplexSparseMatrix, params: Option) -> Result<(), StrError> { // get CSC matrix // (or convert from COO if CSC is not available and COO is available) let csc = mat.get_csc_or_from_coo()?; - // check CSC matrix - if csc.nrow != csc.ncol { - return Err("the matrix must be square"); - } - if csc.symmetry.triangular() { - return Err("for UMFPACK, the matrix must not be triangular"); - } - - // check already initialized data + // check if self.initialized { - if csc.symmetry != self.initialized_symmetry { - return Err("subsequent factorizations must use the same matrix (symmetry differs)"); + if csc.symmetric != self.initialized_sym { + return Err("subsequent factorizations must use the same matrix (symmetric differs)"); } if csc.nrow != self.initialized_ndim { return Err("subsequent factorizations must use the same matrix (ndim differs)"); @@ -216,7 +208,13 @@ impl ComplexLinSolTrait for ComplexSolverUMFPACK { return Err("subsequent factorizations must use the same matrix (nnz differs)"); } } else { - self.initialized_symmetry = csc.symmetry; + if csc.nrow != csc.ncol { + return Err("the matrix must be square"); + } + if csc.symmetric == Sym::YesLower || csc.symmetric == Sym::YesUpper { + return Err("UMFPACK requires Sym::YesFull for symmetric matrices"); + } + self.initialized_sym = csc.symmetric; self.initialized_ndim = csc.nrow; self.initialized_nnz = csc.col_pointers[csc.ncol] as usize; } @@ -303,7 +301,7 @@ impl ComplexLinSolTrait for ComplexSolverUMFPACK { /// /// # Input /// - /// * `mat` -- the coefficient matrix A; must be square and, if symmetric, [crate::Storage::Full]. + /// * `mat` -- the coefficient matrix A; must be square and, if symmetric, [Sym::YesFull]. /// * `rhs` -- the right-hand side vector with know values an dimension equal to mat.nrow /// * `verbose` -- shows messages /// @@ -326,8 +324,8 @@ impl ComplexLinSolTrait for ComplexSolverUMFPACK { // check already factorized data let (nrow, ncol, nnz, sym) = csc.get_info(); - if sym != self.initialized_symmetry { - return Err("solve must use the same matrix (symmetry differs)"); + if sym != self.initialized_sym { + return Err("solve must use the same matrix (symmetric differs)"); } if nrow != self.initialized_ndim || ncol != self.initialized_ndim { return Err("solve must use the same matrix (ndim differs)"); @@ -410,7 +408,7 @@ impl ComplexLinSolTrait for ComplexSolverUMFPACK { #[cfg(test)] mod tests { use super::*; - use crate::{ComplexCooMatrix, ComplexSparseMatrix, Ordering, Samples, Scaling, Storage}; + use crate::{ComplexCooMatrix, ComplexSparseMatrix, Ordering, Samples, Scaling, Sym}; use num_complex::Complex64; use russell_lab::{complex_approx_eq, complex_vec_approx_eq, cpx, ComplexVector}; @@ -427,7 +425,7 @@ mod tests { assert!(!solver.factorized); // COO to CSC errors - let coo = ComplexCooMatrix::new(1, 1, 1, None).unwrap(); + let coo = ComplexCooMatrix::new(1, 1, 1, Sym::No).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), @@ -445,27 +443,27 @@ mod tests { let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), - Some("for UMFPACK, the matrix must not be triangular") + Some("UMFPACK requires Sym::YesFull for symmetric matrices") ); // check already factorized data - let mut coo = ComplexCooMatrix::new(2, 2, 2, None).unwrap(); + let mut coo = ComplexCooMatrix::new(2, 2, 2, Sym::No).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); coo.put(1, 1, cpx!(2.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); // ... factorize once => OK solver.factorize(&mut mat, None).unwrap(); - // ... change matrix (symmetry) - let mut coo = ComplexCooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); + // ... change matrix (symmetric) + let mut coo = ComplexCooMatrix::new(2, 2, 2, Sym::YesFull).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); coo.put(1, 1, cpx!(2.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), - Some("subsequent factorizations must use the same matrix (symmetry differs)") + Some("subsequent factorizations must use the same matrix (symmetric differs)") ); // ... change matrix (ndim) - let mut coo = ComplexCooMatrix::new(1, 1, 1, None).unwrap(); + let mut coo = ComplexCooMatrix::new(1, 1, 1, Sym::No).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( @@ -473,7 +471,7 @@ mod tests { Some("subsequent factorizations must use the same matrix (ndim differs)") ); // ... change matrix (nnz) - let mut coo = ComplexCooMatrix::new(2, 2, 1, None).unwrap(); + let mut coo = ComplexCooMatrix::new(2, 2, 1, Sym::No).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); assert_eq!( @@ -514,7 +512,7 @@ mod tests { #[test] fn factorize_fails_on_singular_matrix() { let mut solver = ComplexSolverUMFPACK::new().unwrap(); - let mut coo = ComplexCooMatrix::new(2, 2, 2, None).unwrap(); + let mut coo = ComplexCooMatrix::new(2, 2, 2, Sym::No).unwrap(); coo.put(0, 0, cpx!(1.0, 0.0)).unwrap(); coo.put(1, 1, cpx!(0.0, 0.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); @@ -523,7 +521,7 @@ mod tests { #[test] fn solve_handles_errors() { - let mut coo = ComplexCooMatrix::new(2, 2, 2, None).unwrap(); + let mut coo = ComplexCooMatrix::new(2, 2, 2, Sym::No).unwrap(); coo.put(0, 0, cpx!(123.0, 1.0)).unwrap(); coo.put(1, 1, cpx!(456.0, 2.0)).unwrap(); let mut mat = ComplexSparseMatrix::from_coo(coo); @@ -547,19 +545,19 @@ mod tests { solver.solve(&mut x, &mut mat, &rhs, false), Err("the dimension of the right-hand side vector is incorrect") ); - // wrong symmetry + // wrong symmetric let rhs = ComplexVector::new(2); - let mut coo_wrong = ComplexCooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); + let mut coo_wrong = ComplexCooMatrix::new(2, 2, 2, Sym::YesFull).unwrap(); coo_wrong.put(0, 0, cpx!(123.0, 1.0)).unwrap(); coo_wrong.put(1, 1, cpx!(456.0, 2.0)).unwrap(); let mut mat_wrong = ComplexSparseMatrix::from_coo(coo_wrong); mat_wrong.get_csc_or_from_coo().unwrap(); // make sure to convert to CSC (because we're not calling factorize on this wrong matrix) assert_eq!( solver.solve(&mut x, &mut mat_wrong, &rhs, false), - Err("solve must use the same matrix (symmetry differs)") + Err("solve must use the same matrix (symmetric differs)") ); // wrong ndim - let mut coo_wrong = ComplexCooMatrix::new(1, 1, 1, None).unwrap(); + let mut coo_wrong = ComplexCooMatrix::new(1, 1, 1, Sym::No).unwrap(); coo_wrong.put(0, 0, cpx!(123.0, 1.0)).unwrap(); let mut mat_wrong = ComplexSparseMatrix::from_coo(coo_wrong); mat_wrong.get_csc_or_from_coo().unwrap(); // make sure to convert to CSC (because we're not calling factorize on this wrong matrix) @@ -568,7 +566,7 @@ mod tests { Err("solve must use the same matrix (ndim differs)") ); // wrong nnz - let mut coo_wrong = ComplexCooMatrix::new(2, 2, 3, None).unwrap(); + let mut coo_wrong = ComplexCooMatrix::new(2, 2, 3, Sym::No).unwrap(); coo_wrong.put(0, 0, cpx!(123.0, 1.0)).unwrap(); coo_wrong.put(1, 1, cpx!(456.0, 2.0)).unwrap(); coo_wrong.put(0, 1, cpx!(100.0, 1.0)).unwrap(); diff --git a/russell_sparse/src/coo_matrix.rs b/russell_sparse/src/coo_matrix.rs index e07c8262..9400bf0b 100644 --- a/russell_sparse/src/coo_matrix.rs +++ b/russell_sparse/src/coo_matrix.rs @@ -1,4 +1,4 @@ -use super::{Storage, Sym}; +use super::Sym; use crate::to_i32; use crate::StrError; use num_traits::{Num, NumCast}; @@ -22,8 +22,8 @@ pub struct NumCooMatrix where T: AddAssign + MulAssign + Num + NumCast + Copy + DeserializeOwned + Serialize, { - /// Defines the symmetry and storage: lower-triangular, upper-triangular, full-matrix - pub(crate) symmetry: Sym, + /// Indicates whether the matrix is symmetric or not. If symmetric, indicates the representation too. + pub(crate) symmetric: Sym, /// Holds the number of rows (must fit i32) pub(crate) nrow: usize, @@ -84,7 +84,8 @@ where /// * `ncol` -- (≥ 1) Is the number of columns of the sparse matrix (must be fit i32) /// * `max_nnz` -- (≥ 1) Maximum number of entries ≥ nnz (number of non-zeros), /// including entries with repeated indices. (must be fit i32) - /// * `symmetry` -- Defines the symmetry/storage, if any + /// * `symmetric` -- indicates whether the matrix is symmetric or not. + /// If symmetric, indicates the representation too. /// /// # Examples /// @@ -99,7 +100,7 @@ where /// // . -1 -3 2 . /// // . . 1 . . /// // . 4 2 . 1 - /// let mut coo = CooMatrix::new(5, 5, 13, None)?; + /// let mut coo = CooMatrix::new(5, 5, 13, Sym::No)?; /// coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate /// coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate /// coo.put(1, 0, 3.0)?; @@ -167,7 +168,7 @@ where /// Ok(()) /// } /// ``` - pub fn new(nrow: usize, ncol: usize, max_nnz: usize, symmetry: Option) -> Result { + pub fn new(nrow: usize, ncol: usize, max_nnz: usize, symmetric: Sym) -> Result { if nrow < 1 { return Err("nrow must be ≥ 1"); } @@ -178,7 +179,7 @@ where return Err("max_nnz must be ≥ 1"); } Ok(NumCooMatrix { - symmetry: if let Some(v) = symmetry { v } else { Sym::No }, + symmetric, nrow, ncol, nnz: 0, @@ -198,7 +199,7 @@ where /// * `row_indices` -- (len = nnz) Is the array of row indices /// * `col_indices` -- (len = nnz) Is the array of columns indices /// * `values` -- (len = nnz) Is the array of non-zero values - /// * `symmetry` -- Defines the symmetry/storage, if any + /// * `symmetric` -- Defines the symmetry/storage, if any /// /// # Examples /// @@ -220,8 +221,7 @@ where /// let values = vec![ /// 1.0, /*dup*/ 1.0, 3.0, 3.0, -1.0, 4.0, 4.0, -3.0, 1.0, 2.0, 2.0, 6.0, 1.0, /// ]; - /// let sym = None; - /// let coo = CooMatrix::from(nrow, ncol, row_indices, col_indices, values, sym)?; + /// let coo = CooMatrix::from(nrow, ncol, row_indices, col_indices, values, Sym::No)?; /// /// // covert to dense /// let a = coo.as_dense(); @@ -242,7 +242,7 @@ where row_indices: Vec, col_indices: Vec, values: Vec, - symmetry: Option, + symmetric: Sym, ) -> Result { if nrow < 1 { return Err("nrow must be ≥ 1"); @@ -271,7 +271,7 @@ where } } Ok(NumCooMatrix { - symmetry: if let Some(v) = symmetry { v } else { Sym::No }, + symmetric, nrow, ncol, nnz, @@ -299,7 +299,7 @@ where /// /// fn main() -> Result<(), StrError> { /// let (nrow, ncol, nnz) = (3, 3, 4); - /// let mut coo = CooMatrix::new(nrow, ncol, nnz, None)?; + /// let mut coo = CooMatrix::new(nrow, ncol, nnz, Sym::No)?; /// coo.put(0, 0, 1.0)?; /// coo.put(1, 1, 2.0)?; /// coo.put(2, 2, 3.0)?; @@ -325,15 +325,14 @@ where if self.nnz >= self.max_nnz { return Err("COO matrix: max number of items has been reached"); } - if self.symmetry != Sym::No { - if self.symmetry.lower() { - if j > i { - return Err("COO matrix: j > i is incorrect for lower triangular storage"); - } - } else if self.symmetry.upper() { - if j < i { - return Err("COO matrix: j < i is incorrect for upper triangular storage"); - } + if self.symmetric == Sym::YesLower { + if j > i { + return Err("COO matrix: j > i is incorrect for lower triangular storage"); + } + } + if self.symmetric == Sym::YesUpper { + if j < i { + return Err("COO matrix: j < i is incorrect for upper triangular storage"); } } @@ -359,7 +358,7 @@ where /// /// fn main() -> Result<(), StrError> { /// let (nrow, ncol, max_nnz) = (3, 3, 10); - /// let mut coo = CooMatrix::new(nrow, ncol, max_nnz, None)?; + /// let mut coo = CooMatrix::new(nrow, ncol, max_nnz, Sym::No)?; /// coo.put(0, 0, 1.0)?; /// coo.put(1, 1, 2.0)?; /// coo.put(2, 2, 3.0)?; @@ -393,7 +392,7 @@ where /// // define (4 x 4) sparse matrix with 6+1 non-zero values /// // (with an extra ij-repeated entry) /// let (nrow, ncol, max_nnz) = (4, 4, 10); - /// let mut coo = CooMatrix::new(nrow, ncol, max_nnz, None)?; + /// let mut coo = CooMatrix::new(nrow, ncol, max_nnz, Sym::No)?; /// coo.put(0, 0, 0.5)?; // (0, 0, a00/2) << duplicate /// coo.put(0, 0, 0.5)?; // (0, 0, a00/2) << duplicate /// coo.put(0, 1, 2.0)?; @@ -437,7 +436,7 @@ where /// // define (4 x 4) sparse matrix with 6+1 non-zero values /// // (with an extra ij-repeated entry) /// let (nrow, ncol, max_nnz) = (4, 4, 10); - /// let mut coo = CooMatrix::new(nrow, ncol, max_nnz, None)?; + /// let mut coo = CooMatrix::new(nrow, ncol, max_nnz, Sym::No)?; /// coo.put(0, 0, 0.5)?; // (0, 0, a00/2) << duplicate /// coo.put(0, 0, 0.5)?; // (0, 0, a00/2) << duplicate /// coo.put(0, 1, 2.0)?; @@ -464,7 +463,7 @@ where if m != self.nrow || n != self.ncol { return Err("wrong matrix dimensions"); } - let mirror_required = self.symmetry.triangular(); + let mirror_required = self.symmetric.triangular(); a.fill(T::zero()); for p in 0..self.nnz { let i = self.indices_i[p] as usize; @@ -506,7 +505,7 @@ where /// fn main() -> Result<(), StrError> { /// // set sparse matrix (3 x 3) with 6 non-zeros /// let (nrow, ncol, max_nnz) = (3, 3, 6); - /// let mut coo = CooMatrix::new(nrow, ncol, max_nnz, None)?; + /// let mut coo = CooMatrix::new(nrow, ncol, max_nnz, Sym::No)?; /// coo.put(0, 0, 1.0)?; /// coo.put(1, 0, 2.0)?; /// coo.put(1, 1, 3.0)?; @@ -545,7 +544,7 @@ where if v.dim() != self.nrow { return Err("v vector is incompatible"); } - let mirror_required = self.symmetry.triangular(); + let mirror_required = self.symmetric.triangular(); v.fill(T::zero()); for p in 0..self.nnz { let i = self.indices_i[p] as usize; @@ -575,7 +574,7 @@ where if other.ncol != self.ncol { return Err("matrices must have the same ncol"); } - if other.symmetry != self.symmetry { + if other.symmetric != self.symmetric { return Err("matrices must have the same symmetry"); } self.reset(); @@ -603,7 +602,7 @@ where if other.ncol != self.ncol { return Err("matrices must have the same ncol"); } - if other.symmetry != self.symmetry { + if other.symmetric != self.symmetric { return Err("matrices must have the same symmetry"); } for p in 0..other.nnz { @@ -625,7 +624,7 @@ where /// use russell_sparse::StrError; /// /// fn main() -> Result<(), StrError> { - /// let coo = CooMatrix::new(1, 2, 3, None)?; + /// let coo = CooMatrix::new(1, 2, 3, Sym::No)?; /// let (nrow, ncol, nnz, sym) = coo.get_info(); /// assert_eq!(nrow, 1); /// assert_eq!(ncol, 2); @@ -635,17 +634,7 @@ where /// } /// ``` pub fn get_info(&self) -> (usize, usize, usize, Sym) { - (self.nrow, self.ncol, self.nnz, self.symmetry) - } - - /// Returns the storage corresponding to the symmetry type (if any) - pub fn get_storage(&self) -> Storage { - Sym::storage(self.symmetry) - } - - /// Returns whether the symmetry flag corresponds to a symmetric matrix or not - pub fn get_symmetric(&self) -> bool { - self.symmetry != Sym::No + (self.nrow, self.ncol, self.nnz, self.symmetric) } /// Get an access to the row indices @@ -692,24 +681,30 @@ where #[cfg(test)] mod tests { use super::NumCooMatrix; - use crate::{Samples, Storage, Sym}; + use crate::{Samples, Sym}; use num_complex::Complex64; use russell_lab::{complex_vec_approx_eq, cpx, vec_approx_eq, ComplexVector, NumMatrix, NumVector}; #[test] fn new_captures_errors() { - assert_eq!(NumCooMatrix::::new(0, 1, 3, None).err(), Some("nrow must be ≥ 1")); - assert_eq!(NumCooMatrix::::new(1, 0, 3, None).err(), Some("ncol must be ≥ 1")); assert_eq!( - NumCooMatrix::::new(1, 1, 0, None).err(), + NumCooMatrix::::new(0, 1, 3, Sym::No).err(), + Some("nrow must be ≥ 1") + ); + assert_eq!( + NumCooMatrix::::new(1, 0, 3, Sym::No).err(), + Some("ncol must be ≥ 1") + ); + assert_eq!( + NumCooMatrix::::new(1, 1, 0, Sym::No).err(), Some("max_nnz must be ≥ 1") ); } #[test] fn new_works() { - let coo = NumCooMatrix::::new(1, 1, 3, None).unwrap(); - assert_eq!(coo.symmetry, Sym::No); + let coo = NumCooMatrix::::new(1, 1, 3, Sym::No).unwrap(); + assert_eq!(coo.symmetric, Sym::No); assert_eq!(coo.nrow, 1); assert_eq!(coo.ncol, 1); assert_eq!(coo.nnz, 0); @@ -722,21 +717,21 @@ mod tests { #[test] #[rustfmt::skip] fn from_captures_errors(){ - assert_eq!(NumCooMatrix::::from(0, 1, vec![ 0], vec![ 0], vec![0.0], None).err(), Some("nrow must be ≥ 1")); - assert_eq!(NumCooMatrix::::from(1, 0, vec![ 0], vec![ 0], vec![0.0], None).err(), Some("ncol must be ≥ 1")); - assert_eq!(NumCooMatrix::::from(1, 1, vec![ ], vec![ 0], vec![0.0], None).err(), Some("nnz must be ≥ 1")); - assert_eq!(NumCooMatrix::::from(1, 1, vec![ 0], vec![ ], vec![0.0], None).err(), Some("col_indices.len() must be = nnz")); - assert_eq!(NumCooMatrix::::from(1, 1, vec![ 0], vec![ 0], vec![ ], None).err(), Some("values.len() must be = nnz")); - assert_eq!(NumCooMatrix::::from(1, 1, vec![-1], vec![ 0], vec![0.0], None).err(), Some("row index is out-of-range")); - assert_eq!(NumCooMatrix::::from(1, 1, vec![ 1], vec![ 0], vec![0.0], None).err(), Some("row index is out-of-range")); - assert_eq!(NumCooMatrix::::from(1, 1, vec![ 0], vec![-1], vec![0.0], None).err(), Some("col index is out-of-range")); - assert_eq!(NumCooMatrix::::from(1, 1, vec![ 0], vec![ 1], vec![0.0], None).err(), Some("col index is out-of-range")); + assert_eq!(NumCooMatrix::::from(0, 1, vec![ 0], vec![ 0], vec![0.0], Sym::No).err(), Some("nrow must be ≥ 1")); + assert_eq!(NumCooMatrix::::from(1, 0, vec![ 0], vec![ 0], vec![0.0], Sym::No).err(), Some("ncol must be ≥ 1")); + assert_eq!(NumCooMatrix::::from(1, 1, vec![ ], vec![ 0], vec![0.0], Sym::No).err(), Some("nnz must be ≥ 1")); + assert_eq!(NumCooMatrix::::from(1, 1, vec![ 0], vec![ ], vec![0.0], Sym::No).err(), Some("col_indices.len() must be = nnz")); + assert_eq!(NumCooMatrix::::from(1, 1, vec![ 0], vec![ 0], vec![ ], Sym::No).err(), Some("values.len() must be = nnz")); + assert_eq!(NumCooMatrix::::from(1, 1, vec![-1], vec![ 0], vec![0.0], Sym::No).err(), Some("row index is out-of-range")); + assert_eq!(NumCooMatrix::::from(1, 1, vec![ 1], vec![ 0], vec![0.0], Sym::No).err(), Some("row index is out-of-range")); + assert_eq!(NumCooMatrix::::from(1, 1, vec![ 0], vec![-1], vec![0.0], Sym::No).err(), Some("col index is out-of-range")); + assert_eq!(NumCooMatrix::::from(1, 1, vec![ 0], vec![ 1], vec![0.0], Sym::No).err(), Some("col index is out-of-range")); } #[test] fn from_works() { - let coo = NumCooMatrix::::from(1, 1, vec![0], vec![0], vec![123.0], None).unwrap(); - assert_eq!(coo.symmetry, Sym::No); + let coo = NumCooMatrix::::from(1, 1, vec![0], vec![0], vec![123.0], Sym::No).unwrap(); + assert_eq!(coo.symmetric, Sym::No); assert_eq!(coo.nrow, 1); assert_eq!(coo.ncol, 1); assert_eq!(coo.nnz, 1); @@ -744,14 +739,13 @@ mod tests { assert_eq!(coo.indices_i, &[0]); assert_eq!(coo.indices_j, &[0]); assert_eq!(coo.values, &[123.0]); - let sym = Some(Sym::new_general_full()); - let coo = NumCooMatrix::::from(1, 1, vec![0], vec![0], vec![123.0], sym).unwrap(); - assert_eq!(coo.symmetry, Sym::General(Storage::Full)); + let coo = NumCooMatrix::::from(1, 1, vec![0], vec![0], vec![123.0], Sym::YesFull).unwrap(); + assert_eq!(coo.symmetric, Sym::YesFull); } #[test] fn get_info_works() { - let coo = NumCooMatrix::::new(1, 2, 10, None).unwrap(); + let coo = NumCooMatrix::::new(1, 2, 10, Sym::No).unwrap(); let (nrow, ncol, nnz, sym) = coo.get_info(); assert_eq!(nrow, 1); assert_eq!(ncol, 2); @@ -761,7 +755,7 @@ mod tests { #[test] fn put_fails_on_wrong_values() { - let mut coo = NumCooMatrix::::new(1, 1, 1, None).unwrap(); + let mut coo = NumCooMatrix::::new(1, 1, 1, Sym::No).unwrap(); assert_eq!( coo.put(1, 0, 0).err(), Some("COO matrix: index of row is outside range") @@ -775,14 +769,12 @@ mod tests { coo.put(0, 0, 0).err(), Some("COO matrix: max number of items has been reached") ); - let sym = Some(Sym::General(Storage::Lower)); - let mut coo = NumCooMatrix::::new(2, 2, 4, sym).unwrap(); + let mut coo = NumCooMatrix::::new(2, 2, 4, Sym::YesLower).unwrap(); assert_eq!( coo.put(0, 1, 0).err(), Some("COO matrix: j > i is incorrect for lower triangular storage") ); - let sym = Some(Sym::General(Storage::Upper)); - let mut coo = NumCooMatrix::::new(2, 2, 4, sym).unwrap(); + let mut coo = NumCooMatrix::::new(2, 2, 4, Sym::YesUpper).unwrap(); assert_eq!( coo.put(1, 0, 0).err(), Some("COO matrix: j < i is incorrect for upper triangular storage") @@ -791,7 +783,7 @@ mod tests { #[test] fn put_works() { - let mut coo = NumCooMatrix::::new(3, 3, 5, None).unwrap(); + let mut coo = NumCooMatrix::::new(3, 3, 5, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); assert_eq!(coo.nnz, 1); coo.put(0, 1, 2.0).unwrap(); @@ -806,7 +798,7 @@ mod tests { #[test] fn reset_works() { - let mut coo = NumCooMatrix::::new(2, 2, 4, None).unwrap(); + let mut coo = NumCooMatrix::::new(2, 2, 4, Sym::No).unwrap(); assert_eq!(coo.nnz, 0); coo.put(0, 0, 1.0).unwrap(); coo.put(0, 1, 4.0).unwrap(); @@ -819,7 +811,7 @@ mod tests { #[test] fn to_dense_fails_on_wrong_dims() { - let mut coo = NumCooMatrix::::new(1, 1, 1, None).unwrap(); + let mut coo = NumCooMatrix::::new(1, 1, 1, Sym::No).unwrap(); coo.put(0, 0, 123.0).unwrap(); let mut a_2x1 = NumMatrix::::new(2, 1); let mut a_1x2 = NumMatrix::::new(1, 2); @@ -829,7 +821,7 @@ mod tests { #[test] fn to_dense_works() { - let mut coo = NumCooMatrix::::new(3, 3, 5, None).unwrap(); + let mut coo = NumCooMatrix::::new(3, 3, 5, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(0, 1, 2.0).unwrap(); coo.put(1, 0, 3.0).unwrap(); @@ -854,11 +846,11 @@ mod tests { assert_eq!(bb.get(0, 0), 1.0); assert_eq!(bb.get(1, 0), 3.0); // empty matrix - let empty = NumCooMatrix::::new(2, 2, 3, None).unwrap(); + let empty = NumCooMatrix::::new(2, 2, 3, Sym::No).unwrap(); let mat = empty.as_dense(); assert_eq!(mat.as_data(), &[0.0, 0.0, 0.0, 0.0]); // single component matrix - let mut single = NumCooMatrix::::new(1, 1, 1, None).unwrap(); + let mut single = NumCooMatrix::::new(1, 1, 1, Sym::No).unwrap(); single.put(0, 0, 123.0).unwrap(); let mat = single.as_dense(); assert_eq!(mat.as_data(), &[123.0]); @@ -868,7 +860,7 @@ mod tests { fn to_dense_with_duplicates_works() { // allocate a square matrix let (nrow, ncol, nnz) = (5, 5, 13); - let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, None).unwrap(); + let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); // << (0, 0, a00/2) coo.put(0, 0, 1.0).unwrap(); // << (0, 0, a00/2) coo.put(1, 0, 3.0).unwrap(); @@ -898,8 +890,7 @@ mod tests { #[test] fn to_dense_symmetric_lower_works() { - let sym = Some(Sym::General(Storage::Lower)); - let mut coo = NumCooMatrix::::new(3, 3, 4, sym).unwrap(); + let mut coo = NumCooMatrix::::new(3, 3, 4, Sym::YesLower).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(1, 0, 2.0).unwrap(); coo.put(1, 1, 3.0).unwrap(); @@ -916,8 +907,7 @@ mod tests { #[test] fn to_dense_symmetric_upper_works() { - let sym = Some(Sym::General(Storage::Upper)); - let mut coo = NumCooMatrix::::new(3, 3, 4, sym).unwrap(); + let mut coo = NumCooMatrix::::new(3, 3, 4, Sym::YesUpper).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(0, 1, 2.0).unwrap(); coo.put(1, 1, 3.0).unwrap(); @@ -934,7 +924,7 @@ mod tests { #[test] fn mat_vec_mul_captures_errors() { - let mut coo = NumCooMatrix::::new(2, 2, 1, None).unwrap(); + let mut coo = NumCooMatrix::::new(2, 2, 1, Sym::No).unwrap(); coo.put(0, 0, 123).unwrap(); let u = NumVector::::new(3); let mut v = NumVector::::new(coo.nrow); @@ -949,7 +939,7 @@ mod tests { // 1.0 2.0 3.0 // 0.1 0.2 0.3 // 10.0 20.0 30.0 - let mut coo = NumCooMatrix::::new(3, 3, 9, None).unwrap(); + let mut coo = NumCooMatrix::::new(3, 3, 9, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(0, 1, 2.0).unwrap(); coo.put(0, 2, 3.0).unwrap(); @@ -970,7 +960,7 @@ mod tests { vec_approx_eq(v.as_data(), correct_v, 1e-15); // single component matrix - let mut single = NumCooMatrix::::new(1, 1, 1, None).unwrap(); + let mut single = NumCooMatrix::::new(1, 1, 1, Sym::No).unwrap(); single.put(0, 0, 123.0).unwrap(); let u = NumVector::from(&[2.0]); let mut v = NumVector::::new(1); @@ -986,8 +976,7 @@ mod tests { // 3 1 1 7 // 2 1 5 1 8 let (nrow, ncol, nnz) = (5, 5, 15); - let sym = Some(Sym::General(Storage::Lower)); - let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, sym).unwrap(); + let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, Sym::YesLower).unwrap(); coo.put(0, 0, 2.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); coo.put(2, 2, 9.0).unwrap(); @@ -1022,8 +1011,7 @@ mod tests { // 3 1 1 7 1 // 2 1 5 1 8 let (nrow, ncol, nnz) = (5, 5, 25); - let sym = Some(Sym::General(Storage::Full)); - let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, sym).unwrap(); + let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, Sym::YesFull).unwrap(); coo.put(0, 0, 2.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); coo.put(2, 2, 9.0).unwrap(); @@ -1066,8 +1054,7 @@ mod tests { // -1 2 -1 => -1 2 // -1 2 -1 2 let (nrow, ncol, nnz) = (3, 3, 5); - let sym = Some(Sym::PositiveDefinite(Storage::Lower)); - let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, sym).unwrap(); + let mut coo = NumCooMatrix::::new(nrow, ncol, nnz, Sym::YesLower).unwrap(); coo.put(0, 0, 2.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); coo.put(2, 2, 2.0).unwrap(); @@ -1104,14 +1091,13 @@ mod tests { #[test] fn assign_capture_errors() { - let sym = Some(Sym::General(Storage::Full)); let nnz_a = 1; let nnz_b = 2; // wrong: must be ≤ nnz_a - let mut a_1x2 = NumCooMatrix::::new(1, 2, nnz_a, None).unwrap(); - let b_2x1 = NumCooMatrix::::new(2, 1, nnz_b, None).unwrap(); - let b_1x3 = NumCooMatrix::::new(1, 3, nnz_b, None).unwrap(); - let b_1x2_sym = NumCooMatrix::::new(1, 2, nnz_b, sym).unwrap(); - let mut b_1x2 = NumCooMatrix::::new(1, 2, nnz_b, None).unwrap(); + let mut a_1x2 = NumCooMatrix::::new(1, 2, nnz_a, Sym::No).unwrap(); + let b_2x1 = NumCooMatrix::::new(2, 1, nnz_b, Sym::No).unwrap(); + let b_1x3 = NumCooMatrix::::new(1, 3, nnz_b, Sym::No).unwrap(); + let b_1x2_sym = NumCooMatrix::::new(1, 2, nnz_b, Sym::YesFull).unwrap(); + let mut b_1x2 = NumCooMatrix::::new(1, 2, nnz_b, Sym::No).unwrap(); a_1x2.put(0, 0, 123).unwrap(); b_1x2.put(0, 0, 456).unwrap(); b_1x2.put(0, 1, 654).unwrap(); @@ -1130,8 +1116,8 @@ mod tests { #[test] fn assign_works() { let nnz = 2; - let mut a = NumCooMatrix::::new(3, 2, nnz, None).unwrap(); - let mut b = NumCooMatrix::::new(3, 2, nnz, None).unwrap(); + let mut a = NumCooMatrix::::new(3, 2, nnz, Sym::No).unwrap(); + let mut b = NumCooMatrix::::new(3, 2, nnz, Sym::No).unwrap(); a.put(2, 1, 1000.0).unwrap(); b.put(0, 0, 10.0).unwrap(); b.put(2, 1, 20.0).unwrap(); @@ -1156,14 +1142,13 @@ mod tests { #[test] fn augment_capture_errors() { - let sym = Some(Sym::General(Storage::Full)); let nnz_a = 1; let nnz_b = 1; - let mut a_1x2 = NumCooMatrix::::new(1, 2, nnz_a /* + nnz_b */, None).unwrap(); - let b_2x1 = NumCooMatrix::::new(2, 1, nnz_b, None).unwrap(); - let b_1x3 = NumCooMatrix::::new(1, 3, nnz_b, None).unwrap(); - let b_1x2_sym = NumCooMatrix::::new(1, 2, nnz_b, sym).unwrap(); - let mut b_1x2 = NumCooMatrix::::new(1, 2, nnz_b, None).unwrap(); + let mut a_1x2 = NumCooMatrix::::new(1, 2, nnz_a /* + nnz_b */, Sym::No).unwrap(); + let b_2x1 = NumCooMatrix::::new(2, 1, nnz_b, Sym::No).unwrap(); + let b_1x3 = NumCooMatrix::::new(1, 3, nnz_b, Sym::No).unwrap(); + let b_1x2_sym = NumCooMatrix::::new(1, 2, nnz_b, Sym::YesFull).unwrap(); + let mut b_1x2 = NumCooMatrix::::new(1, 2, nnz_b, Sym::No).unwrap(); a_1x2.put(0, 0, 123).unwrap(); b_1x2.put(0, 0, 456).unwrap(); assert_eq!(a_1x2.augment(2, &b_2x1).err(), Some("matrices must have the same nrow")); @@ -1182,8 +1167,8 @@ mod tests { fn augment_works() { let nnz_a = 1; let nnz_b = 2; - let mut a = NumCooMatrix::::new(3, 2, nnz_a + nnz_b, None).unwrap(); - let mut b = NumCooMatrix::::new(3, 2, nnz_b, None).unwrap(); + let mut a = NumCooMatrix::::new(3, 2, nnz_a + nnz_b, Sym::No).unwrap(); + let mut b = NumCooMatrix::::new(3, 2, nnz_b, Sym::No).unwrap(); a.put(2, 1, 1000.0).unwrap(); b.put(0, 0, 10.0).unwrap(); b.put(2, 1, 20.0).unwrap(); @@ -1210,17 +1195,11 @@ mod tests { fn getters_are_correct() { let (coo, _, _, _) = Samples::rectangular_1x2(false, false); assert_eq!(coo.get_info(), (1, 2, 2, Sym::No)); - assert_eq!(coo.get_storage(), Storage::Full); - assert_eq!(coo.get_symmetric(), false); assert_eq!(coo.get_row_indices(), &[0, 0]); assert_eq!(coo.get_col_indices(), &[0, 1]); assert_eq!(coo.get_values(), &[10.0, 20.0]); - let sym = Some(Sym::new_general_full()); - let coo = NumCooMatrix::::new(2, 2, 2, sym).unwrap(); - assert_eq!(coo.get_symmetric(), true); - - let mut coo = NumCooMatrix::::new(2, 1, 2, None).unwrap(); + let mut coo = NumCooMatrix::::new(2, 1, 2, Sym::No).unwrap(); coo.put(0, 0, 123.0).unwrap(); coo.put(1, 0, 456.0).unwrap(); assert_eq!(coo.get_values_mut(), &[123.0, 456.0]); @@ -1240,10 +1219,10 @@ mod tests { let json = serde_json::to_string(&coo).unwrap(); assert_eq!( json, - r#"{"symmetry":"No","nrow":1,"ncol":1,"nnz":1,"max_nnz":1,"indices_i":[0],"indices_j":[0],"values":[123.0]}"# + r#"{"symmetric":"No","nrow":1,"ncol":1,"nnz":1,"max_nnz":1,"indices_i":[0],"indices_j":[0],"values":[123.0]}"# ); let from_json: NumCooMatrix = serde_json::from_str(&json).unwrap(); - assert_eq!(from_json.symmetry, coo.symmetry); + assert_eq!(from_json.symmetric, coo.symmetric); assert_eq!(from_json.nrow, coo.nrow); assert_eq!(from_json.ncol, coo.ncol); assert_eq!(from_json.nnz, coo.nnz); diff --git a/russell_sparse/src/csc_matrix.rs b/russell_sparse/src/csc_matrix.rs index a80264f1..cfa0631f 100644 --- a/russell_sparse/src/csc_matrix.rs +++ b/russell_sparse/src/csc_matrix.rs @@ -54,8 +54,8 @@ pub struct NumCscMatrix where T: AddAssign + MulAssign + Num + NumCast + Copy + DeserializeOwned + Serialize, { - /// Defines the symmetry and storage: lower-triangular, upper-triangular, full-matrix - pub(crate) symmetry: Sym, + /// Indicates whether the matrix is symmetric or not. If symmetric, indicates the representation too. + pub(crate) symmetric: Sym, /// Holds the number of rows (must fit i32) pub(crate) nrow: usize, @@ -125,6 +125,8 @@ where /// to the number of non-zero values (sorted) /// * `row_indices` -- (len = nnz) row indices (sorted) /// * `values` -- the non-zero components of the matrix + /// * `symmetric` -- indicates whether the matrix is symmetric or not. + /// If symmetric, indicates the representation too. /// /// The following conditions must be satisfied (nnz is the number of non-zeros /// and nnz_dup is the number of non-zeros with possible duplicates): @@ -173,8 +175,7 @@ where /// 6.0, 1.0, // j = 4, count = 10, 11, /// // 12 /// ]; - /// let sym = None; - /// let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, sym)?; + /// let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, Sym::No)?; /// /// // covert to dense /// let a = csc.as_dense(); @@ -195,7 +196,7 @@ where col_pointers: Vec, row_indices: Vec, values: Vec, - symmetry: Option, + symmetric: Sym, ) -> Result { if nrow < 1 { return Err("nrow must be ≥ 1"); @@ -242,7 +243,7 @@ where } } Ok(NumCscMatrix { - symmetry: if let Some(v) = symmetry { v } else { Sym::No }, + symmetric, nrow, ncol, col_pointers, @@ -275,7 +276,7 @@ where /// // . . 1 . . /// // . 4 2 . 1 /// let (nrow, ncol, nnz) = (5, 5, 13); - /// let mut coo = CooMatrix::new(nrow, ncol, nnz, None)?; + /// let mut coo = CooMatrix::new(nrow, ncol, nnz, Sym::No)?; /// coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate /// coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate /// coo.put(1, 0, 3.0)?; @@ -329,7 +330,7 @@ where return Err("COO to CSC requires nnz > 0"); } let mut csc = NumCscMatrix { - symmetry: coo.symmetry, + symmetric: coo.symmetric, nrow: coo.nrow, ncol: coo.ncol, col_pointers: vec![0; coo.ncol + 1], @@ -354,7 +355,7 @@ where /// may have been summed up. The final nnz is available as `nnz = col_pointers[ncol]`. pub fn update_from_coo(&mut self, coo: &NumCooMatrix) -> Result<(), StrError> { // check dimensions - if coo.symmetry != self.symmetry { + if coo.symmetric != self.symmetric { return Err("coo.symmetry must be equal to csc.symmetry"); } if coo.nrow != self.nrow { @@ -517,7 +518,7 @@ where // allocate the CSC arrays let mut csc = NumCscMatrix { - symmetry: csr.symmetry, + symmetric: csr.symmetric, nrow: csr.nrow, ncol: csr.ncol, col_pointers: vec![0; ncol + 1], @@ -611,8 +612,7 @@ where /// 6.0, 1.0, // j = 4, count = 10, 11, /// // 12 /// ]; - /// let sym = None; - /// let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, sym)?; + /// let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, Sym::No)?; /// /// // covert to dense /// let a = csc.as_dense(); @@ -675,9 +675,7 @@ where /// 6.0, 1.0, // j = 4, count = 10, 11, /// // 12 /// ]; - /// let symmetry = None; - /// let csc = CscMatrix::new(nrow, ncol, - /// col_pointers, row_indices, values, symmetry)?; + /// let csc = CscMatrix::new(nrow, ncol, col_pointers, row_indices, values, Sym::No)?; /// /// // covert to dense /// let a = csc.as_dense(); @@ -697,7 +695,7 @@ where if m != self.nrow || n != self.ncol { return Err("wrong matrix dimensions"); } - let mirror_required = self.symmetry.triangular(); + let mirror_required = self.symmetric.triangular(); a.fill(T::zero()); for j in 0..self.ncol { for p in self.col_pointers[j]..self.col_pointers[j + 1] { @@ -732,7 +730,7 @@ where if v.dim() != self.nrow { return Err("v vector is incompatible"); } - let mirror_required = self.symmetry.triangular(); + let mirror_required = self.symmetric.triangular(); v.fill(T::zero()); for j in 0..self.ncol { for p in self.col_pointers[j]..self.col_pointers[j + 1] { @@ -765,7 +763,7 @@ where /// let row_indices = vec![0, 0]; /// let values = vec![10.0, 20.0]; /// let csc = CscMatrix::new(1, 2, - /// col_pointers, row_indices, values, None)?; + /// col_pointers, row_indices, values, Sym::No)?; /// let (nrow, ncol, nnz, sym) = csc.get_info(); /// assert_eq!(nrow, 1); /// assert_eq!(ncol, 2); @@ -784,7 +782,7 @@ where self.nrow, self.ncol, self.col_pointers[self.ncol] as usize, - self.symmetry, + self.symmetric, ) } @@ -838,50 +836,50 @@ where #[cfg(test)] mod tests { use super::NumCscMatrix; - use crate::{CooMatrix, Samples, Storage, Sym}; + use crate::{CooMatrix, Samples, Sym}; use num_complex::Complex64; use russell_lab::{complex_vec_approx_eq, cpx, vec_approx_eq, ComplexVector, Matrix, Vector}; #[test] fn new_captures_errors() { assert_eq!( - NumCscMatrix::::new(0, 1, vec![0], vec![], vec![], None).err(), + NumCscMatrix::::new(0, 1, vec![0], vec![], vec![], Sym::No).err(), Some("nrow must be ≥ 1") ); assert_eq!( - NumCscMatrix::::new(1, 0, vec![0], vec![], vec![], None).err(), + NumCscMatrix::::new(1, 0, vec![0], vec![], vec![], Sym::No).err(), Some("ncol must be ≥ 1") ); assert_eq!( - NumCscMatrix::::new(1, 1, vec![0], vec![], vec![], None).err(), + NumCscMatrix::::new(1, 1, vec![0], vec![], vec![], Sym::No).err(), Some("col_pointers.len() must be = ncol + 1") ); assert_eq!( - NumCscMatrix::::new(1, 1, vec![0, 0], vec![], vec![], None).err(), + NumCscMatrix::::new(1, 1, vec![0, 0], vec![], vec![], Sym::No).err(), Some("nnz = col_pointers[ncol] must be ≥ 1") ); assert_eq!( - NumCscMatrix::::new(1, 1, vec![0, 1], vec![], vec![], None).err(), + NumCscMatrix::::new(1, 1, vec![0, 1], vec![], vec![], Sym::No).err(), Some("row_indices.len() must be ≥ nnz") ); assert_eq!( - NumCscMatrix::::new(1, 1, vec![0, 1], vec![0], vec![], None).err(), + NumCscMatrix::::new(1, 1, vec![0, 1], vec![0], vec![], Sym::No).err(), Some("values.len() must be ≥ nnz") ); assert_eq!( - NumCscMatrix::::new(1, 1, vec![-1, 1], vec![0], vec![0.0], None).err(), + NumCscMatrix::::new(1, 1, vec![-1, 1], vec![0], vec![0.0], Sym::No).err(), Some("col pointers must be ≥ 0") ); assert_eq!( - NumCscMatrix::::new(1, 1, vec![2, 1], vec![0], vec![0.0], None).err(), + NumCscMatrix::::new(1, 1, vec![2, 1], vec![0], vec![0.0], Sym::No).err(), Some("col pointers must be sorted in ascending order") ); assert_eq!( - NumCscMatrix::::new(1, 1, vec![0, 1], vec![-1], vec![0.0], None).err(), + NumCscMatrix::::new(1, 1, vec![0, 1], vec![-1], vec![0.0], Sym::No).err(), Some("row indices must be ≥ 0") ); assert_eq!( - NumCscMatrix::::new(1, 1, vec![0, 1], vec![2], vec![0.0], None).err(), + NumCscMatrix::::new(1, 1, vec![0, 1], vec![2], vec![0.0], Sym::No).err(), Some("row indices must be < nrow") ); // ┌ ┐ @@ -894,7 +892,7 @@ mod tests { let row_indices = vec![1, 0]; // << incorrect, should be [0, 1] let col_pointers = vec![0, 2]; assert_eq!( - NumCscMatrix::::new(2, 1, col_pointers, row_indices, values, None).err(), + NumCscMatrix::::new(2, 1, col_pointers, row_indices, values, Sym::No).err(), Some("row indices must be sorted in ascending order (within their column)") ); } @@ -902,8 +900,8 @@ mod tests { #[test] fn new_works() { let (_, csc_correct, _, _) = Samples::rectangular_1x2(false, false); - let csc = NumCscMatrix::::new(1, 2, vec![0, 1, 2], vec![0, 0], vec![10.0, 20.0], None).unwrap(); - assert_eq!(csc.symmetry, Sym::No); + let csc = NumCscMatrix::::new(1, 2, vec![0, 1, 2], vec![0, 0], vec![10.0, 20.0], Sym::No).unwrap(); + assert_eq!(csc.symmetric, Sym::No); assert_eq!(csc.nrow, 1); assert_eq!(csc.ncol, 2); assert_eq!(&csc.col_pointers, &csc_correct.col_pointers); @@ -913,7 +911,7 @@ mod tests { #[test] fn from_coo_captures_errors() { - let coo = CooMatrix::new(1, 1, 1, None).unwrap(); + let coo = CooMatrix::new(1, 1, 1, Sym::No).unwrap(); assert_eq!( NumCscMatrix::::from_coo(&coo).err(), Some("COO to CSC requires nnz > 0") @@ -1010,12 +1008,12 @@ mod tests { fn update_from_coo_captures_errors() { let (coo, _, _, _) = Samples::rectangular_1x2(false, false, ); let mut csc = NumCscMatrix::::from_coo(&coo).unwrap(); - let yes = Sym::General(Storage::Lower); + let yes = Sym::YesLower; let no = Sym::No; - assert_eq!(csc.update_from_coo(&CooMatrix { symmetry: yes, nrow: 1, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.symmetry must be equal to csc.symmetry")); - assert_eq!(csc.update_from_coo(&CooMatrix { symmetry: no, nrow: 2, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.nrow must be equal to csc.nrow")); - assert_eq!(csc.update_from_coo(&CooMatrix { symmetry: no, nrow: 1, ncol: 1, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.ncol must be equal to csc.ncol")); - assert_eq!(csc.update_from_coo(&CooMatrix { symmetry: no, nrow: 1, ncol: 2, nnz: 3, max_nnz: 3, indices_i: vec![0,0,0], indices_j: vec![0,0,0], values: vec![0.0,0.0,0.0] }).err(), Some("coo.nnz must be equal to nnz(dup) = csc.row_indices.len() = csc.values.len()")); + assert_eq!(csc.update_from_coo(&CooMatrix { symmetric: yes, nrow: 1, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.symmetry must be equal to csc.symmetry")); + assert_eq!(csc.update_from_coo(&CooMatrix { symmetric: no, nrow: 2, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.nrow must be equal to csc.nrow")); + assert_eq!(csc.update_from_coo(&CooMatrix { symmetric: no, nrow: 1, ncol: 1, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.ncol must be equal to csc.ncol")); + assert_eq!(csc.update_from_coo(&CooMatrix { symmetric: no, nrow: 1, ncol: 2, nnz: 3, max_nnz: 3, indices_i: vec![0,0,0], indices_j: vec![0,0,0], values: vec![0.0,0.0,0.0] }).err(), Some("coo.nnz must be equal to nnz(dup) = csc.row_indices.len() = csc.values.len()")); } #[test] @@ -1252,7 +1250,7 @@ mod tests { assert_eq!(csc.get_values(), &[10.0, 20.0]); // mutable let mut csc = NumCscMatrix:: { - symmetry: Sym::No, + symmetric: Sym::No, nrow: 1, ncol: 2, values: vec![10.0, 20.0], @@ -1288,10 +1286,10 @@ mod tests { let json = serde_json::to_string(&csc).unwrap(); assert_eq!( json, - r#"{"symmetry":"No","nrow":5,"ncol":5,"col_pointers":[0,2,5,9,10,12],"row_indices":[0,1,0,2,4,1,2,3,4,2,1,4,0],"values":[2.0,3.0,3.0,-1.0,4.0,4.0,-3.0,1.0,2.0,2.0,6.0,1.0,0.0]}"# + r#"{"symmetric":"No","nrow":5,"ncol":5,"col_pointers":[0,2,5,9,10,12],"row_indices":[0,1,0,2,4,1,2,3,4,2,1,4,0],"values":[2.0,3.0,3.0,-1.0,4.0,4.0,-3.0,1.0,2.0,2.0,6.0,1.0,0.0]}"# ); let from_json: NumCscMatrix = serde_json::from_str(&json).unwrap(); - assert_eq!(from_json.symmetry, csc.symmetry); + assert_eq!(from_json.symmetric, csc.symmetric); assert_eq!(from_json.nrow, csc.nrow); assert_eq!(from_json.ncol, csc.ncol); assert_eq!(from_json.col_pointers, csc.col_pointers); diff --git a/russell_sparse/src/csr_matrix.rs b/russell_sparse/src/csr_matrix.rs index 805b654f..8dbeb4d2 100644 --- a/russell_sparse/src/csr_matrix.rs +++ b/russell_sparse/src/csr_matrix.rs @@ -54,8 +54,8 @@ pub struct NumCsrMatrix where T: AddAssign + MulAssign + Num + NumCast + Copy + DeserializeOwned + Serialize, { - /// Defines the symmetry and storage: lower-triangular, upper-triangular, full-matrix - pub(crate) symmetry: Sym, + /// Indicates whether the matrix is symmetric or not. If symmetric, indicates the representation too. + pub(crate) symmetric: Sym, /// Holds the number of rows (must fit i32) pub(crate) nrow: usize, @@ -121,6 +121,8 @@ where /// to the number of non-zero values (sorted) /// * `col_indices` -- (len = nnz) column indices (sorted) /// * `values` -- the non-zero components of the matrix + /// * `symmetric` -- indicates whether the matrix is symmetric or not. + /// If symmetric, indicates the representation too. /// /// The following conditions must be satisfied (nnz is the number of non-zeros /// and nnz_dup is the number of non-zeros with possible duplicates): @@ -169,8 +171,7 @@ where /// 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 /// // count = 12 /// ]; - /// let symmetry = None; - /// let csr = CsrMatrix::new(nrow, ncol, row_pointers, col_indices, values, symmetry)?; + /// let csr = CsrMatrix::new(nrow, ncol, row_pointers, col_indices, values, Sym::No)?; /// /// // covert to dense /// let a = csr.as_dense(); @@ -191,7 +192,7 @@ where row_pointers: Vec, col_indices: Vec, values: Vec, - symmetry: Option, + symmetric: Sym, ) -> Result { if nrow < 1 { return Err("nrow must be ≥ 1"); @@ -238,7 +239,7 @@ where } } Ok(NumCsrMatrix { - symmetry: if let Some(v) = symmetry { v } else { Sym::No }, + symmetric, nrow, ncol, row_pointers, @@ -270,7 +271,7 @@ where /// // . . 1 . . /// // . 4 2 . 1 /// let (nrow, ncol, nnz) = (5, 5, 13); - /// let mut coo = CooMatrix::new(nrow, ncol, nnz, None)?; + /// let mut coo = CooMatrix::new(nrow, ncol, nnz, Sym::No)?; /// coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate /// coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate /// coo.put(1, 0, 3.0)?; @@ -324,7 +325,7 @@ where return Err("COO to CSR requires nnz > 0"); } let mut csr = NumCsrMatrix { - symmetry: coo.symmetry, + symmetric: coo.symmetric, nrow: coo.nrow, ncol: coo.ncol, row_pointers: vec![0; coo.nrow + 1], @@ -348,7 +349,7 @@ where /// may have been summed up. The final nnz is available as `nnz = row_pointers[nrow]`. pub fn update_from_coo(&mut self, coo: &NumCooMatrix) -> Result<(), StrError> { // check dimensions - if coo.symmetry != self.symmetry { + if coo.symmetric != self.symmetric { return Err("coo.symmetry must be equal to csr.symmetry"); } if coo.nrow != self.nrow { @@ -493,7 +494,7 @@ where // allocate the CSR arrays let mut csr = NumCsrMatrix { - symmetry: csc.symmetry, + symmetric: csc.symmetric, ncol: csc.ncol, nrow: csc.nrow, row_pointers: vec![0; nrow + 1], @@ -585,9 +586,7 @@ where /// 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 /// // count = 12 /// ]; - /// let symmetry = None; - /// let csr = CsrMatrix::new(nrow, ncol, - /// row_pointers, col_indices, values, symmetry)?; + /// let csr = CsrMatrix::new(nrow, ncol, row_pointers, col_indices, values, Sym::No)?; /// /// // covert to dense /// let a = csr.as_dense(); @@ -650,9 +649,7 @@ where /// 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 /// // count = 12 /// ]; - /// let symmetry = None; - /// let csr = CsrMatrix::new(nrow, ncol, - /// row_pointers, col_indices, values, symmetry)?; + /// let csr = CsrMatrix::new(nrow, ncol, row_pointers, col_indices, values, Sym::No)?; /// /// // covert to dense /// let a = csr.as_dense(); @@ -672,7 +669,7 @@ where if m != self.nrow || n != self.ncol { return Err("wrong matrix dimensions"); } - let mirror_required = self.symmetry.triangular(); + let mirror_required = self.symmetric.triangular(); a.fill(T::zero()); for i in 0..self.nrow { for p in self.row_pointers[i]..self.row_pointers[i + 1] { @@ -707,7 +704,7 @@ where if v.dim() != self.nrow { return Err("v vector is incompatible"); } - let mirror_required = self.symmetry.triangular(); + let mirror_required = self.symmetric.triangular(); v.fill(T::zero()); for i in 0..self.nrow { for p in self.row_pointers[i]..self.row_pointers[i + 1] { @@ -740,7 +737,7 @@ where /// let col_indices = vec![0, 1]; /// let values = vec![10.0, 20.0]; /// let csr = CsrMatrix::new(1, 2, - /// row_pointers, col_indices, values, None)?; + /// row_pointers, col_indices, values, Sym::No)?; /// let (nrow, ncol, nnz, sym) = csr.get_info(); /// assert_eq!(nrow, 1); /// assert_eq!(ncol, 2); @@ -759,7 +756,7 @@ where self.nrow, self.ncol, self.row_pointers[self.nrow] as usize, - self.symmetry, + self.symmetric, ) } @@ -813,50 +810,50 @@ where #[cfg(test)] mod tests { use super::NumCsrMatrix; - use crate::{CooMatrix, Samples, Storage, Sym}; + use crate::{CooMatrix, Samples, Sym}; use num_complex::Complex64; use russell_lab::{complex_vec_approx_eq, cpx, vec_approx_eq, ComplexVector, Matrix, Vector}; #[test] fn new_captures_errors() { assert_eq!( - NumCsrMatrix::::new(0, 1, vec![0], vec![], vec![], None).err(), + NumCsrMatrix::::new(0, 1, vec![0], vec![], vec![], Sym::No).err(), Some("nrow must be ≥ 1") ); assert_eq!( - NumCsrMatrix::::new(1, 0, vec![0], vec![], vec![], None).err(), + NumCsrMatrix::::new(1, 0, vec![0], vec![], vec![], Sym::No).err(), Some("ncol must be ≥ 1") ); assert_eq!( - NumCsrMatrix::::new(1, 1, vec![0], vec![], vec![], None).err(), + NumCsrMatrix::::new(1, 1, vec![0], vec![], vec![], Sym::No).err(), Some("row_pointers.len() must be = nrow + 1") ); assert_eq!( - NumCsrMatrix::::new(1, 1, vec![0, 0], vec![], vec![], None).err(), + NumCsrMatrix::::new(1, 1, vec![0, 0], vec![], vec![], Sym::No).err(), Some("nnz = row_pointers[nrow] must be ≥ 1") ); assert_eq!( - NumCsrMatrix::::new(1, 1, vec![0, 1], vec![], vec![], None).err(), + NumCsrMatrix::::new(1, 1, vec![0, 1], vec![], vec![], Sym::No).err(), Some("col_indices.len() must be ≥ nnz") ); assert_eq!( - NumCsrMatrix::::new(1, 1, vec![0, 1], vec![0], vec![], None).err(), + NumCsrMatrix::::new(1, 1, vec![0, 1], vec![0], vec![], Sym::No).err(), Some("values.len() must be ≥ nnz") ); assert_eq!( - NumCsrMatrix::::new(1, 1, vec![-1, 1], vec![0], vec![0.0], None).err(), + NumCsrMatrix::::new(1, 1, vec![-1, 1], vec![0], vec![0.0], Sym::No).err(), Some("row pointers must be ≥ 0") ); assert_eq!( - NumCsrMatrix::::new(1, 1, vec![2, 1], vec![0], vec![0.0], None).err(), + NumCsrMatrix::::new(1, 1, vec![2, 1], vec![0], vec![0.0], Sym::No).err(), Some("row pointers must be sorted in ascending order") ); assert_eq!( - NumCsrMatrix::::new(1, 1, vec![0, 1], vec![-1], vec![0.0], None).err(), + NumCsrMatrix::::new(1, 1, vec![0, 1], vec![-1], vec![0.0], Sym::No).err(), Some("column indices must be ≥ 0") ); assert_eq!( - NumCsrMatrix::::new(1, 1, vec![0, 1], vec![2], vec![0.0], None).err(), + NumCsrMatrix::::new(1, 1, vec![0, 1], vec![2], vec![0.0], Sym::No).err(), Some("column indices must be < ncol") ); // ┌ ┐ @@ -868,7 +865,7 @@ mod tests { let col_indices = vec![1, 0]; // << incorrect, should be [0, 1] let row_pointers = vec![0, 2]; assert_eq!( - NumCsrMatrix::::new(1, 2, row_pointers, col_indices, values, None).err(), + NumCsrMatrix::::new(1, 2, row_pointers, col_indices, values, Sym::No).err(), Some("column indices must be sorted in ascending order (within their row)") ); } @@ -876,8 +873,8 @@ mod tests { #[test] fn new_works() { let (_, _, csr_correct, _) = Samples::rectangular_1x2(false, false); - let csr = NumCsrMatrix::::new(1, 2, vec![0, 2], vec![0, 1], vec![10.0, 20.0], None).unwrap(); - assert_eq!(csr.symmetry, Sym::No); + let csr = NumCsrMatrix::::new(1, 2, vec![0, 2], vec![0, 1], vec![10.0, 20.0], Sym::No).unwrap(); + assert_eq!(csr.symmetric, Sym::No); assert_eq!(csr.nrow, 1); assert_eq!(csr.ncol, 2); assert_eq!(&csr.row_pointers, &csr_correct.row_pointers); @@ -887,7 +884,7 @@ mod tests { #[test] fn from_coo_captures_errors() { - let coo = CooMatrix::new(1, 1, 1, None).unwrap(); + let coo = CooMatrix::new(1, 1, 1, Sym::No).unwrap(); assert_eq!( NumCsrMatrix::::from_coo(&coo).err(), Some("COO to CSR requires nnz > 0") @@ -999,12 +996,12 @@ mod tests { fn update_from_coo_captures_errors() { let (coo, _, _, _) = Samples::rectangular_1x2(false, false); let mut csr = NumCsrMatrix::::from_coo(&coo).unwrap(); - let yes = Sym::General(Storage::Lower); + let yes = Sym::YesLower; let no = Sym::No; - assert_eq!(csr.update_from_coo(&CooMatrix { symmetry: yes, nrow: 1, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.symmetry must be equal to csr.symmetry")); - assert_eq!(csr.update_from_coo(&CooMatrix { symmetry: no, nrow: 2, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.nrow must be equal to csr.nrow")); - assert_eq!(csr.update_from_coo(&CooMatrix { symmetry: no, nrow: 1, ncol: 1, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.ncol must be equal to csr.ncol")); - assert_eq!(csr.update_from_coo(&CooMatrix { symmetry: no, nrow: 1, ncol: 2, nnz: 3, max_nnz: 3, indices_i: vec![0,0,0], indices_j: vec![0,0,0], values: vec![0.0,0.0,0.0] }).err(), Some("coo.nnz must be equal to nnz(dup) = self.col_indices.len() = csr.values.len()")); + assert_eq!(csr.update_from_coo(&CooMatrix { symmetric: yes, nrow: 1, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.symmetry must be equal to csr.symmetry")); + assert_eq!(csr.update_from_coo(&CooMatrix { symmetric: no, nrow: 2, ncol: 2, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.nrow must be equal to csr.nrow")); + assert_eq!(csr.update_from_coo(&CooMatrix { symmetric: no, nrow: 1, ncol: 1, nnz: 1, max_nnz: 1, indices_i: vec![0], indices_j: vec![0], values: vec![0.0] }).err(), Some("coo.ncol must be equal to csr.ncol")); + assert_eq!(csr.update_from_coo(&CooMatrix { symmetric: no, nrow: 1, ncol: 2, nnz: 3, max_nnz: 3, indices_i: vec![0,0,0], indices_j: vec![0,0,0], values: vec![0.0,0.0,0.0] }).err(), Some("coo.nnz must be equal to nnz(dup) = self.col_indices.len() = csr.values.len()")); } #[test] @@ -1241,7 +1238,7 @@ mod tests { assert_eq!(csr.get_values(), &[10.0, 20.0]); // mutable let mut csr = NumCsrMatrix:: { - symmetry: Sym::No, + symmetric: Sym::No, nrow: 1, ncol: 2, values: vec![10.0, 20.0], @@ -1275,10 +1272,10 @@ mod tests { let json = serde_json::to_string(&csr).unwrap(); assert_eq!( json, - r#"{"symmetry":"No","nrow":5,"ncol":5,"row_pointers":[0,2,5,8,9,12],"col_indices":[0,1,0,2,4,1,2,3,2,1,2,4,0],"values":[2.0,3.0,3.0,4.0,6.0,-1.0,-3.0,2.0,1.0,4.0,2.0,1.0,0.0]}"# + r#"{"symmetric":"No","nrow":5,"ncol":5,"row_pointers":[0,2,5,8,9,12],"col_indices":[0,1,0,2,4,1,2,3,2,1,2,4,0],"values":[2.0,3.0,3.0,4.0,6.0,-1.0,-3.0,2.0,1.0,4.0,2.0,1.0,0.0]}"# ); let from_json: NumCsrMatrix = serde_json::from_str(&json).unwrap(); - assert_eq!(from_json.symmetry, csr.symmetry); + assert_eq!(from_json.symmetric, csr.symmetric); assert_eq!(from_json.nrow, csr.nrow); assert_eq!(from_json.ncol, csr.ncol); assert_eq!(from_json.row_pointers, csr.row_pointers); diff --git a/russell_sparse/src/enums.rs b/russell_sparse/src/enums.rs index c8d702fd..00246d68 100644 --- a/russell_sparse/src/enums.rs +++ b/russell_sparse/src/enums.rs @@ -1,4 +1,3 @@ -use crate::StrError; use serde::{Deserialize, Serialize}; /// Specifies the underlying library that does all the magic @@ -15,37 +14,27 @@ pub enum Genie { Umfpack, } -/// Specifies how the matrix components are stored -#[derive(Clone, Copy, Debug, Eq, PartialEq, Deserialize, Serialize)] -pub enum Storage { - /// Lower triangular storage for symmetric matrix (e.g., for MUMPS) - Lower, - - /// Upper triangular storage for symmetric matrix - Upper, - - /// Full matrix storage for symmetric or unsymmetric matrix (e.g., for UMFPACK) - Full, -} - /// Specifies the type of matrix symmetry #[derive(Clone, Copy, Debug, Eq, PartialEq, Deserialize, Serialize)] pub enum Sym { /// Unknown symmetry (possibly unsymmetric) No, - /// General symmetric - General(Storage), + /// Symmetric with full representation (i.e., not triangular) + YesFull, + + /// Symmetric with lower-triangle representation + YesLower, - /// Symmetric and positive-definite - PositiveDefinite(Storage), + /// Symmetric with upper-triangle representation + YesUpper, } /// Holds options to handle a MatrixMarket when the matrix is specified as being symmetric /// -/// **Note:** This is ignored if not the matrix is not specified as symmetric. +/// **Note:** This is ignored if the matrix is not symmetric. #[derive(Clone, Copy, Debug, Eq, PartialEq, Deserialize, Serialize)] -pub enum MMsymOption { +pub enum MMsym { /// Leave the storage as lower triangular (if symmetric) /// /// **Note:** Lower triangular is the standard MatrixMarket format. @@ -162,27 +151,12 @@ impl Genie { } } - /// Returns which storage is required by the solver - /// - /// ```text - /// MUMPS : Storage::Lower - /// UMFPACK : Storage::Full - /// ```` - pub fn storage(&self) -> Storage { - match self { - Genie::Mumps => Storage::Lower, - Genie::Umfpack => Storage::Full, - } - } - - /// Returns the solver's required symmetry/storage configuration - pub fn symmetry(&self, symmetric: bool, positive_definite: bool) -> Sym { - let storage = self.storage(); - if symmetric || positive_definite { - if positive_definite { - Sym::PositiveDefinite(storage) - } else { - Sym::General(storage) + /// Returns the solver's required symmetry-representation + pub fn symmetry(&self, symmetric: bool) -> Sym { + if symmetric { + match self { + Genie::Mumps => Sym::YesLower, + Genie::Umfpack => Sym::YesFull, } } else { Sym::No @@ -191,108 +165,12 @@ impl Genie { } impl Sym { - /// Returns a new general symmetry flag with lower storage - pub fn new_general_lower() -> Self { - Sym::General(Storage::Lower) - } - - /// Returns a new general symmetry flag with upper storage - pub fn new_general_upper() -> Self { - Sym::General(Storage::Upper) - } - - /// Returns a new general symmetry flag with full storage - pub fn new_general_full() -> Self { - Sym::General(Storage::Full) - } - - /// Returns a new positive-definite symmetry flag with lower storage - pub fn new_pos_def_lower() -> Self { - Sym::PositiveDefinite(Storage::Lower) - } - - /// Returns a new positive-definite symmetry flag with upper storage - pub fn new_pos_def_upper() -> Self { - Sym::PositiveDefinite(Storage::Upper) - } - - /// Returns a new positive-definite symmetry flag with full storage - pub fn new_pos_def_full() -> Self { - Sym::PositiveDefinite(Storage::Full) - } - - /// Returns which type of storage is used, if symmetric - pub fn storage(symmetry: Sym) -> Storage { - if symmetry.lower() { - Storage::Lower - } else if symmetry.upper() { - Storage::Upper - } else { - Storage::Full - } - } - - /// Returns true if the storage is triangular (lower or upper) + /// Returns true if the representation is Lower or Upper pub fn triangular(&self) -> bool { match self { - Sym::No => false, - Sym::General(storage) => *storage != Storage::Full, - Sym::PositiveDefinite(storage) => *storage != Storage::Full, - } - } - - /// Returns true if the storage is lower triangular - pub fn lower(&self) -> bool { - match self { - Sym::No => false, - Sym::General(storage) => *storage == Storage::Lower, - Sym::PositiveDefinite(storage) => *storage == Storage::Lower, - } - } - - /// Returns true if the storage is upper triangular - pub fn upper(&self) -> bool { - match self { - Sym::No => false, - Sym::General(storage) => *storage == Storage::Upper, - Sym::PositiveDefinite(storage) => *storage == Storage::Upper, - } - } - - /// Returns status flags indicating the type of symmetry, if any - /// - /// # Input - /// - /// * `must_be_lower` -- makes sure that the storage is Lower - /// * `must_be_upper` -- makes sure that the storage is Upper - /// - /// # Output - /// - /// Returns `(general_symmetric, positive_definite)` where: - /// - /// * `general_symmetric` -- 1 if true, 0 otherwise - /// * `positive_definite` -- 1 if true, 0 otherwise - pub fn status(&self, must_be_lower: bool, must_be_upper: bool) -> Result<(i32, i32), StrError> { - match self { - Sym::No => Ok((0, 0)), - Sym::General(storage) => { - if must_be_lower && *storage != Storage::Lower { - return Err("if the matrix is general symmetric, the required storage is lower triangular"); - } - if must_be_upper && *storage != Storage::Upper { - return Err("if the matrix is general symmetric, the required storage is upper triangular"); - } - Ok((1, 0)) - } - Sym::PositiveDefinite(storage) => { - if must_be_lower && *storage != Storage::Lower { - return Err("if the matrix is positive-definite, the required storage is lower triangular"); - } - if must_be_upper && *storage != Storage::Upper { - return Err("if the matrix is positive-definite, the required storage is upper triangular"); - } - Ok((0, 1)) - } + Sym::YesLower => true, + Sym::YesUpper => true, + _ => false, } } } @@ -377,34 +255,24 @@ mod tests { let from_json: Genie = serde_json::from_str(&json).unwrap(); assert_eq!(from_json, genie); - let storage = Storage::Full; - let copy = storage; - let clone = storage.clone(); - assert_eq!(format!("{:?}", storage), "Full"); - assert_eq!(copy, Storage::Full); - assert_eq!(clone, Storage::Full); - let json = serde_json::to_string(&storage).unwrap(); - let from_json: Storage = serde_json::from_str(&json).unwrap(); - assert_eq!(from_json, storage); - - let symmetry = Sym::PositiveDefinite(Storage::Lower); + let symmetry = Sym::YesLower; let copy = symmetry; let clone = symmetry.clone(); - assert_eq!(format!("{:?}", symmetry), "PositiveDefinite(Lower)"); - assert_eq!(copy, Sym::PositiveDefinite(Storage::Lower)); - assert_eq!(clone, Sym::PositiveDefinite(Storage::Lower)); + assert_eq!(format!("{:?}", symmetry), "YesLower"); + assert_eq!(copy, Sym::YesLower); + assert_eq!(clone, Sym::YesLower); let json = serde_json::to_string(&symmetry).unwrap(); let from_json: Sym = serde_json::from_str(&json).unwrap(); assert_eq!(from_json, symmetry); - let handling = MMsymOption::LeaveAsLower; + let handling = MMsym::LeaveAsLower; let copy = handling; let clone = handling.clone(); assert_eq!(format!("{:?}", handling), "LeaveAsLower"); - assert_eq!(copy, MMsymOption::LeaveAsLower); - assert_eq!(clone, MMsymOption::LeaveAsLower); + assert_eq!(copy, MMsym::LeaveAsLower); + assert_eq!(clone, MMsym::LeaveAsLower); let json = serde_json::to_string(&handling).unwrap(); - let from_json: MMsymOption = serde_json::from_str(&json).unwrap(); + let from_json: MMsym = serde_json::from_str(&json).unwrap(); assert_eq!(from_json, handling); let ordering = Ordering::Amd; @@ -489,106 +357,22 @@ mod tests { assert_eq!(Genie::from("Mumps"), Genie::Mumps); assert_eq!(Genie::from("Umfpack"), Genie::Umfpack); - let l = Storage::Lower; - let f = Storage::Full; - - let gl = Sym::General(l); - let gf = Sym::General(f); - - let pl = Sym::PositiveDefinite(l); - let pf = Sym::PositiveDefinite(f); - let genie = Genie::Mumps; assert_eq!(genie.to_string(), "mumps"); - assert_eq!(genie.storage(), l); - assert_eq!(genie.symmetry(false, false), Sym::No); - assert_eq!(genie.symmetry(true, false), gl); - assert_eq!(genie.symmetry(false, true), pl); - assert_eq!(genie.symmetry(true, true), pl); + assert_eq!(genie.symmetry(false), Sym::No); + assert_eq!(genie.symmetry(true), Sym::YesLower); let genie = Genie::Umfpack; assert_eq!(genie.to_string(), "umfpack"); - assert_eq!(genie.storage(), f); - assert_eq!(genie.symmetry(false, false), Sym::No); - assert_eq!(genie.symmetry(true, false), gf); - assert_eq!(genie.symmetry(false, true), pf); - assert_eq!(genie.symmetry(true, true), pf); + assert_eq!(genie.symmetry(false,), Sym::No); + assert_eq!(genie.symmetry(true), Sym::YesFull); } #[test] fn symmetry_functions_work() { - let l = Storage::Lower; - let u = Storage::Upper; - let f = Storage::Full; - - let gl = Sym::General(l); - let gu = Sym::General(u); - let gf = Sym::General(f); - - let pl = Sym::PositiveDefinite(l); - let pu = Sym::PositiveDefinite(u); - let pf = Sym::PositiveDefinite(f); - - assert_eq!(Sym::new_general_lower(), gl); - assert_eq!(Sym::new_general_upper(), gu); - assert_eq!(Sym::new_general_full(), gf); - assert_eq!(Sym::storage(gl), Storage::Lower); - assert_eq!(Sym::storage(gu), Storage::Upper); - assert_eq!(Sym::storage(gf), Storage::Full); - assert_eq!(gl.triangular(), true); - assert_eq!(gu.triangular(), true); - assert_eq!(gf.triangular(), false); - assert_eq!(gl.lower(), true); - assert_eq!(gu.lower(), false); - assert_eq!(gf.lower(), false); - assert_eq!(gl.upper(), false); - assert_eq!(gu.upper(), true); - assert_eq!(gf.upper(), false); - - assert_eq!(gl.status(true, false), Ok((1, 0))); - assert_eq!(gl.status(false, false), Ok((1, 0))); - assert_eq!( - gl.status(false, true), - Err("if the matrix is general symmetric, the required storage is upper triangular") - ); - - assert_eq!(gu.status(false, true), Ok((1, 0))); - assert_eq!(gu.status(false, false), Ok((1, 0))); - assert_eq!( - gu.status(true, false), - Err("if the matrix is general symmetric, the required storage is lower triangular") - ); - - assert_eq!(Sym::new_pos_def_lower(), pl); - assert_eq!(Sym::new_pos_def_upper(), pu); - assert_eq!(Sym::new_pos_def_full(), pf); - assert_eq!(Sym::storage(pl), Storage::Lower); - assert_eq!(Sym::storage(pu), Storage::Upper); - assert_eq!(Sym::storage(pf), Storage::Full); - assert_eq!(pl.triangular(), true); - assert_eq!(pu.triangular(), true); - assert_eq!(pf.triangular(), false); - assert_eq!(pl.lower(), true); - assert_eq!(pu.lower(), false); - assert_eq!(pf.lower(), false); - assert_eq!(pl.upper(), false); - assert_eq!(pu.upper(), true); - assert_eq!(pf.upper(), false); - - assert_eq!(pl.status(true, false), Ok((0, 1))); - assert_eq!(pl.status(false, false), Ok((0, 1))); - assert_eq!( - pl.status(false, true), - Err("if the matrix is positive-definite, the required storage is upper triangular") - ); - - assert_eq!(pu.status(false, true), Ok((0, 1))); - assert_eq!(pu.status(false, false), Ok((0, 1))); - assert_eq!( - pu.status(true, false), - Err("if the matrix is positive-definite, the required storage is lower triangular") - ); - - assert_eq!(Sym::storage(Sym::No), Storage::Full); + assert_eq!(Sym::No.triangular(), false); + assert_eq!(Sym::YesFull.triangular(), false); + assert_eq!(Sym::YesLower.triangular(), true); + assert_eq!(Sym::YesUpper.triangular(), true); } } diff --git a/russell_sparse/src/lib.rs b/russell_sparse/src/lib.rs index d62765f9..ea401f66 100644 --- a/russell_sparse/src/lib.rs +++ b/russell_sparse/src/lib.rs @@ -80,7 +80,7 @@ //! // │ 4 5 6 │ //! // └ ┘ //! let (nrow, ncol, nnz) = (3, 3, 6); -//! let mut coo = CooMatrix::new(nrow, ncol, nnz, None)?; +//! let mut coo = CooMatrix::new(nrow, ncol, nnz, Sym::No)?; //! coo.put(0, 0, 1.0)?; //! coo.put(0, 2, 2.0)?; //! coo.put(1, 2, 3.0)?; @@ -130,7 +130,7 @@ //! let mut solver = LinSolver::new(Genie::Umfpack)?; //! //! // allocate the coefficient matrix -//! let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?; +//! let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; //! coo.put(0, 0, 0.2)?; //! coo.put(0, 1, 0.2)?; //! coo.put(1, 0, 0.5)?; @@ -189,7 +189,7 @@ //! // . -1 -3 2 . //! // . . 1 . . //! // . 4 2 . 1 -//! let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?; +//! let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; //! coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate //! coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate //! coo.put(1, 0, 3.0)?; diff --git a/russell_sparse/src/lin_sol_params.rs b/russell_sparse/src/lin_sol_params.rs index 2ae61b84..0601f86c 100644 --- a/russell_sparse/src/lin_sol_params.rs +++ b/russell_sparse/src/lin_sol_params.rs @@ -9,6 +9,9 @@ pub struct LinSolParams { /// Defines the scaling strategy pub scaling: Scaling, + /// Indicates that the coefficient matrix is positive-definite (only considered if the matrix is symmetric) + pub positive_definite: bool, + /// Requests that the determinant be computed /// /// **Note:** The determinant will be available after `factorize` @@ -55,6 +58,7 @@ impl LinSolParams { LinSolParams { ordering: Ordering::Auto, scaling: Scaling::Auto, + positive_definite: false, compute_determinant: false, compute_error_estimates: false, compute_condition_numbers: false, diff --git a/russell_sparse/src/lin_solver.rs b/russell_sparse/src/lin_solver.rs index 60a66fac..4ed8f47c 100644 --- a/russell_sparse/src/lin_solver.rs +++ b/russell_sparse/src/lin_solver.rs @@ -109,7 +109,7 @@ impl<'a> LinSolver<'a> { /// let nnz = 5; // number of non-zero values /// /// // allocate the coefficient matrix - /// let mut mat = SparseMatrix::new_coo(ndim, ndim, nnz, None)?; + /// let mut mat = SparseMatrix::new_coo(ndim, ndim, nnz, Sym::No)?; /// mat.put(0, 0, 0.2)?; /// mat.put(0, 1, 0.2)?; /// mat.put(1, 0, 0.5)?; diff --git a/russell_sparse/src/read_matrix_market.rs b/russell_sparse/src/read_matrix_market.rs index 4659f33a..d2d128b8 100644 --- a/russell_sparse/src/read_matrix_market.rs +++ b/russell_sparse/src/read_matrix_market.rs @@ -1,4 +1,4 @@ -use super::{CooMatrix, MMsymOption, Storage, Sym}; +use super::{CooMatrix, MMsym, Sym}; use crate::StrError; use std::ffi::OsStr; use std::fs::File; @@ -264,7 +264,7 @@ impl MatrixMarketData { /// /// fn main() -> Result<(), StrError> { /// let name = "./data/matrix_market/ok_simple_general.mtx"; -/// let coo = read_matrix_market(name, MMsymOption::LeaveAsLower)?; +/// let coo = read_matrix_market(name, MMsym::LeaveAsLower)?; /// let (nrow, ncol, nnz, sym) = coo.get_info(); /// assert_eq!(nrow, 3); /// assert_eq!(ncol, 3); @@ -303,12 +303,12 @@ impl MatrixMarketData { /// /// fn main() -> Result<(), StrError> { /// let name = "./data/matrix_market/ok_simple_symmetric.mtx"; -/// let coo = read_matrix_market(name, MMsymOption::LeaveAsLower)?; +/// let coo = read_matrix_market(name, MMsym::LeaveAsLower)?; /// let (nrow, ncol, nnz, sym) = coo.get_info(); /// assert_eq!(nrow, 3); /// assert_eq!(ncol, 3); /// assert_eq!(nnz, 4); -/// assert_eq!(sym, Sym::General(Storage::Lower)); +/// assert_eq!(sym, Sym::YesLower); /// let a = coo.as_dense(); /// let correct = "┌ ┐\n\ /// │ 1 2 0 │\n\ @@ -319,7 +319,7 @@ impl MatrixMarketData { /// Ok(()) /// } /// ``` -pub fn read_matrix_market

(full_path: &P, symmetric_handling: MMsymOption) -> Result +pub fn read_matrix_market

(full_path: &P, symmetric_handling: MMsym) -> Result where P: AsRef + ?Sized, { @@ -349,27 +349,27 @@ where } // symmetry option - let symmetry = if data.symmetric { + let sym = if data.symmetric { if data.m != data.n { - return Err("MatrixMarket data is invalid: the number of rows must be equal the number of columns for symmetric matrices"); + return Err("MatrixMarket data is invalid: the number of rows must equal the number of columns for symmetric matrices"); } match symmetric_handling { - MMsymOption::LeaveAsLower => Some(Sym::General(Storage::Lower)), - MMsymOption::SwapToUpper => Some(Sym::General(Storage::Upper)), - MMsymOption::MakeItFull => Some(Sym::General(Storage::Full)), + MMsym::LeaveAsLower => Sym::YesLower, + MMsym::SwapToUpper => Sym::YesUpper, + MMsym::MakeItFull => Sym::YesFull, } } else { - None + Sym::No }; // set max number of entries let mut max = data.nnz; - if data.symmetric && symmetric_handling == MMsymOption::MakeItFull { + if data.symmetric && symmetric_handling == MMsym::MakeItFull { max = 2 * data.nnz; } // allocate triplet - let mut coo = CooMatrix::new(data.m as usize, data.n as usize, max as usize, symmetry).unwrap(); + let mut coo = CooMatrix::new(data.m as usize, data.n as usize, max as usize, sym).unwrap(); // read and parse triples loop { @@ -379,13 +379,13 @@ where if data.parse_triple(&line)? { if data.symmetric { match symmetric_handling { - MMsymOption::LeaveAsLower => { + MMsym::LeaveAsLower => { coo.put(data.i as usize, data.j as usize, data.aij).unwrap(); } - MMsymOption::SwapToUpper => { + MMsym::SwapToUpper => { coo.put(data.j as usize, data.i as usize, data.aij).unwrap(); } - MMsymOption::MakeItFull => { + MMsym::MakeItFull => { coo.put(data.i as usize, data.j as usize, data.aij).unwrap(); if data.i != data.j { coo.put(data.j as usize, data.i as usize, data.aij).unwrap(); @@ -414,7 +414,7 @@ where #[cfg(test)] mod tests { use super::{read_matrix_market, MatrixMarketData}; - use crate::{MMsymOption, Storage, Sym}; + use crate::{MMsym, Sym}; use russell_lab::Matrix; #[test] @@ -531,7 +531,7 @@ mod tests { #[test] fn read_matrix_market_handle_wrong_files() { - let h = MMsymOption::LeaveAsLower; + let h = MMsym::LeaveAsLower; assert_eq!(read_matrix_market("__wrong__", h).err(), Some("cannot open file")); assert_eq!( read_matrix_market("./data/matrix_market/bad_empty_file.mtx", h).err(), @@ -555,16 +555,16 @@ mod tests { ); assert_eq!( read_matrix_market("./data/matrix_market/bad_symmetric_rectangular.mtx", h).err(), - Some("MatrixMarket data is invalid: the number of rows must be equal the number of columns for symmetric matrices") + Some("MatrixMarket data is invalid: the number of rows must equal the number of columns for symmetric matrices") ); } #[test] fn read_matrix_market_works() { - let h = MMsymOption::LeaveAsLower; + let h = MMsym::LeaveAsLower; let filepath = "./data/matrix_market/ok_general.mtx".to_string(); let coo = read_matrix_market(&filepath, h).unwrap(); - assert_eq!(coo.symmetry, Sym::No); + assert_eq!(coo.symmetric, Sym::No); assert_eq!((coo.nrow, coo.ncol, coo.nnz, coo.max_nnz), (5, 5, 12, 12)); assert_eq!(coo.indices_i, &[0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4]); assert_eq!(coo.indices_j, &[0, 0, 1, 1, 1, 2, 2, 2, 2, 3, 4, 4]); @@ -576,10 +576,10 @@ mod tests { #[test] fn read_matrix_market_symmetric_lower_works() { - let h = MMsymOption::LeaveAsLower; + let h = MMsym::LeaveAsLower; let filepath = "./data/matrix_market/ok_symmetric.mtx".to_string(); let coo = read_matrix_market(&filepath, h).unwrap(); - assert_eq!(coo.symmetry, Sym::General(Storage::Lower)); + assert_eq!(coo.symmetric, Sym::YesLower); assert_eq!((coo.nrow, coo.ncol, coo.nnz, coo.max_nnz), (5, 5, 15, 15)); assert_eq!(coo.indices_i, &[0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4]); assert_eq!(coo.indices_j, &[0, 1, 2, 3, 4, 0, 0, 0, 0, 1, 1, 1, 2, 2, 3]); @@ -591,10 +591,10 @@ mod tests { #[test] fn read_matrix_market_symmetric_upper_works() { - let h = MMsymOption::SwapToUpper; + let h = MMsym::SwapToUpper; let filepath = "./data/matrix_market/ok_symmetric.mtx".to_string(); let coo = read_matrix_market(&filepath, h).unwrap(); - assert_eq!(coo.symmetry, Sym::General(Storage::Upper)); + assert_eq!(coo.symmetric, Sym::YesUpper); assert_eq!((coo.nrow, coo.ncol, coo.nnz, coo.max_nnz), (5, 5, 15, 15)); assert_eq!(coo.indices_i, &[0, 1, 2, 3, 4, 0, 0, 0, 0, 1, 1, 1, 2, 2, 3]); assert_eq!(coo.indices_j, &[0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4]); @@ -606,10 +606,10 @@ mod tests { #[test] fn read_matrix_market_symmetric_to_full_works() { - let h = MMsymOption::MakeItFull; + let h = MMsym::MakeItFull; let filepath = "./data/matrix_market/ok_symmetric_small.mtx".to_string(); let coo = read_matrix_market(&filepath, h).unwrap(); - assert_eq!(coo.symmetry, Sym::General(Storage::Full)); + assert_eq!(coo.symmetric, Sym::YesFull); assert_eq!((coo.nrow, coo.ncol, coo.nnz, coo.max_nnz), (5, 5, 11, 14)); assert_eq!(coo.indices_i, &[0, 1, 0, 2, 1, 3, 2, 3, 4, 1, 4, 0, 0, 0]); assert_eq!(coo.indices_j, &[0, 0, 1, 1, 2, 2, 3, 3, 1, 4, 4, 0, 0, 0]); diff --git a/russell_sparse/src/samples.rs b/russell_sparse/src/samples.rs index ad23f452..f3e9bf48 100644 --- a/russell_sparse/src/samples.rs +++ b/russell_sparse/src/samples.rs @@ -1,5 +1,5 @@ use crate::{ComplexCooMatrix, ComplexCscMatrix, ComplexCsrMatrix}; -use crate::{CooMatrix, CscMatrix, CsrMatrix, Storage, Sym}; +use crate::{CooMatrix, CscMatrix, CsrMatrix, Sym}; use num_complex::Complex64; use russell_lab::cpx; @@ -17,7 +17,7 @@ impl Samples { /// └ ┘ /// ``` pub fn tiny_1x1() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = None; + let sym = Sym::No; let nrow = 1; let ncol = 1; let max_nnz = 1; @@ -44,7 +44,7 @@ impl Samples { /// └ ┘ /// ``` pub fn complex_tiny_1x1() -> (ComplexCooMatrix, ComplexCscMatrix, ComplexCsrMatrix, Complex64) { - let sym = None; + let sym = Sym::No; let nrow = 1; let ncol = 1; let max_nnz = 1; @@ -70,9 +70,9 @@ impl Samples { /// -1 2 -1 => -1 2 /// -1 2 -1 2 /// ``` - pub fn positive_definite_3x3() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { + pub fn positive_definite_3x3_lower() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { let (nrow, ncol, nnz) = (3, 3, 6); - let sym = Some(Sym::PositiveDefinite(Storage::Lower)); + let sym = Sym::YesLower; let mut coo = CooMatrix::new(nrow, ncol, nnz, sym).unwrap(); coo.put(1, 0, -0.5).unwrap(); // duplicate coo.put(0, 0, 2.0).unwrap(); @@ -118,7 +118,7 @@ impl Samples { /// ``` pub fn complex_symmetric_3x3_lower() -> (ComplexCooMatrix, ComplexCscMatrix, ComplexCsrMatrix, Complex64) { let (nrow, ncol, nnz) = (3, 3, 6); - let sym = Some(Sym::General(Storage::Lower)); + let sym = Sym::YesLower; let mut coo = ComplexCooMatrix::new(nrow, ncol, nnz, sym).unwrap(); coo.put(1, 0, cpx!(-0.5, -0.5)).unwrap(); // duplicate coo.put(0, 0, cpx!(2.0, 1.0)).unwrap(); @@ -166,7 +166,7 @@ impl Samples { /// ``` pub fn complex_symmetric_3x3_full() -> (ComplexCooMatrix, ComplexCscMatrix, ComplexCsrMatrix, Complex64) { let (nrow, ncol, nnz) = (3, 3, 8); - let sym = Some(Sym::General(Storage::Full)); + let sym = Sym::YesFull; let mut coo = ComplexCooMatrix::new(nrow, ncol, nnz, sym).unwrap(); coo.put(1, 0, cpx!(-0.5, -0.5)).unwrap(); // duplicate coo.put(0, 0, cpx!(2.0, 1.0)).unwrap(); @@ -218,7 +218,7 @@ impl Samples { /// ``` pub fn lower_symmetric_5x5() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { let (nrow, ncol, nnz) = (5, 5, 18); - let sym = Some(Sym::PositiveDefinite(Storage::Lower)); + let sym = Sym::YesLower; let mut coo = CooMatrix::new(nrow, ncol, nnz, sym).unwrap(); coo.put(1, 1, 2.0).unwrap(); coo.put(4, 2, 2.5).unwrap(); // duplicate @@ -302,7 +302,7 @@ impl Samples { shuffle_coo_entries: bool, duplicate_coo_entries: bool, ) -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = None; + let sym = Sym::No; let nrow = 3; let ncol = 3; let max_nnz = 10; // more nnz than needed => OK @@ -401,7 +401,7 @@ impl Samples { /// let x_correct = &[1.0, 2.0, 3.0, 4.0, 5.0]; /// ``` pub fn umfpack_unsymmetric_5x5() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = None; + let sym = Sym::No; let nrow = 5; let ncol = 5; let max_nnz = 13; @@ -471,7 +471,7 @@ impl Samples { /// Reference: /// pub fn mkl_unsymmetric_5x5() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = None; + let sym = Sym::No; let nrow = 5; let ncol = 5; let mut coo = CooMatrix::new(nrow, ncol, 13, sym).unwrap(); @@ -540,7 +540,7 @@ impl Samples { shuffle_coo_entries: bool, duplicate_coo_entries: bool, ) -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = None; + let sym = Sym::No; let nrow = 5; let ncol = 5; let max_nnz = 11; // more nnz than needed => OK @@ -655,7 +655,7 @@ impl Samples { /// x_correct = vec![-979.0 / 3.0, 983.0, 1961.0 / 12.0, 398.0, 123.0 / 2.0]; /// ``` pub fn mkl_positive_definite_5x5_lower() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = Some(Sym::PositiveDefinite(Storage::Lower)); + let sym = Sym::YesLower; let nrow = 5; let ncol = 5; let mut coo = CooMatrix::new(nrow, ncol, 9, sym).unwrap(); @@ -729,7 +729,7 @@ impl Samples { /// x_correct = vec![-979.0 / 3.0, 983.0, 1961.0 / 12.0, 398.0, 123.0 / 2.0]; /// ``` pub fn mkl_positive_definite_5x5_upper() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = Some(Sym::PositiveDefinite(Storage::Upper)); + let sym = Sym::YesUpper; let nrow = 5; let ncol = 5; let mut coo = CooMatrix::new(nrow, ncol, 9, sym).unwrap(); @@ -806,7 +806,7 @@ impl Samples { shuffle_coo_entries: bool, duplicate_coo_entries: bool, ) -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = Some(Sym::General(Storage::Lower)); + let sym = Sym::YesLower; let nrow = 5; let ncol = 5; let max_nnz = 13; @@ -928,7 +928,7 @@ impl Samples { shuffle_coo_entries: bool, duplicate_coo_entries: bool, ) -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = Some(Sym::General(Storage::Upper)); + let sym = Sym::YesUpper; let nrow = 5; let ncol = 5; let max_nnz = 15; @@ -1046,7 +1046,7 @@ impl Samples { /// x_correct = vec![-979.0 / 3.0, 983.0, 1961.0 / 12.0, 398.0, 123.0 / 2.0]; /// ``` pub fn mkl_symmetric_5x5_full() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = Some(Sym::General(Storage::Full)); + let sym = Sym::YesFull; let nrow = 5; let ncol = 5; let max_nnz = 13; @@ -1114,7 +1114,7 @@ impl Samples { shuffle_coo_entries: bool, duplicate_coo_entries: bool, ) -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = None; + let sym = Sym::No; let nrow = 1; let ncol = 2; let max_nnz = 10; @@ -1174,7 +1174,7 @@ impl Samples { /// └ ┘ /// ``` pub fn rectangular_1x7() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = None; + let sym = Sym::No; let nrow = 1; let ncol = 7; let mut coo = CooMatrix::new(nrow, ncol, 4, sym).unwrap(); @@ -1221,7 +1221,7 @@ impl Samples { /// └ ┘ /// ``` pub fn rectangular_7x1() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = None; + let sym = Sym::No; let nrow = 7; let ncol = 1; let mut coo = CooMatrix::new(nrow, ncol, 3, sym).unwrap(); @@ -1261,7 +1261,7 @@ impl Samples { /// 15 -6 . 3 /// ``` pub fn rectangular_3x4() -> (CooMatrix, CscMatrix, CsrMatrix, f64) { - let sym = None; + let sym = Sym::No; let nrow = 3; let ncol = 4; let mut coo = CooMatrix::new(nrow, ncol, 9, sym).unwrap(); @@ -1316,7 +1316,7 @@ impl Samples { /// 1 . . /// ``` pub fn complex_rectangular_4x3() -> (ComplexCooMatrix, ComplexCscMatrix, ComplexCsrMatrix, f64) { - let sym = None; + let sym = Sym::No; let nrow = 4; let ncol = 3; let mut coo = ComplexCooMatrix::new(nrow, ncol, 7, sym).unwrap(); @@ -1417,7 +1417,7 @@ mod tests { assert_eq!(csc.row_indices.len(), nnz); assert_eq!(csc.values.len(), nnz); // CSC vs COO - assert_eq!(csc.symmetry, coo.symmetry); + assert_eq!(csc.symmetric, coo.symmetric); assert_eq!(csc.nrow, coo.nrow); assert_eq!(csc.ncol, coo.ncol); // CSR @@ -1430,7 +1430,7 @@ mod tests { assert_eq!(csr.col_indices.len(), nnz); assert_eq!(csr.values.len(), nnz); // CSR vs COO - assert_eq!(csr.symmetry, coo.symmetry); + assert_eq!(csr.symmetric, coo.symmetric); assert_eq!(csr.nrow, coo.nrow); assert_eq!(csr.ncol, coo.ncol); } @@ -1475,7 +1475,7 @@ mod tests { let a = Matrix::from(correct); let mut ai = Matrix::new(3, 3); let correct_det = mat_inverse(&mut ai, &a).unwrap(); - let (coo, csc, csr, det) = Samples::positive_definite_3x3(); + let (coo, csc, csr, det) = Samples::positive_definite_3x3_lower(); approx_eq(det, correct_det, 1e-15); mat_approx_eq(&coo.as_dense(), correct, 1e-15); mat_approx_eq(&csc.as_dense(), correct, 1e-15); diff --git a/russell_sparse/src/solver_mumps.rs b/russell_sparse/src/solver_mumps.rs index d92df2e6..f225b040 100644 --- a/russell_sparse/src/solver_mumps.rs +++ b/russell_sparse/src/solver_mumps.rs @@ -123,8 +123,8 @@ pub struct SolverMUMPS { /// Indicates whether the sparse matrix has been factorized or not factorized: bool, - /// Holds the symmetry type used in the initialize - initialized_symmetry: Sym, + /// Holds the symmetric flag saved in initialize + initialized_sym: Sym, /// Holds the matrix dimension saved in initialize initialized_ndim: usize, @@ -199,7 +199,7 @@ impl SolverMUMPS { solver, initialized: false, factorized: false, - initialized_symmetry: Sym::No, + initialized_sym: Sym::No, initialized_ndim: 0, initialized_nnz: 0, effective_ordering: -1, @@ -227,7 +227,7 @@ impl LinSolTrait for SolverMUMPS { /// /// * `mat` -- the coefficient matrix A (one-base **COO** only, not CSC and not CSR). /// Also, the matrix must be square (`nrow = ncol`) and, if symmetric, - /// the symmetry/storage must [crate::Storage::Lower]. + /// the symmetric flag must be [Sym::YesLower] /// * `params` -- configuration parameters; None => use default /// /// # Notes @@ -239,24 +239,16 @@ impl LinSolTrait for SolverMUMPS { /// 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. - /// 4. For symmetric matrices, `MUMPS` requires that the symmetry/storage be [crate::Storage::Lower]. + /// 4. For symmetric matrices, `MUMPS` requires [Sym::YesLower]. /// 5. The COO matrix must be one-based. fn factorize(&mut self, mat: &mut SparseMatrix, params: Option) -> Result<(), StrError> { // get COO matrix let coo = mat.get_coo()?; - // check the COO matrix - if coo.nrow != coo.ncol { - return Err("the COO matrix must be square"); - } - if coo.nnz < 1 { - return Err("the COO matrix must have at least one non-zero value"); - } - - // check already initialized data + // check if self.initialized { - if coo.symmetry != self.initialized_symmetry { - return Err("subsequent factorizations must use the same matrix (symmetry differs)"); + if coo.symmetric != self.initialized_sym { + return Err("subsequent factorizations must use the same matrix (symmetric differs)"); } if coo.nrow != self.initialized_ndim { return Err("subsequent factorizations must use the same matrix (ndim differs)"); @@ -265,7 +257,16 @@ impl LinSolTrait for SolverMUMPS { return Err("subsequent factorizations must use the same matrix (nnz differs)"); } } else { - self.initialized_symmetry = coo.symmetry; + if coo.nrow != coo.ncol { + return Err("the COO matrix must be square"); + } + if coo.nnz < 1 { + return Err("the COO matrix must have at least one non-zero value"); + } + if coo.symmetric == Sym::YesFull || coo.symmetric == Sym::YesUpper { + return Err("MUMPS requires Sym::YesLower for symmetric matrices"); + } + self.initialized_sym = coo.symmetric; self.initialized_ndim = coo.nrow; self.initialized_nnz = coo.nnz; self.fortran_indices_i = vec![0; coo.nnz]; @@ -304,10 +305,9 @@ impl LinSolTrait for SolverMUMPS { let compute_determinant = if par.compute_determinant { 1 } else { 0 }; let verbose = if par.verbose { 1 } else { 0 }; - // extract the symmetry flags and check the storage type - let (general_symmetric, positive_definite) = coo.symmetry.status(true, false)?; - // matrix config + let general_symmetric = if coo.symmetric == Sym::YesLower { 1 } else { 0 }; + let positive_definite = if par.positive_definite { 1 } else { 0 }; let ndim = to_i32(coo.nrow); let nnz = to_i32(coo.nnz); @@ -377,7 +377,7 @@ impl LinSolTrait for SolverMUMPS { /// /// # Input /// - /// * `mat` -- the coefficient matrix A; must be square and, if symmetric, [crate::Storage::Lower]. + /// * `mat` -- the coefficient matrix A; must be square and, if symmetric, [Sym::YesLower]. /// * `rhs` -- the right-hand side vector with know values an dimension equal to mat.nrow /// * `verbose` -- shows messages /// @@ -393,8 +393,8 @@ impl LinSolTrait for SolverMUMPS { // check already factorized data let (nrow, ncol, nnz, sym) = coo.get_info(); - if sym != self.initialized_symmetry { - return Err("solve must use the same matrix (symmetry differs)"); + if sym != self.initialized_sym { + return Err("solve must use the same matrix (symmetric differs)"); } if nrow != self.initialized_ndim || ncol != self.initialized_ndim { return Err("solve must use the same matrix (ndim differs)"); @@ -616,7 +616,7 @@ pub(crate) fn handle_mumps_error_code(err: i32) -> StrError { #[cfg(test)] mod tests { use super::*; - use crate::{CooMatrix, LinSolParams, LinSolTrait, Ordering, Samples, Scaling, SparseMatrix, Storage, Sym}; + use crate::{CooMatrix, LinSolParams, LinSolTrait, Ordering, Samples, Scaling, SparseMatrix, Sym}; use russell_lab::{approx_eq, vec_approx_eq, Vector}; use serial_test::serial; @@ -645,7 +645,7 @@ mod tests { solver.factorize(&mut mat, None).err(), Some("the COO matrix must be square") ); - let coo = CooMatrix::new(1, 1, 1, None).unwrap(); + let coo = CooMatrix::new(1, 1, 1, Sym::No).unwrap(); let mut mat = SparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), @@ -655,27 +655,27 @@ mod tests { let mut mat = SparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), - Some("if the matrix is general symmetric, the required storage is lower triangular") + Some("MUMPS requires Sym::YesLower for symmetric matrices") ); // check already factorized data - let mut coo = CooMatrix::new(2, 2, 2, None).unwrap(); + let mut coo = CooMatrix::new(2, 2, 2, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); // ... factorize once => OK solver.factorize(&mut mat, None).unwrap(); - // ... change matrix (symmetry) - let mut coo = CooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); + // ... change matrix (symmetric) + let mut coo = CooMatrix::new(2, 2, 2, Sym::YesFull).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), - Some("subsequent factorizations must use the same matrix (symmetry differs)") + Some("subsequent factorizations must use the same matrix (symmetric differs)") ); // ... change matrix (ndim) - let mut coo = CooMatrix::new(1, 1, 1, None).unwrap(); + let mut coo = CooMatrix::new(1, 1, 1, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); assert_eq!( @@ -683,7 +683,7 @@ mod tests { Some("subsequent factorizations must use the same matrix (ndim differs)") ); // ... change matrix (nnz) - let mut coo = CooMatrix::new(2, 2, 1, None).unwrap(); + let mut coo = CooMatrix::new(2, 2, 1, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); assert_eq!( @@ -695,7 +695,7 @@ mod tests { #[test] #[serial] fn factorize_fails_on_singular_matrix() { - let mut mat_singular = SparseMatrix::new_coo(5, 5, 2, None).unwrap(); + let mut mat_singular = SparseMatrix::new_coo(5, 5, 2, Sym::No).unwrap(); mat_singular.put(0, 0, 1.0).unwrap(); mat_singular.put(4, 4, 1.0).unwrap(); let mut solver = SolverMUMPS::new().unwrap(); @@ -708,7 +708,7 @@ mod tests { #[test] #[serial] fn solve_handles_errors() { - let mut coo = CooMatrix::new(2, 2, 2, None).unwrap(); + let mut coo = CooMatrix::new(2, 2, 2, Sym::No).unwrap(); coo.put(0, 0, 123.0).unwrap(); coo.put(1, 1, 456.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); @@ -732,19 +732,19 @@ mod tests { solver.solve(&mut x, &mut mat, &rhs, false), Err("the dimension of the right-hand side vector is incorrect") ); - // wrong symmetry + // wrong symmetric let rhs = Vector::new(2); - let mut coo_wrong = CooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); + let mut coo_wrong = CooMatrix::new(2, 2, 2, Sym::YesFull).unwrap(); coo_wrong.put(0, 0, 123.0).unwrap(); coo_wrong.put(1, 1, 456.0).unwrap(); let mut mat_wrong = SparseMatrix::from_coo(coo_wrong); mat_wrong.get_csc_or_from_coo().unwrap(); // make sure to convert to CSC (because we're not calling factorize on this wrong matrix) assert_eq!( solver.solve(&mut x, &mut mat_wrong, &rhs, false), - Err("solve must use the same matrix (symmetry differs)") + Err("solve must use the same matrix (symmetric differs)") ); // wrong ndim - let mut coo_wrong = CooMatrix::new(1, 1, 1, None).unwrap(); + let mut coo_wrong = CooMatrix::new(1, 1, 1, Sym::No).unwrap(); coo_wrong.put(0, 0, 123.0).unwrap(); let mut mat_wrong = SparseMatrix::from_coo(coo_wrong); mat_wrong.get_csc_or_from_coo().unwrap(); // make sure to convert to CSC (because we're not calling factorize on this wrong matrix) @@ -753,7 +753,7 @@ mod tests { Err("solve must use the same matrix (ndim differs)") ); // wrong nnz - let mut coo_wrong = CooMatrix::new(2, 2, 3, None).unwrap(); + let mut coo_wrong = CooMatrix::new(2, 2, 3, Sym::No).unwrap(); coo_wrong.put(0, 0, 123.0).unwrap(); coo_wrong.put(1, 1, 123.0).unwrap(); coo_wrong.put(0, 1, 100.0).unwrap(); @@ -829,6 +829,36 @@ mod tests { vec_approx_eq(x.as_data(), x_correct, 1e-10); } + #[test] + #[serial] + fn solve_works_symmetric() { + // allocate x and rhs + let mut x = Vector::new(5); + let rhs = Vector::from(&[1.0, 2.0, 3.0, 4.0, 5.0]); + let x_correct = &[-979.0 / 3.0, 983.0, 1961.0 / 12.0, 398.0, 123.0 / 2.0]; + + // allocate a new solver + let mut solver = SolverMUMPS::new().unwrap(); + assert!(!solver.factorized); + + // sample matrix + let (coo, _, _, _) = Samples::mkl_symmetric_5x5_lower(false, true); + let mut mat = SparseMatrix::from_coo(coo); + + // factorize works + solver.factorize(&mut mat, None).unwrap(); + assert!(solver.factorized); + + // solve works + solver.solve(&mut x, &mat, &rhs, false).unwrap(); + vec_approx_eq(x.as_data(), x_correct, 1e-11); + + // calling solve again works + let mut x_again = Vector::new(5); + solver.solve(&mut x_again, &mat, &rhs, false).unwrap(); + vec_approx_eq(x_again.as_data(), x_correct, 1e-11); + } + #[test] fn ordering_and_scaling_works() { assert_eq!(mumps_ordering(Ordering::Amd), MUMPS_ORDERING_AMD); diff --git a/russell_sparse/src/solver_umfpack.rs b/russell_sparse/src/solver_umfpack.rs index a003244d..3c0e568f 100644 --- a/russell_sparse/src/solver_umfpack.rs +++ b/russell_sparse/src/solver_umfpack.rs @@ -74,8 +74,8 @@ pub struct SolverUMFPACK { /// Indicates whether the sparse matrix has been factorized or not factorized: bool, - /// Holds the symmetry type used in initialize - initialized_symmetry: Sym, + /// Holds the symmetric flag saved in initialize + initialized_sym: Sym, /// Holds the matrix dimension saved in initialize initialized_ndim: usize, @@ -139,7 +139,7 @@ impl SolverUMFPACK { solver, initialized: false, factorized: false, - initialized_symmetry: Sym::No, + initialized_sym: Sym::No, initialized_ndim: 0, initialized_nnz: 0, effective_strategy: -1, @@ -164,7 +164,7 @@ impl LinSolTrait for SolverUMFPACK { /// /// * `mat` -- the coefficient matrix A (**COO** or **CSC**, but not CSR). /// Also, the matrix must be square (`nrow = ncol`) and, if symmetric, - /// the symmetry/storage must [crate::Storage::Full]. + /// the symmetric flag must be [Sym::YesFull] /// * `params` -- configuration parameters; None => use default /// /// # Notes @@ -176,24 +176,16 @@ impl LinSolTrait for SolverUMFPACK { /// 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. - /// 4. For symmetric matrices, `UMFPACK` requires that the symmetry/storage be [crate::Storage::Full]. + /// 4. For symmetric matrices, `UMFPACK` requires [Sym::YesFull] fn factorize(&mut self, mat: &mut SparseMatrix, params: Option) -> Result<(), StrError> { // get CSC matrix // (or convert from COO if CSC is not available and COO is available) let csc = mat.get_csc_or_from_coo()?; - // check CSC matrix - if csc.nrow != csc.ncol { - return Err("the matrix must be square"); - } - if csc.symmetry.triangular() { - return Err("for UMFPACK, the matrix must not be triangular"); - } - - // check already initialized data + // check if self.initialized { - if csc.symmetry != self.initialized_symmetry { - return Err("subsequent factorizations must use the same matrix (symmetry differs)"); + if csc.symmetric != self.initialized_sym { + return Err("subsequent factorizations must use the same matrix (symmetric differs)"); } if csc.nrow != self.initialized_ndim { return Err("subsequent factorizations must use the same matrix (ndim differs)"); @@ -202,7 +194,13 @@ impl LinSolTrait for SolverUMFPACK { return Err("subsequent factorizations must use the same matrix (nnz differs)"); } } else { - self.initialized_symmetry = csc.symmetry; + if csc.nrow != csc.ncol { + return Err("the matrix must be square"); + } + if csc.symmetric == Sym::YesLower || csc.symmetric == Sym::YesUpper { + return Err("UMFPACK requires Sym::YesFull for symmetric matrices"); + } + self.initialized_sym = csc.symmetric; self.initialized_ndim = csc.nrow; self.initialized_nnz = csc.col_pointers[csc.ncol] as usize; } @@ -288,7 +286,7 @@ impl LinSolTrait for SolverUMFPACK { /// /// # Input /// - /// * `mat` -- the coefficient matrix A; must be square and, if symmetric, [crate::Storage::Full]. + /// * `mat` -- the coefficient matrix A; must be square and, if symmetric, [Sym::YesFull]. /// * `rhs` -- the right-hand side vector with know values an dimension equal to mat.nrow /// * `verbose` -- shows messages /// @@ -305,8 +303,8 @@ impl LinSolTrait for SolverUMFPACK { // check already factorized data let (nrow, ncol, nnz, sym) = csc.get_info(); - if sym != self.initialized_symmetry { - return Err("solve must use the same matrix (symmetry differs)"); + if sym != self.initialized_sym { + return Err("solve must use the same matrix (symmetric differs)"); } if nrow != self.initialized_ndim || ncol != self.initialized_ndim { return Err("solve must use the same matrix (ndim differs)"); @@ -465,7 +463,7 @@ pub(crate) fn handle_umfpack_error_code(err: i32) -> StrError { #[cfg(test)] mod tests { use super::*; - use crate::{CooMatrix, LinSolParams, LinSolTrait, Ordering, Samples, Scaling, SparseMatrix, Storage, Sym}; + use crate::{CooMatrix, LinSolParams, LinSolTrait, Ordering, Samples, Scaling, SparseMatrix, Sym}; use russell_lab::{approx_eq, vec_approx_eq, Vector}; #[test] @@ -481,7 +479,7 @@ mod tests { assert!(!solver.factorized); // COO to CSC errors - let coo = CooMatrix::new(1, 1, 1, None).unwrap(); + let coo = CooMatrix::new(1, 1, 1, Sym::No).unwrap(); let mut mat = SparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), @@ -499,27 +497,27 @@ mod tests { let mut mat = SparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), - Some("for UMFPACK, the matrix must not be triangular") + Some("UMFPACK requires Sym::YesFull for symmetric matrices") ); // check already factorized data - let mut coo = CooMatrix::new(2, 2, 2, None).unwrap(); + let mut coo = CooMatrix::new(2, 2, 2, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); // ... factorize once => OK solver.factorize(&mut mat, None).unwrap(); - // ... change matrix (symmetry) - let mut coo = CooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); + // ... change matrix (symmetric) + let mut coo = CooMatrix::new(2, 2, 2, Sym::YesFull).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(1, 1, 2.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); assert_eq!( solver.factorize(&mut mat, None).err(), - Some("subsequent factorizations must use the same matrix (symmetry differs)") + Some("subsequent factorizations must use the same matrix (symmetric differs)") ); // ... change matrix (ndim) - let mut coo = CooMatrix::new(1, 1, 1, None).unwrap(); + let mut coo = CooMatrix::new(1, 1, 1, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); assert_eq!( @@ -527,7 +525,7 @@ mod tests { Some("subsequent factorizations must use the same matrix (ndim differs)") ); // ... change matrix (nnz) - let mut coo = CooMatrix::new(2, 2, 1, None).unwrap(); + let mut coo = CooMatrix::new(2, 2, 1, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); assert_eq!( @@ -566,7 +564,7 @@ mod tests { #[test] fn factorize_fails_on_singular_matrix() { let mut solver = SolverUMFPACK::new().unwrap(); - let mut coo = CooMatrix::new(2, 2, 2, None).unwrap(); + let mut coo = CooMatrix::new(2, 2, 2, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(1, 1, 0.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); @@ -575,7 +573,7 @@ mod tests { #[test] fn solve_handles_errors() { - let mut coo = CooMatrix::new(2, 2, 2, None).unwrap(); + let mut coo = CooMatrix::new(2, 2, 2, Sym::No).unwrap(); coo.put(0, 0, 123.0).unwrap(); coo.put(1, 1, 456.0).unwrap(); let mut mat = SparseMatrix::from_coo(coo); @@ -599,19 +597,19 @@ mod tests { solver.solve(&mut x, &mut mat, &rhs, false), Err("the dimension of the right-hand side vector is incorrect") ); - // wrong symmetry + // wrong symmetric let rhs = Vector::new(2); - let mut coo_wrong = CooMatrix::new(2, 2, 2, Some(Sym::General(Storage::Full))).unwrap(); + let mut coo_wrong = CooMatrix::new(2, 2, 2, Sym::YesFull).unwrap(); coo_wrong.put(0, 0, 123.0).unwrap(); coo_wrong.put(1, 1, 456.0).unwrap(); let mut mat_wrong = SparseMatrix::from_coo(coo_wrong); mat_wrong.get_csc_or_from_coo().unwrap(); // make sure to convert to CSC (because we're not calling factorize on this wrong matrix) assert_eq!( solver.solve(&mut x, &mut mat_wrong, &rhs, false), - Err("solve must use the same matrix (symmetry differs)") + Err("solve must use the same matrix (symmetric differs)") ); // wrong ndim - let mut coo_wrong = CooMatrix::new(1, 1, 1, None).unwrap(); + let mut coo_wrong = CooMatrix::new(1, 1, 1, Sym::No).unwrap(); coo_wrong.put(0, 0, 123.0).unwrap(); let mut mat_wrong = SparseMatrix::from_coo(coo_wrong); mat_wrong.get_csc_or_from_coo().unwrap(); // make sure to convert to CSC (because we're not calling factorize on this wrong matrix) @@ -620,7 +618,7 @@ mod tests { Err("solve must use the same matrix (ndim differs)") ); // wrong nnz - let mut coo_wrong = CooMatrix::new(2, 2, 3, None).unwrap(); + let mut coo_wrong = CooMatrix::new(2, 2, 3, Sym::No).unwrap(); coo_wrong.put(0, 0, 123.0).unwrap(); coo_wrong.put(1, 1, 123.0).unwrap(); coo_wrong.put(0, 1, 100.0).unwrap(); @@ -656,6 +654,30 @@ mod tests { assert_eq!(stats.output.effective_scaling, "Sum"); } + #[test] + fn solve_works_symmetric() { + let mut solver = SolverUMFPACK::new().unwrap(); + let (coo, _, _, _) = Samples::mkl_symmetric_5x5_full(); + let mut mat = SparseMatrix::from_coo(coo); + let mut x = Vector::new(5); + let rhs = Vector::from(&[1.0, 2.0, 3.0, 4.0, 5.0]); + let x_correct = &[-979.0 / 3.0, 983.0, 1961.0 / 12.0, 398.0, 123.0 / 2.0]; + solver.factorize(&mut mat, None).unwrap(); + solver.solve(&mut x, &mut mat, &rhs, false).unwrap(); + vec_approx_eq(x.as_data(), x_correct, 1e-10); + + // calling solve again works + let mut x_again = Vector::new(5); + solver.solve(&mut x_again, &mut mat, &rhs, false).unwrap(); + vec_approx_eq(x_again.as_data(), x_correct, 1e-10); + + // update stats + let mut stats = StatsLinSol::new(); + solver.update_stats(&mut stats); + assert_eq!(stats.output.effective_ordering, "Amd"); + assert_eq!(stats.output.effective_scaling, "Sum"); + } + #[test] fn ordering_and_scaling_works() { assert_eq!(umfpack_ordering(Ordering::Amd), UMFPACK_ORDERING_AMD); diff --git a/russell_sparse/src/sparse_matrix.rs b/russell_sparse/src/sparse_matrix.rs index e25fe539..47772213 100644 --- a/russell_sparse/src/sparse_matrix.rs +++ b/russell_sparse/src/sparse_matrix.rs @@ -58,10 +58,11 @@ where /// * `ncol` -- (≥ 1) Is the number of columns of the sparse matrix (must be fit i32) /// * `max_nnz` -- (≥ 1) Maximum number of entries ≥ nnz (number of non-zeros), /// including entries with repeated indices. (must be fit i32) - /// * `sym` -- Defines the symmetry/storage, if any - pub fn new_coo(nrow: usize, ncol: usize, max_nnz: usize, sym: Option) -> Result { + /// * `symmetric` -- indicates whether the matrix is symmetric or not. + /// If symmetric, indicates the representation too. + pub fn new_coo(nrow: usize, ncol: usize, max_nnz: usize, symmetric: Sym) -> Result { Ok(NumSparseMatrix { - coo: Some(NumCooMatrix::new(nrow, ncol, max_nnz, sym)?), + coo: Some(NumCooMatrix::new(nrow, ncol, max_nnz, symmetric)?), csc: None, csr: None, }) @@ -79,6 +80,8 @@ where /// to the number of non-zero values (sorted) /// * `row_indices` -- (len = nnz) row indices (sorted) /// * `values` -- the non-zero components of the matrix + /// * `symmetric` -- indicates whether the matrix is symmetric or not. + /// If symmetric, indicates the representation too. /// /// The following conditions must be satisfied (nnz is the number of non-zeros /// and nnz_dup is the number of non-zeros with possible duplicates): @@ -98,7 +101,7 @@ where col_pointers: Vec, row_indices: Vec, values: Vec, - symmetry: Option, + symmetric: Sym, ) -> Result { Ok(NumSparseMatrix { coo: None, @@ -108,7 +111,7 @@ where col_pointers, row_indices, values, - symmetry, + symmetric, )?), csr: None, }) @@ -126,6 +129,8 @@ where /// to the number of non-zero values (sorted) /// * `col_indices` -- (len = nnz) column indices (sorted) /// * `values` -- the non-zero components of the matrix + /// * `symmetric` -- indicates whether the matrix is symmetric or not. + /// If symmetric, indicates the representation too. /// /// The following conditions must be satisfied (nnz is the number of non-zeros /// and nnz_dup is the number of non-zeros with possible duplicates): @@ -145,7 +150,7 @@ where row_pointers: Vec, col_indices: Vec, values: Vec, - symmetry: Option, + symmetric: Sym, ) -> Result { Ok(NumSparseMatrix { coo: None, @@ -156,7 +161,7 @@ where row_pointers, col_indices, values, - symmetry, + symmetric, )?), }) } @@ -443,21 +448,21 @@ mod tests { #[test] fn new_functions_work() { // COO - NumSparseMatrix::::new_coo(1, 1, 1, None).unwrap(); + NumSparseMatrix::::new_coo(1, 1, 1, Sym::No).unwrap(); assert_eq!( - NumSparseMatrix::::new_coo(0, 1, 1, None).err(), + NumSparseMatrix::::new_coo(0, 1, 1, Sym::No).err(), Some("nrow must be ≥ 1") ); // CSC - NumSparseMatrix::::new_csc(1, 1, vec![0, 1], vec![0], vec![0.0], None).unwrap(); + NumSparseMatrix::::new_csc(1, 1, vec![0, 1], vec![0], vec![0.0], Sym::No).unwrap(); assert_eq!( - NumSparseMatrix::::new_csc(0, 1, vec![0, 1], vec![0], vec![0.0], None).err(), + NumSparseMatrix::::new_csc(0, 1, vec![0, 1], vec![0], vec![0.0], Sym::No).err(), Some("nrow must be ≥ 1") ); // CSR - NumSparseMatrix::::new_csr(1, 1, vec![0, 1], vec![0], vec![0.0], None).unwrap(); + NumSparseMatrix::::new_csr(1, 1, vec![0, 1], vec![0], vec![0.0], Sym::No).unwrap(); assert_eq!( - NumSparseMatrix::::new_csr(0, 1, vec![0, 1], vec![0], vec![0.0], None).err(), + NumSparseMatrix::::new_csr(0, 1, vec![0, 1], vec![0], vec![0.0], Sym::No).err(), Some("nrow must be ≥ 1") ); } @@ -514,8 +519,8 @@ mod tests { fn setters_work() { // test matrices let (coo, csc, csr, _) = Samples::rectangular_1x2(false, false); - let mut other = NumSparseMatrix::::new_coo(1, 1, 1, None).unwrap(); - let mut wrong = NumSparseMatrix::::new_coo(1, 1, 3, None).unwrap(); + let mut other = NumSparseMatrix::::new_coo(1, 1, 1, Sym::No).unwrap(); + let mut wrong = NumSparseMatrix::::new_coo(1, 1, 3, Sym::No).unwrap(); other.put(0, 0, 2.0).unwrap(); wrong.put(0, 0, 1.0).unwrap(); wrong.put(0, 0, 2.0).unwrap(); @@ -525,7 +530,7 @@ mod tests { assert_eq!(coo_mat.get_coo_mut().unwrap().get_info(), (1, 2, 2, Sym::No)); assert_eq!(coo_mat.get_csc_mut().err(), Some("CSC matrix is not available")); assert_eq!(coo_mat.get_csr_mut().err(), Some("CSR matrix is not available")); - let mut empty = NumSparseMatrix::::new_coo(1, 1, 1, None).unwrap(); + let mut empty = NumSparseMatrix::::new_coo(1, 1, 1, Sym::No).unwrap(); assert_eq!(empty.get_csc_or_from_coo().err(), Some("COO to CSC requires nnz > 0")); assert_eq!(empty.get_csr_or_from_coo().err(), Some("COO to CSR requires nnz > 0")); // CSC @@ -581,7 +586,7 @@ mod tests { Some("COO matrix is not available to augment") ); // COO - let mut coo = NumSparseMatrix::::new_coo(2, 2, 1, None).unwrap(); + let mut coo = NumSparseMatrix::::new_coo(2, 2, 1, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); assert_eq!( coo.put(1, 1, 2.0).err(), @@ -590,7 +595,7 @@ mod tests { coo.reset().unwrap(); coo.put(1, 1, 2.0).unwrap(); // COO (assign) - let mut this = NumSparseMatrix::::new_coo(1, 1, 1, None).unwrap(); + let mut this = NumSparseMatrix::::new_coo(1, 1, 1, Sym::No).unwrap(); this.put(0, 0, 8000.0).unwrap(); this.assign(4.0, &other).unwrap(); assert_eq!( @@ -605,7 +610,7 @@ mod tests { ); assert_eq!(this.assign(2.0, &csc_mat).err(), Some("COO matrix is not available")); // COO (augment) - let mut this = NumSparseMatrix::::new_coo(1, 1, 1 + 1, None).unwrap(); + let mut this = NumSparseMatrix::::new_coo(1, 1, 1 + 1, Sym::No).unwrap(); this.put(0, 0, 100.0).unwrap(); this.augment(4.0, &other).unwrap(); assert_eq!( @@ -662,7 +667,7 @@ mod tests { let json = serde_json::to_string(&mat).unwrap(); assert_eq!( json, - r#"{"coo":{"symmetry":"No","nrow":1,"ncol":1,"nnz":1,"max_nnz":1,"indices_i":[0],"indices_j":[0],"values":[123.0]},"csc":null,"csr":null}"# + r#"{"coo":{"symmetric":"No","nrow":1,"ncol":1,"nnz":1,"max_nnz":1,"indices_i":[0],"indices_j":[0],"values":[123.0]},"csc":null,"csr":null}"# ); let from_json: NumSparseMatrix = serde_json::from_str(&json).unwrap(); let (json_nrow, json_ncol, json_nnz, json_sym) = from_json.get_coo().unwrap().get_info(); diff --git a/russell_sparse/src/verify_lin_sys.rs b/russell_sparse/src/verify_lin_sys.rs index dca60b51..9bec6ee2 100644 --- a/russell_sparse/src/verify_lin_sys.rs +++ b/russell_sparse/src/verify_lin_sys.rs @@ -32,7 +32,7 @@ impl VerifyLinSys { /// fn main() -> Result<(), StrError> { /// // set sparse matrix (3 x 3) with 4 non-zeros /// let (nrow, ncol, nnz) = (3, 3, 4); - /// let mut coo = SparseMatrix::new_coo(nrow, ncol, nnz, None)?; + /// let mut coo = SparseMatrix::new_coo(nrow, ncol, nnz, Sym::No)?; /// coo.put(0, 0, 1.0)?; /// coo.put(0, 2, 4.0)?; /// coo.put(1, 1, 2.0)?; @@ -152,14 +152,14 @@ impl VerifyLinSys { #[cfg(test)] mod tests { use super::VerifyLinSys; - use crate::{ComplexSparseMatrix, Samples, SparseMatrix}; + use crate::{ComplexSparseMatrix, Samples, SparseMatrix, Sym}; use num_complex::Complex64; use russell_lab::{approx_eq, cpx, ComplexVector, Vector}; #[test] fn from_captures_errors() { // real - let coo = SparseMatrix::new_coo(2, 1, 1, None).unwrap(); + let coo = SparseMatrix::new_coo(2, 1, 1, Sym::No).unwrap(); let x = Vector::new(1); let rhs = Vector::new(2); assert_eq!(VerifyLinSys::from(&coo, &x, &rhs).err(), Some("matrix is empty")); @@ -174,7 +174,7 @@ mod tests { Some("rhs.dim() must be equal to nrow") ); // complex - let coo = ComplexSparseMatrix::new_coo(2, 1, 1, None).unwrap(); + let coo = ComplexSparseMatrix::new_coo(2, 1, 1, Sym::No).unwrap(); let x = ComplexVector::new(1); let rhs = ComplexVector::new(2); assert_eq!( @@ -198,7 +198,7 @@ mod tests { // 1 3 -2 // 3 5 6 // 2 4 3 - let mut coo = SparseMatrix::new_coo(3, 3, 9, None).unwrap(); + let mut coo = SparseMatrix::new_coo(3, 3, 9, Sym::No).unwrap(); coo.put(0, 0, 1.0).unwrap(); coo.put(0, 1, 3.0).unwrap(); coo.put(0, 2, -2.0).unwrap(); diff --git a/russell_sparse/src/write_matrix_market.rs b/russell_sparse/src/write_matrix_market.rs index 36c98a07..22f6f1af 100644 --- a/russell_sparse/src/write_matrix_market.rs +++ b/russell_sparse/src/write_matrix_market.rs @@ -182,7 +182,7 @@ mod tests { // 2 -1 2 sym // -1 2 -1 => -1 2 // -1 2 -1 2 - let (_, csc, _, _) = Samples::positive_definite_3x3(); + let (_, csc, _, _) = Samples::positive_definite_3x3_lower(); let full_path = "/tmp/russell_sparse/test_write_matrix_market_csc.mtx"; csc_write_matrix_market(&csc, full_path, false).unwrap(); let contents = fs::read_to_string(full_path).map_err(|_| "cannot open file").unwrap(); @@ -229,7 +229,7 @@ mod tests { // 2 -1 2 sym // -1 2 -1 => -1 2 // -1 2 -1 2 - let (_, _, csr, _) = Samples::positive_definite_3x3(); + let (_, _, csr, _) = Samples::positive_definite_3x3_lower(); let full_path = "/tmp/russell_sparse/test_write_matrix_market_csr.mtx"; csr_write_matrix_market(&csr, full_path, false).unwrap(); let contents = fs::read_to_string(full_path).map_err(|_| "cannot open file").unwrap(); diff --git a/russell_sparse/tests/test_complex_coo_matrix.rs b/russell_sparse/tests/test_complex_coo_matrix.rs index 6943e2f1..dedf4b55 100644 --- a/russell_sparse/tests/test_complex_coo_matrix.rs +++ b/russell_sparse/tests/test_complex_coo_matrix.rs @@ -5,8 +5,7 @@ use russell_sparse::StrError; #[test] fn test_complex_coo_matrix() -> Result<(), StrError> { - let sym = Some(Sym::new_general_lower()); - let mut coo = ComplexCooMatrix::new(3, 3, 4, sym)?; + let mut coo = ComplexCooMatrix::new(3, 3, 4, Sym::YesLower)?; coo.put(0, 0, cpx!(1.0, 0.1))?; coo.put(1, 0, cpx!(2.0, 0.2))?; coo.put(1, 1, cpx!(3.0, 0.3))?; diff --git a/russell_sparse/tests/test_complex_mumps.rs b/russell_sparse/tests/test_complex_mumps.rs index 0f881481..4536d8da 100644 --- a/russell_sparse/tests/test_complex_mumps.rs +++ b/russell_sparse/tests/test_complex_mumps.rs @@ -9,7 +9,7 @@ use serial_test::serial; fn test_complex_mumps() -> Result<(), StrError> { let n = 10; let d = (n as f64) / 10.0; - let mut coo = ComplexCooMatrix::new(n, n, n, None)?; + let mut coo = ComplexCooMatrix::new(n, n, n, Sym::No)?; let mut x_correct = ComplexVector::new(n); let mut rhs = ComplexVector::new(n); for k in 0..n { diff --git a/russell_sparse/tests/test_complex_umfpack.rs b/russell_sparse/tests/test_complex_umfpack.rs index 07dc2063..1cc34cea 100644 --- a/russell_sparse/tests/test_complex_umfpack.rs +++ b/russell_sparse/tests/test_complex_umfpack.rs @@ -7,7 +7,7 @@ use russell_sparse::StrError; fn test_complex_umfpack() -> Result<(), StrError> { let n = 10; let d = (n as f64) / 10.0; - let mut coo = ComplexCooMatrix::new(n, n, n, None)?; + let mut coo = ComplexCooMatrix::new(n, n, n, Sym::No)?; let mut x_correct = ComplexVector::new(n); let mut rhs = ComplexVector::new(n); for k in 0..n { diff --git a/russell_sparse/tests/test_mumps.rs b/russell_sparse/tests/test_mumps.rs index 58d770d1..eedc9a2e 100644 --- a/russell_sparse/tests/test_mumps.rs +++ b/russell_sparse/tests/test_mumps.rs @@ -8,7 +8,7 @@ use serial_test::serial; fn test_complex_umfpack() -> Result<(), StrError> { let n = 10; let d = (n as f64) / 10.0; - let mut coo = CooMatrix::new(n, n, n, None)?; + let mut coo = CooMatrix::new(n, n, n, Sym::No)?; let mut x_correct = Vector::new(n); let mut rhs = Vector::new(n); for k in 0..n { diff --git a/russell_sparse/tests/test_nonlinear_system.rs b/russell_sparse/tests/test_nonlinear_system.rs index 39d60289..036006fc 100644 --- a/russell_sparse/tests/test_nonlinear_system.rs +++ b/russell_sparse/tests/test_nonlinear_system.rs @@ -67,7 +67,7 @@ fn check_jacobian() { } } let nnz = neq * neq; - let mut jj_tri = CooMatrix::new(neq, neq, nnz, None).unwrap(); + let mut jj_tri = CooMatrix::new(neq, neq, nnz, Sym::No).unwrap(); calc_jacobian(&mut jj_tri, &uu).unwrap(); let mut jj_ana = Matrix::new(neq, neq); jj_tri.to_dense(&mut jj_ana).unwrap(); @@ -77,7 +77,7 @@ fn check_jacobian() { fn solve_nonlinear_system(genie: Genie) -> Result<(), StrError> { let (neq, nnz) = (4, 16); let mut solver = LinSolver::new(genie)?; - let mut jj = SparseMatrix::new_coo(neq, neq, nnz, None).unwrap(); + let mut jj = SparseMatrix::new_coo(neq, neq, nnz, Sym::No).unwrap(); let mut rr = Vector::new(neq); let mut uu = Vector::from(&[0.0, 0.0, 0.0, 0.0]); let mut mdu = Vector::new(neq); diff --git a/russell_sparse/tests/test_umfpack.rs b/russell_sparse/tests/test_umfpack.rs index 76761946..e318f1b7 100644 --- a/russell_sparse/tests/test_umfpack.rs +++ b/russell_sparse/tests/test_umfpack.rs @@ -6,7 +6,7 @@ use russell_sparse::StrError; fn test_complex_umfpack() -> Result<(), StrError> { let n = 10; let d = (n as f64) / 10.0; - let mut coo = CooMatrix::new(n, n, n, None)?; + let mut coo = CooMatrix::new(n, n, n, Sym::No)?; let mut x_correct = Vector::new(n); let mut rhs = Vector::new(n); for k in 0..n {