[sparse] [Important] Remove the one_based flag. Store Fortran indices…

… in MUMPS struct directly

README.md

@@ -298,7 +298,7 @@ fn main() -> Result<(), StrError> {
     //  . -1 -3  2  .
     //  .  .  1  .  .
     //  . 4  2  .  1
-    let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None, false)?;
+    let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?;
     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)?;

russell_ode/src/euler_backward.rs

@@ -1,7 +1,7 @@
 use crate::StrError;
 use crate::{OdeSolverTrait, Params, System, Workspace};
 use russell_lab::{vec_copy, vec_rms_scaled, vec_update, Vector};
-use russell_sparse::{CooMatrix, Genie, LinSolver, SparseMatrix};
+use russell_sparse::{CooMatrix, LinSolver, SparseMatrix};
 
 /// Implements the backward Euler (implicit) solver
 pub(crate) struct EulerBackward<'a, F, J, A>
@@ -51,15 +51,14 @@ where
         };
         let nnz = jac_nnz + ndim; // +ndim corresponds to the diagonal I matrix
         let symmetry = Some(system.jac_symmetry);
-        let one_based = params.newton.genie == Genie::Mumps;
         EulerBackward {
             params,
             system,
             k: Vector::new(ndim),
             w: Vector::new(ndim),
             r: Vector::new(ndim),
             dy: Vector::new(ndim),
-            kk: SparseMatrix::new_coo(ndim, ndim, nnz, symmetry, one_based).unwrap(),
+            kk: SparseMatrix::new_coo(ndim, ndim, nnz, symmetry).unwrap(),
             solver: LinSolver::new(params.newton.genie).unwrap(),
         }
     }

russell_ode/src/radau5.rs

@@ -111,7 +111,6 @@ where
     pub fn new(params: Params, system: &'a System<F, J, A>) -> Self {
         let ndim = system.ndim;
         let symmetry = Some(system.jac_symmetry);
-        let one_based = params.newton.genie == Genie::Mumps;
         let mass_nnz = match system.mass_matrix.as_ref() {
             Some(mass) => mass.get_info().2,
             None => ndim,
@@ -126,9 +125,9 @@ where
         Radau5 {
             params,
             system,
-            jj: SparseMatrix::new_coo(ndim, ndim, jac_nnz, symmetry, one_based).unwrap(),
-            kk_real: SparseMatrix::new_coo(ndim, ndim, nnz, symmetry, one_based).unwrap(),
-            kk_comp: ComplexSparseMatrix::new_coo(ndim, ndim, nnz, symmetry, one_based).unwrap(),
+            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(),
             solver_real: LinSolver::new(params.newton.genie).unwrap(),
             solver_comp: ComplexLinSolver::new(params.newton.genie).unwrap(),
             reuse_jacobian: false,

russell_ode/src/samples.rs

@@ -200,7 +200,7 @@ impl Samples {
 
         // mass matrix
         let mass_nnz = 5;
-        system.init_mass_matrix(mass_nnz, false).unwrap();
+        system.init_mass_matrix(mass_nnz).unwrap();
         system.mass_put(0, 0, 1.0).unwrap();
         system.mass_put(0, 1, 1.0).unwrap();
         system.mass_put(1, 0, 1.0).unwrap();
@@ -751,7 +751,7 @@ impl Samples {
         const C2: f64 = 2e-6;
         const C3: f64 = 3e-6;
         let mass_nnz = 9;
-        system.init_mass_matrix(mass_nnz, false).unwrap();
+        system.init_mass_matrix(mass_nnz).unwrap();
         system.mass_put(0, 0, -C1).unwrap();
         system.mass_put(0, 1, C1).unwrap();
         system.mass_put(1, 0, C1).unwrap();
@@ -836,7 +836,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, false).unwrap();
+        let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
         (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();
 
         // compute the numerical Jacobian matrix
@@ -864,7 +864,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, false).unwrap();
+        let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
         (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();
 
         // compute the numerical Jacobian matrix
@@ -892,7 +892,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, false).unwrap();
+        let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
         (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();
 
         // compute the numerical Jacobian matrix
@@ -920,7 +920,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, false).unwrap();
+        let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
         (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();
 
         // compute the numerical Jacobian matrix
@@ -940,7 +940,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, false).unwrap();
+        let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
         (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();
 
         // compute the numerical Jacobian matrix
@@ -960,7 +960,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, false).unwrap();
+        let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
         (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();
 
         // compute the numerical Jacobian matrix
@@ -980,7 +980,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, false).unwrap();
+        let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
         (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();
 
         // compute the numerical Jacobian matrix
@@ -1000,7 +1000,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, false).unwrap();
+        let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
         (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();
 
         // compute the numerical Jacobian matrix
@@ -1020,7 +1020,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, false).unwrap();
+        let mut jj = CooMatrix::new(system.ndim, system.ndim, system.jac_nnz, symmetry).unwrap();
         (system.jacobian)(&mut jj, data.x0, &data.y0, multiplier, &mut args).unwrap();
 
         // compute the numerical Jacobian matrix

russell_ode/src/system.rs

@@ -131,15 +131,13 @@ where
     /// # Input
     ///
     /// * `max_nnz` -- Max number of non-zero values
-    /// * `one_based` -- Flag indicating whether the Sparse matrix can be used with Fortran code (e.g., MUMPS) or not.
-    ///   Make sure that this flag is the same used for the Jacobian matrix.
-    pub fn init_mass_matrix(&mut self, max_nnz: usize, one_based: bool) -> Result<(), StrError> {
+    pub fn init_mass_matrix(&mut self, max_nnz: usize) -> Result<(), StrError> {
         let sym = if self.jac_symmetry == Symmetry::No {
             None
         } else {
             Some(self.jac_symmetry)
         };
-        self.mass_matrix = Some(CooMatrix::new(self.ndim, self.ndim, max_nnz, sym, one_based).unwrap());
+        self.mass_matrix = Some(CooMatrix::new(self.ndim, self.ndim, max_nnz, sym).unwrap());
         Ok(())
     }
 
@@ -253,7 +251,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, false).unwrap();
+        let mut jj = CooMatrix::new(2, 2, 2, None).unwrap();
         let m = 1.0;
         assert_eq!(
             (ode.jacobian)(&mut jj, x, &y, m, &mut args),
@@ -303,7 +301,7 @@ mod tests {
         // analytical_solution:
         // y[0] = f64::cos(x * x / 2.0) - 2.0 * f64::sin(x * x / 2.0);
         // y[1] = 2.0 * f64::cos(x * x / 2.0) + f64::sin(x * x / 2.0);
-        ode.init_mass_matrix(2, false).unwrap(); // diagonal mass matrix => OK, but not needed
+        ode.init_mass_matrix(2).unwrap(); // diagonal mass matrix => OK, but not needed
         ode.mass_put(0, 0, 1.0).unwrap();
         ode.mass_put(1, 1, 1.0).unwrap();
         // call system function
@@ -312,7 +310,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, false).unwrap();
+        let mut jj = CooMatrix::new(2, 2, 2, None).unwrap();
         let m = 1.0;
         (ode.jacobian)(&mut jj, x, &y, m, &mut args).unwrap();
         // check

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, false)?;
+    let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?;
     coo.put(0, 0, 0.2)?;
     coo.put(0, 1, 0.2)?;
     coo.put(1, 0, 0.5)?;

russell_sparse/examples/complex_system.rs

@@ -43,8 +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 one_based = genie == Genie::Mumps;
-    let mut coo = ComplexSparseMatrix::new_coo(ndim, ndim, nnz, None, one_based)?;
+    let mut coo = ComplexSparseMatrix::new_coo(ndim, ndim, nnz, None)?;
 
     // first column
     coo.put(0, 0, cpx!(19.73, 0.00))?;

russell_sparse/examples/doc_coo_from_arrays.rs

@@ -16,7 +16,7 @@ fn main() -> Result<(), StrError> {
         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, false)?;
+    let coo = CooMatrix::from(nrow, ncol, row_indices, col_indices, values, symmetry)?;
 
     // covert to dense
     let a = coo.as_dense();

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, true)?;
+    let mut coo = CooMatrix::new(5, 5, 13, None)?;
     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)?;

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, false)?;
+    let mut coo = CooMatrix::new(nrow, ncol, nnz, None)?;
     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)?;

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, false)?;
+    let mut coo = CooMatrix::new(nrow, ncol, nnz, None)?;
     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)?;

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, false)?;
+    let mut mat = SparseMatrix::new_coo(ndim, ndim, nnz, None)?;
     mat.put(0, 0, 0.2)?;
     mat.put(0, 1, 0.2)?;
     mat.put(1, 0, 0.5)?;

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, false)?;
+    let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?;
     coo.put(0, 0, 0.2)?;
     coo.put(0, 1, 0.2)?;
     coo.put(1, 0, 0.5)?;

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, false)?;
+    let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?;
     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)?;

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, false)?;
+    let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?;
     coo.put(0, 0, 0.2)?;
     coo.put(0, 1, 0.2)?;
     coo.put(1, 0, 0.5)?;

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, true)?;
+    let mut coo = SparseMatrix::new_coo(ndim, ndim, nnz, None)?;
     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)?;

russell_sparse/examples/nonlinear_system_4eqs.rs

@@ -29,9 +29,8 @@ fn main() -> Result<(), StrError> {
     let mut solver = LinSolver::new(genie)?;
 
     // allocate Jacobian matrix (J) as SparseMatrix
-    let one_based = if genie == Genie::Mumps { true } else { false };
     let (neq, nnz) = (4, 16);
-    let mut jj = SparseMatrix::new_coo(neq, neq, nnz, None, one_based).unwrap();
+    let mut jj = SparseMatrix::new_coo(neq, neq, nnz, None).unwrap();
 
     // allocate residual (rr), vector of unknowns (uu), and minus-uu (mdu)
     let mut rr = Vector::new(neq);

russell_sparse/src/bin/mem_check.rs

@@ -20,8 +20,7 @@ fn test_solver(genie: Genie) {
         }
     };
 
-    let one_based = genie == Genie::Mumps;
-    let (coo, _, _, _) = Samples::umfpack_unsymmetric_5x5(one_based);
+    let (coo, _, _, _) = Samples::umfpack_unsymmetric_5x5();
     let mut mat = SparseMatrix::from_coo(coo);
 
     match solver.actual.factorize(&mut mat, None) {
@@ -73,7 +72,7 @@ fn test_complex_solver(genie: Genie) {
     };
 
     let coo = match genie {
-        Genie::Mumps => Samples::complex_symmetric_3x3_lower(true).0,
+        Genie::Mumps => Samples::complex_symmetric_3x3_lower().0,
         Genie::Umfpack => Samples::complex_symmetric_3x3_full().0,
         Genie::IntelDss => panic!("TODO"),
     };
@@ -129,8 +128,7 @@ fn test_solver_singular(genie: Genie) {
         }
     };
 
-    let one_based = if genie == Genie::Mumps { true } else { false };
-    let mut coo_singular = match SparseMatrix::new_coo(ndim, ndim, nnz, None, one_based) {
+    let mut coo_singular = match SparseMatrix::new_coo(ndim, ndim, nnz, None) {
         Ok(v) => v,
         Err(e) => {
             println!("FAIL(new CooMatrix): {}", e);

russell_sparse/src/bin/solve_matrix_market.rs

@@ -70,10 +70,10 @@ fn main() -> Result<(), StrError> {
     let genie = Genie::from(&opt.genie);
 
     // select the symmetric handling option
-    let (handling, one_based) = match genie {
-        Genie::Mumps => (MMsymOption::LeaveAsLower, true),
-        Genie::Umfpack => (MMsymOption::MakeItFull, false),
-        Genie::IntelDss => (MMsymOption::SwapToUpper, false),
+    let handling = match genie {
+        Genie::Mumps => MMsymOption::LeaveAsLower,
+        Genie::Umfpack => MMsymOption::MakeItFull,
+        Genie::IntelDss => MMsymOption::SwapToUpper,
     };
 
     // configuration parameters
@@ -99,7 +99,7 @@ fn main() -> Result<(), StrError> {
 
     // read the matrix
     let mut sw = Stopwatch::new();
-    let coo = read_matrix_market(&opt.matrix_market_file, handling, one_based)?;
+    let coo = read_matrix_market(&opt.matrix_market_file, handling)?;
     stats.time_nanoseconds.read_matrix = sw.stop();
 
     // save the COO matrix as a generic SparseMatrix