Implement coo and csr sparse matrices

.vscode/launch.json

@@ -159,15 +159,15 @@
     {
       "type": "lldb",
       "request": "launch",
-      "name": "Debug executable 'solve_mm_build'",
+      "name": "Debug executable 'solve_matrix_market_build'",
       "cargo": {
         "args": [
           "build",
-          "--bin=solve_mm_build",
+          "--bin=solve_matrix_market_build",
           "--package=russell_sparse"
         ],
         "filter": {
-          "name": "solve_mm_build",
+          "name": "solve_matrix_market_build",
           "kind": "bin"
         }
       },
@@ -177,16 +177,16 @@
     {
       "type": "lldb",
       "request": "launch",
-      "name": "Debug unit tests in executable 'solve_mm_build'",
+      "name": "Debug unit tests in executable 'solve_matrix_market_build'",
       "cargo": {
         "args": [
           "test",
           "--no-run",
-          "--bin=solve_mm_build",
+          "--bin=solve_matrix_market_build",
           "--package=russell_sparse"
         ],
         "filter": {
-          "name": "solve_mm_build",
+          "name": "solve_matrix_market_build",
           "kind": "bin"
         }
       },

.vscode/settings.json

@@ -30,6 +30,7 @@
     "lsol",
     "nelt",
     "odyad",
+    "PERMUT",
     "Schur",
     "Teukolsky",
     "tgamma",

README.md

@@ -121,32 +121,31 @@ fn main() -> Result<(), StrError> {
 ### Solve a sparse linear system
 
 ```rust
-use russell_lab::{Matrix, Vector, StrError};
-use russell_sparse::{ConfigSolver, Solver, SparseTriplet};
+use russell_lab::{Matrix, Vector};
+use russell_sparse::prelude::*;
+use russell_sparse::StrError;
 
 fn main() -> Result<(), StrError> {
-
     // allocate a square matrix
-    let neq = 5; // number of equations
-    let nnz = 13; // number of non-zeros
-    let mut trip = SparseTriplet::new(neq, nnz)?;
-    trip.put(0, 0,  1.0)?; // << (0, 0, a00/2)
-    trip.put(0, 0,  1.0)?; // << (0, 0, a00/2)
-    trip.put(1, 0,  3.0)?;
-    trip.put(0, 1,  3.0)?;
-    trip.put(2, 1, -1.0)?;
-    trip.put(4, 1,  4.0)?;
-    trip.put(1, 2,  4.0)?;
-    trip.put(2, 2, -3.0)?;
-    trip.put(3, 2,  1.0)?;
-    trip.put(4, 2,  2.0)?;
-    trip.put(2, 3,  2.0)?;
-    trip.put(1, 4,  6.0)?;
-    trip.put(4, 4,  1.0)?;
+    let (nrow, ncol, nnz) = (5, 5, 13);
+    let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?;
+    coo.put(0, 0, 1.0)?; // << (0, 0, a00/2)
+    coo.put(0, 0, 1.0)?; // << (0, 0, a00/2)
+    coo.put(1, 0, 3.0)?;
+    coo.put(0, 1, 3.0)?;
+    coo.put(2, 1, -1.0)?;
+    coo.put(4, 1, 4.0)?;
+    coo.put(1, 2, 4.0)?;
+    coo.put(2, 2, -3.0)?;
+    coo.put(3, 2, 1.0)?;
+    coo.put(4, 2, 2.0)?;
+    coo.put(2, 3, 2.0)?;
+    coo.put(1, 4, 6.0)?;
+    coo.put(4, 4, 1.0)?;
 
     // print matrix
-    let mut a = Matrix::new(neq, neq);
-    trip.to_matrix(&mut a)?;
+    let mut a = Matrix::new(nrow, ncol);
+    coo.to_matrix(&mut a)?;
     let correct = "┌                ┐\n\
                    │  2  3  0  0  0 │\n\
                    │  3  0  4  0  6 │\n\
@@ -157,13 +156,13 @@ fn main() -> Result<(), StrError> {
     assert_eq!(format!("{}", a), correct);
 
     // allocate x and rhs
-    let mut x = Vector::new(neq);
+    let mut x = Vector::new(nrow);
     let rhs = Vector::from(&[8.0, 45.0, -3.0, 3.0, 19.0]);
 
     // initialize, factorize, and solve
     let config = ConfigSolver::new();
-    let mut solver = Solver::new(config, neq, nnz, None)?;
-    solver.factorize(&trip)?;
+    let mut solver = Solver::new(config, nrow, nnz, None)?;
+    solver.factorize(&coo)?;
     solver.solve(&mut x, &rhs)?;
     let correct = "┌          ┐\n\
                    │ 1.000000 │\n\

russell_sparse/Cargo.toml

@@ -16,6 +16,8 @@ russell_chk = { path = "../russell_chk", version = "0.5.0" }
 russell_lab = { path = "../russell_lab", version = "0.5.0" }
 russell_openblas = { path = "../russell_openblas", version = "0.5.0" }
 structopt = "0.3"
+serde = { version = "1.0", features = ["derive"] }
+serde_json = "1.0"
 
 [build-dependencies]
 cc = "1.0"

russell_sparse/README.md

@@ -65,45 +65,47 @@ export OPENBLAS_NUM_THREADS=1
 
 ```rust
 use russell_lab::{Matrix, Vector};
-use russell_sparse::{ConfigSolver, Solver, SparseTriplet, StrError};
+use russell_sparse::prelude::*;
+use russell_sparse::StrError;
 
 fn main() -> Result<(), StrError> {
     // allocate a square matrix
-    let neq = 3; // number of equations
+    let nrow = 3; // number of equations
+    let ncol = nrow; // number of equations
     let nnz = 5; // number of non-zeros
-    let mut trip = SparseTriplet::new(neq, nnz)?;
-    trip.put(0, 0, 0.2)?;
-    trip.put(0, 1, 0.2)?;
-    trip.put(1, 0, 0.5)?;
-    trip.put(1, 1, -0.25)?;
-    trip.put(2, 2, 0.25)?;
-    
+    let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?;
+    coo.put(0, 0, 0.2)?;
+    coo.put(0, 1, 0.2)?;
+    coo.put(1, 0, 0.5)?;
+    coo.put(1, 1, -0.25)?;
+    coo.put(2, 2, 0.25)?;
+
     // print matrix
-    let mut a = Matrix::new(neq, neq);
-    trip.to_matrix(&mut a)?;
+    let mut a = Matrix::new(nrow, ncol);
+    coo.to_matrix(&mut a)?;
     let correct = "┌                   ┐\n\
                    │   0.2   0.2     0 │\n\
                    │   0.5 -0.25     0 │\n\
                    │     0     0  0.25 │\n\
                    └                   ┘";
     assert_eq!(format!("{}", a), correct);
-    
+
     // allocate rhs
     let rhs1 = Vector::from(&[1.0, 1.0, 1.0]);
     let rhs2 = Vector::from(&[2.0, 2.0, 2.0]);
-    
+
     // calculate solution
     let config = ConfigSolver::new();
-    let (mut solver, x1) = Solver::compute(config, &trip, &rhs1)?;
+    let (mut solver, x1) = Solver::compute(config, &coo, &rhs1)?;
     let correct1 = "┌   ┐\n\
                     │ 3 │\n\
                     │ 2 │\n\
                     │ 4 │\n\
                     └   ┘";
     assert_eq!(format!("{}", x1), correct1);
-    
+
     // solve again
-    let mut x2 = Vector::new(neq);
+    let mut x2 = Vector::new(nrow);
     solver.solve(&mut x2, &rhs2)?;
     let correct2 = "┌   ┐\n\
                     │ 6 │\n\
@@ -121,9 +123,9 @@ We wrap two direct sparse solvers: UMFPACK (aka **UMF**) and MUMPS (aka **MMP**)
 
 ## Tools
 
-This crate includes a tool named `solve_mm_build` to study the performance of the available sparse solvers (currently MMP and UMF). The `_build` suffix is to disable the coverage tool.
+This crate includes a tool named `solve_matrix_market_build` to study the performance of the available sparse solvers (currently MMP and UMF). The `_build` suffix is to disable the coverage tool.
 
-`solve_mm_build` reads a [Matrix Market file](https://math.nist.gov/MatrixMarket/formats.html) and solves the linear system:
+`solve_matrix_market_build` reads a [Matrix Market file](https://math.nist.gov/MatrixMarket/formats.html) and solves the linear system:
 
 ```text
 a ⋅ x = rhs
@@ -136,13 +138,13 @@ The data directory contains an example of Matrix Market file named `bfwb62.mtx`
 Run the command:
 
 ```bash
-cargo run --release --bin solve_mm_build -- ~/Downloads/matrix-market/bfwb62.mtx
+cargo run --release --bin solve_matrix_market_build -- ~/Downloads/matrix-market/bfwb62.mtx
 ```
 
 Or
 
 ```bash
-cargo run --release --bin solve_mm_build -- --help
+cargo run --release --bin solve_matrix_market_build -- --help
 ```
 
 for more options.

russell_sparse/c_code/constants.h

@@ -20,74 +20,74 @@ const MUMPS_INT MUMPS_JOB_ANALYZE = 1;
 const MUMPS_INT MUMPS_JOB_FACTORIZE = 2;
 const MUMPS_INT MUMPS_JOB_SOLVE = 3;
 
-const MUMPS_INT MUMPS_PAR_HOST_ALSO_WORKS = 1;      // section 5.1.4, page 26
-const MUMPS_INT MUMPS_ICNTL5_ASSEMBLED_MATRIX = 0;  // section 5.2.2, page 27
-const MUMPS_INT MUMPS_ICNTL18_CENTRALIZED = 0;      // section 5.2.2, page 27
-const MUMPS_INT MUMPS_ICNTL6_PERMUT_AUTO = 7;       // section 5.3, page 32
-const MUMPS_INT MUMPS_ICNTL28_SEQUENTIAL = 1;       // section 5.4, page 33
+const MUMPS_INT MUMPS_PAR_HOST_ALSO_WORKS = 1;     // section 5.1.4, page 26
+const MUMPS_INT MUMPS_ICNTL5_ASSEMBLED_MATRIX = 0; // section 5.2.2, page 27
+const MUMPS_INT MUMPS_ICNTL18_CENTRALIZED = 0;     // section 5.2.2, page 27
+const MUMPS_INT MUMPS_ICNTL6_PERMUT_AUTO = 7;      // section 5.3, page 32
+const MUMPS_INT MUMPS_ICNTL28_SEQUENTIAL = 1;      // section 5.4, page 33
 
-const double UMF_PRINT_LEVEL_SILENT = 0.0;   // page 116
-const double UMF_PRINT_LEVEL_VERBOSE = 2.0;  // page 116
+const double UMF_PRINT_LEVEL_SILENT = 0.0;  // page 116
+const double UMF_PRINT_LEVEL_VERBOSE = 2.0; // page 116
 
 const MUMPS_INT MMP_SYMMETRY[3] = {
-    0,  // Unsymmetric
-    1,  // Positive-definite symmetric
-    2,  // General symmetric
+    0, // Unsymmetric
+    1, // Positive-definite symmetric
+    2, // General symmetric
 };
 
 const MUMPS_INT MMP_ORDERING[10] = {
-    0,  // Amd
-    2,  // Amf
-    7,  // Auto
-    7,  // Best => Auto
-    7,  // Cholmod => Auto
-    5,  // Metis
-    7,  // No => Auto
-    4,  // Pord
-    6,  // Qamd
-    3,  // Scotch
+    0, // 0: Amd
+    2, // 1: Amf
+    7, // 2: Auto
+    7, // 3: Best => Auto
+    7, // 4: Cholmod => Auto
+    5, // 5: Metis
+    7, // 6: No => Auto
+    4, // 7: Pord
+    6, // 8: Qamd
+    3, // 9: Scotch
 };
 
 const MUMPS_INT MMP_SCALING[9] = {
-    77,  // Auto
-    3,   // Column
-    1,   // Diagonal
-    77,  // Max => Auto
-    0,   // No
-    4,   // RowCol
-    7,   // RowColIter
-    8,   // RowColRig
-    77,  // Sum => Auto
+    77, // 0: Auto
+    3,  // 1: Column
+    1,  // 2: Diagonal
+    77, // 3: Max => Auto
+    0,  // 4: No
+    4,  // 5: RowCol
+    7,  // 6: RowColIter
+    8,  // 7: RowColRig
+    77, // 8: Sum => Auto
 };
 
 const double UMF_SYMMETRY[2] = {
-    UMFPACK_STRATEGY_UNSYMMETRIC,  // Unsymmetric
-    UMFPACK_STRATEGY_SYMMETRIC,    // General symmetric
+    UMFPACK_STRATEGY_UNSYMMETRIC, // Unsymmetric
+    UMFPACK_STRATEGY_SYMMETRIC,   // General symmetric
 };
 
 const double UMF_ORDERING[10] = {
-    UMFPACK_ORDERING_AMD,      // Amd
-    UMFPACK_DEFAULT_ORDERING,  // Amf => Auto
-    UMFPACK_DEFAULT_ORDERING,  // Auto
-    UMFPACK_ORDERING_BEST,     // Best
-    UMFPACK_ORDERING_CHOLMOD,  // Cholmod
-    UMFPACK_ORDERING_METIS,    // Metis
-    UMFPACK_ORDERING_NONE,     // No
-    UMFPACK_DEFAULT_ORDERING,  // Pord => Auto
-    UMFPACK_DEFAULT_ORDERING,  // Qamd => Auto
-    UMFPACK_DEFAULT_ORDERING,  // Scotch => Auto
+    UMFPACK_ORDERING_AMD,     // Amd
+    UMFPACK_DEFAULT_ORDERING, // Amf => Auto
+    UMFPACK_DEFAULT_ORDERING, // Auto
+    UMFPACK_ORDERING_BEST,    // Best
+    UMFPACK_ORDERING_CHOLMOD, // Cholmod
+    UMFPACK_ORDERING_METIS,   // Metis
+    UMFPACK_ORDERING_NONE,    // No
+    UMFPACK_DEFAULT_ORDERING, // Pord => Auto
+    UMFPACK_DEFAULT_ORDERING, // Qamd => Auto
+    UMFPACK_DEFAULT_ORDERING, // Scotch => Auto
 };
 
 const double UMF_SCALING[9] = {
-    UMFPACK_DEFAULT_SCALE,  // Auto
-    UMFPACK_DEFAULT_SCALE,  // Column => Auto
-    UMFPACK_DEFAULT_SCALE,  // Diagonal => Auto
-    UMFPACK_SCALE_MAX,      // Max
-    UMFPACK_SCALE_NONE,     // No
-    UMFPACK_DEFAULT_SCALE,  // RowCol => Auto
-    UMFPACK_DEFAULT_SCALE,  // RowColIter => Auto
-    UMFPACK_DEFAULT_SCALE,  // RowColRig => Auto
-    UMFPACK_SCALE_SUM,      // Sum
+    UMFPACK_DEFAULT_SCALE, // Auto
+    UMFPACK_DEFAULT_SCALE, // Column => Auto
+    UMFPACK_DEFAULT_SCALE, // Diagonal => Auto
+    UMFPACK_SCALE_MAX,     // Max
+    UMFPACK_SCALE_NONE,    // No
+    UMFPACK_DEFAULT_SCALE, // RowCol => Auto
+    UMFPACK_DEFAULT_SCALE, // RowColIter => Auto
+    UMFPACK_DEFAULT_SCALE, // RowColRig => Auto
+    UMFPACK_SCALE_SUM,     // Sum
 };
 
 #endif

russell_sparse/c_code/solver_mmp.h

@@ -9,27 +9,27 @@
 #include "constants.h"
 #include "dmumps_c.h"
 
-#define ICNTL(i) icntl[(i)-1]  // macro to make indices match documentation
-#define INFOG(i) infog[(i)-1]  // macro to make indices match documentation
-#define INFO(i) info[(i)-1]    // macro to make indices match documentation
+#define ICNTL(i) icntl[(i)-1] // macro to make indices match documentation
+#define INFOG(i) infog[(i)-1] // macro to make indices match documentation
+#define INFO(i) info[(i)-1]   // macro to make indices match documentation
 
 static inline void set_mmp_verbose(DMUMPS_STRUC_C *data, int32_t verbose) {
     if (verbose == C_TRUE) {
-        data->ICNTL(1) = 6;  // standard output stream
-        data->ICNTL(2) = 0;  // output stream
-        data->ICNTL(3) = 6;  // standard output stream
-        data->ICNTL(4) = 3;  // errors, warnings, and main statistics printed
+        data->ICNTL(1) = 6; // standard output stream
+        data->ICNTL(2) = 0; // output stream
+        data->ICNTL(3) = 6; // standard output stream
+        data->ICNTL(4) = 3; // errors, warnings, and main statistics printed
     } else {
-        data->ICNTL(1) = -1;  // no output messages
-        data->ICNTL(2) = -1;  // no warnings
-        data->ICNTL(3) = -1;  // no global information
-        data->ICNTL(4) = -1;  // message level
+        data->ICNTL(1) = -1; // no output messages
+        data->ICNTL(2) = -1; // no warnings
+        data->ICNTL(3) = -1; // no global information
+        data->ICNTL(4) = -1; // message level
     }
 }
 
 struct SolverMMP {
-    DMUMPS_STRUC_C data;    // data structure
-    int32_t done_job_init;  // job init successfully
+    DMUMPS_STRUC_C data;   // data structure
+    int32_t done_job_init; // job init successfully
 };
 
 struct SolverMMP *new_solver_mmp() {
@@ -141,7 +141,7 @@ int32_t solver_mmp_initialize(struct SolverMMP *solver,
     solver->data.ICNTL(28) = MUMPS_ICNTL28_SEQUENTIAL;
     solver->data.ICNTL(29) = MUMPS_IGNORED;
 
-    return 0;  // success
+    return 0; // success
 }
 
 int32_t solver_mmp_factorize(struct SolverMMP *solver,