diff --git a/.vscode/launch.json b/.vscode/launch.json index fb7b4828..4c36f330 100644 --- a/.vscode/launch.json +++ b/.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" } }, diff --git a/.vscode/settings.json b/.vscode/settings.json index aa4c8e94..d297b2da 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -30,6 +30,7 @@ "lsol", "nelt", "odyad", + "PERMUT", "Schur", "Teukolsky", "tgamma", diff --git a/README.md b/README.md index 5c80b892..6ace29a1 100644 --- a/README.md +++ b/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\ diff --git a/russell_sparse/Cargo.toml b/russell_sparse/Cargo.toml index e387dc6e..775e8d0d 100644 --- a/russell_sparse/Cargo.toml +++ b/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" diff --git a/russell_sparse/README.md b/russell_sparse/README.md index 1ad11e8d..976896b5 100644 --- a/russell_sparse/README.md +++ b/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. diff --git a/russell_sparse/c_code/constants.h b/russell_sparse/c_code/constants.h index 2b43b51a..bd36665c 100644 --- a/russell_sparse/c_code/constants.h +++ b/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 diff --git a/russell_sparse/c_code/solver_mmp.h b/russell_sparse/c_code/solver_mmp.h index 3c661654..a5722233 100644 --- a/russell_sparse/c_code/solver_mmp.h +++ b/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, diff --git a/russell_sparse/data/figures/doc-example-vismatrix.png b/russell_sparse/data/figures/doc-example-vismatrix.png new file mode 100644 index 00000000..6359ff28 Binary files /dev/null and b/russell_sparse/data/figures/doc-example-vismatrix.png differ diff --git a/russell_sparse/data/matrix_market/bad_symmetric_rectangular.mtx b/russell_sparse/data/matrix_market/bad_symmetric_rectangular.mtx new file mode 100644 index 00000000..77bb89a7 --- /dev/null +++ b/russell_sparse/data/matrix_market/bad_symmetric_rectangular.mtx @@ -0,0 +1,6 @@ +%%MatrixMarket matrix coordinate real symmetric +4 3 4 + 1 1 1.0 + 2 1 2.0 + 2 2 3.0 + 3 2 4.0 \ No newline at end of file diff --git a/russell_sparse/data/matrix_market/ok1.mtx b/russell_sparse/data/matrix_market/ok_general.mtx similarity index 100% rename from russell_sparse/data/matrix_market/ok1.mtx rename to russell_sparse/data/matrix_market/ok_general.mtx diff --git a/russell_sparse/data/matrix_market/bad_rectangular.mtx b/russell_sparse/data/matrix_market/ok_rectangular.mtx similarity index 59% rename from russell_sparse/data/matrix_market/bad_rectangular.mtx rename to russell_sparse/data/matrix_market/ok_rectangular.mtx index cb1babf1..7728b329 100644 --- a/russell_sparse/data/matrix_market/bad_rectangular.mtx +++ b/russell_sparse/data/matrix_market/ok_rectangular.mtx @@ -1,3 +1,7 @@ %%MatrixMarket matrix coordinate real general % comment line -1 2 3 \ No newline at end of file +2 3 4 +0 0 10.0 +0 1 20.0 +1 1 30.0 +1 2 40.0 \ No newline at end of file diff --git a/russell_sparse/data/matrix_market/simple_gen.mtx b/russell_sparse/data/matrix_market/ok_simple_general.mtx similarity index 100% rename from russell_sparse/data/matrix_market/simple_gen.mtx rename to russell_sparse/data/matrix_market/ok_simple_general.mtx diff --git a/russell_sparse/data/matrix_market/simple_sym.mtx b/russell_sparse/data/matrix_market/ok_simple_symmetric.mtx similarity index 100% rename from russell_sparse/data/matrix_market/simple_sym.mtx rename to russell_sparse/data/matrix_market/ok_simple_symmetric.mtx diff --git a/russell_sparse/data/matrix_market/ok2.mtx b/russell_sparse/data/matrix_market/ok_symmetric.mtx similarity index 70% rename from russell_sparse/data/matrix_market/ok2.mtx rename to russell_sparse/data/matrix_market/ok_symmetric.mtx index 023923d2..c724db9c 100644 --- a/russell_sparse/data/matrix_market/ok2.mtx +++ b/russell_sparse/data/matrix_market/ok_symmetric.mtx @@ -10,13 +10,13 @@ 3 3 9.0 4 4 7.0 5 5 8.0 -1 2 1.0 -1 3 1.0 -1 4 3.0 -1 5 2.0 -2 3 2.0 -2 4 1.0 -2 5 1.0 -3 4 1.0 -3 5 5.0 -4 5 1.0 +2 1 1.0 +3 1 1.0 +4 1 3.0 +5 1 2.0 +3 2 2.0 +4 2 1.0 +5 2 1.0 +4 3 1.0 +5 3 5.0 +5 4 1.0 diff --git a/russell_sparse/data/matrix_market/ok3.mtx b/russell_sparse/data/matrix_market/ok_symmetric_small.mtx similarity index 100% rename from russell_sparse/data/matrix_market/ok3.mtx rename to russell_sparse/data/matrix_market/ok_symmetric_small.mtx diff --git a/russell_sparse/data/results/solve-matrix-market-mumps-bfwb62.json b/russell_sparse/data/results/solve-matrix-market-mumps-bfwb62.json new file mode 100644 index 00000000..cd150ddd --- /dev/null +++ b/russell_sparse/data/results/solve-matrix-market-mumps-bfwb62.json @@ -0,0 +1,30 @@ +{ + "platform": "Russell", + "blas_lib": "OpenBLAS", + "solver_name": "MUMPS-local", + "matrix_name": "bfwb62", + "symmetry": "General", + "layout": "Lower", + "nrow": 62, + "ncol": 62, + "nnz": 202, + "time_read_matrix_market_nanosecond": 33245, + "time_read_matrix_market_human": "33.245µs", + "time_factorize_nanosecond": 338231, + "time_factorize_human": "338.231µs", + "time_solve_nanosecond": 114215, + "time_solve_human": "114.215µs", + "time_factorize_and_solve_nanosecond": 452446, + "time_factorize_and_solve_human": "452.446µs", + "requested_ordering": "Auto", + "requested_scaling": "Auto", + "requested_openmp_num_threads": 1, + "effective_ordering": "Amf", + "effective_scaling": "RowColIter", + "effective_openmp_num_threads": 32, + "verify_max_abs_a": 0.0001, + "verify_max_abs_a_times_x": 1.0000000000000004, + "verify_relative_error": 5.550560067119071e-16, + "verify_time_nanosecond": 1838, + "verify_time_human": "1.838µs" +} diff --git a/russell_sparse/data/results/solve-matrix-market-umfpack-bfwb62.json b/russell_sparse/data/results/solve-matrix-market-umfpack-bfwb62.json new file mode 100644 index 00000000..458ef4f6 --- /dev/null +++ b/russell_sparse/data/results/solve-matrix-market-umfpack-bfwb62.json @@ -0,0 +1,30 @@ +{ + "platform": "Russell", + "blas_lib": "OpenBLAS", + "solver_name": "UMFPACK", + "matrix_name": "bfwb62", + "symmetry": "General", + "layout": "Full", + "nrow": 62, + "ncol": 62, + "nnz": 342, + "time_read_matrix_market_nanosecond": 37481, + "time_read_matrix_market_human": "37.481µs", + "time_factorize_nanosecond": 176701, + "time_factorize_human": "176.701µs", + "time_solve_nanosecond": 6181, + "time_solve_human": "6.181µs", + "time_factorize_and_solve_nanosecond": 182882, + "time_factorize_and_solve_human": "182.882µs", + "requested_ordering": "Auto", + "requested_scaling": "Auto", + "requested_openmp_num_threads": 1, + "effective_ordering": "Amd", + "effective_scaling": "Sum", + "effective_openmp_num_threads": 32, + "verify_max_abs_a": 0.0001, + "verify_max_abs_a_times_x": 1.0000000000000004, + "verify_relative_error": 4.4404480536952566e-16, + "verify_time_nanosecond": 1364, + "verify_time_human": "1.364µs" +} diff --git a/russell_sparse/examples/doc_csr_from_small_coo.rs b/russell_sparse/examples/doc_csr_from_small_coo.rs new file mode 100644 index 00000000..04e004d3 --- /dev/null +++ b/russell_sparse/examples/doc_csr_from_small_coo.rs @@ -0,0 +1,56 @@ +use russell_sparse::prelude::*; +use russell_sparse::StrError; + +fn main() -> Result<(), StrError> { + // allocate a square matrix and store as COO matrix + // ┌ ┐ + // │ 2 3 0 0 0 │ + // │ 3 0 4 0 6 │ + // │ 0 -1 -3 2 0 │ + // │ 0 0 1 0 0 │ + // │ 0 4 2 0 1 │ + // └ ┘ + 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) duplicate + coo.put(0, 0, 1.0)?; // << (0, 0, a00/2) duplicate + 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)?; + + // convert to CSR matrix + let csr = CsrMatrix::from(&coo); + let correct_j = &[ + // p + 0, 1, // i = 0, count = 0, 1 + 0, 2, 4, // i = 1, count = 2, 3, 4 + 1, 2, 3, // i = 2, count = 5, 6, 7 + 2, // i = 3, count = 8 + 1, 2, 4, // i = 4, count = 9, 10, 11 + // count = 12 + ]; + let correct_v = &[ + // p + 2.0, 3.0, // i = 0, count = 0, 1 + 3.0, 4.0, 6.0, // i = 1, count = 2, 3, 4 + -1.0, -3.0, 2.0, // i = 2, count = 5, 6, 7 + 1.0, // i = 3, count = 8 + 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 + // count = 12 + ]; + let correct_p = &[0, 2, 5, 8, 9, 12]; + + // check + assert_eq!(&csr.row_pointers, correct_p); + assert_eq!(&csr.col_indices, correct_j); + assert_eq!(&csr.values, correct_v); + Ok(()) +} diff --git a/russell_sparse/examples/doc_csr_from_tiny_coo.rs b/russell_sparse/examples/doc_csr_from_tiny_coo.rs new file mode 100644 index 00000000..3fc1c5d0 --- /dev/null +++ b/russell_sparse/examples/doc_csr_from_tiny_coo.rs @@ -0,0 +1,43 @@ +use russell_sparse::prelude::*; +use russell_sparse::StrError; + +fn main() -> Result<(), StrError> { + // allocate a square matrix and store as COO matrix + // ┌ ┐ + // │ 1 0 2 │ + // │ 0 0 3 │ << the diagonal 0 entry is optional, + // │ 4 5 6 │ but should be saved for Intel DSS + // └ ┘ + let (nrow, ncol, nnz) = (3, 3, 6); + let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?; + coo.put(0, 0, 1.0)?; + coo.put(0, 2, 2.0)?; + coo.put(1, 2, 3.0)?; + coo.put(2, 0, 4.0)?; + coo.put(2, 1, 5.0)?; + coo.put(2, 2, 6.0)?; + + // convert to CSR matrix + let csr = CsrMatrix::from(&coo); + let correct_v = &[ + // p + 1.0, 2.0, // i = 0, count = 0, 1 + 3.0, // i = 1, count = 2 + 4.0, 5.0, 6.0, // i = 2, count = 3, 4, 5 + // count = 6 + ]; + let correct_j = &[ + // p + 0, 2, // i = 0, count = 0, 1 + 2, // i = 1, count = 2 + 0, 1, 2, // i = 2, count = 3, 4, 5 + // count = 6 + ]; + let correct_p = &[0, 2, 3, 6]; + + // check + assert_eq!(&csr.row_pointers, correct_p); + assert_eq!(&csr.col_indices, correct_j); + assert_eq!(&csr.values, correct_v); + Ok(()) +} diff --git a/russell_sparse/examples/doc_csr_to_matrix.rs b/russell_sparse/examples/doc_csr_to_matrix.rs new file mode 100644 index 00000000..9f3810c6 --- /dev/null +++ b/russell_sparse/examples/doc_csr_to_matrix.rs @@ -0,0 +1,51 @@ +use russell_lab::Matrix; +use russell_sparse::prelude::*; +use russell_sparse::StrError; + +fn main() -> Result<(), StrError> { + // allocate a square matrix and store as CSR matrix + // ┌ ┐ + // │ 2 3 0 0 0 │ + // │ 3 0 4 0 6 │ + // │ 0 -1 -3 2 0 │ + // │ 0 0 1 0 0 │ + // │ 0 4 2 0 1 │ + // └ ┘ + let csr = CsrMatrix { + layout: Layout::Full, + nrow: 5, + ncol: 5, + row_pointers: vec![0, 2, 5, 8, 9, 12], + col_indices: vec![ + // p + 0, 1, // i = 0, count = 0, 1 + 0, 2, 4, // i = 1, count = 2, 3, 4 + 1, 2, 3, // i = 2, count = 5, 6, 7 + 2, // i = 3, count = 8 + 1, 2, 4, // i = 4, count = 9, 10, 11 + // count = 12 + ], + values: vec![ + // p + 2.0, 3.0, // i = 0, count = 0, 1 + 3.0, 4.0, 6.0, // i = 1, count = 2, 3, 4 + -1.0, -3.0, 2.0, // i = 2, count = 5, 6, 7 + 1.0, // i = 3, count = 8 + 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 + // count = 12 + ], + }; + + // covert to dense + let mut a = Matrix::new(5, 5); + csr.to_matrix(&mut a)?; + let correct = "┌ ┐\n\ + │ 2 3 0 0 0 │\n\ + │ 3 0 4 0 6 │\n\ + │ 0 -1 -3 2 0 │\n\ + │ 0 0 1 0 0 │\n\ + │ 0 4 2 0 1 │\n\ + └ ┘"; + assert_eq!(format!("{}", a), correct); + Ok(()) +} diff --git a/russell_sparse/examples/doc_solve_small_linear_system.rs b/russell_sparse/examples/doc_solve_small_linear_system.rs new file mode 100644 index 00000000..396a2db5 --- /dev/null +++ b/russell_sparse/examples/doc_solve_small_linear_system.rs @@ -0,0 +1,53 @@ +use russell_lab::{Matrix, Vector}; +use russell_sparse::prelude::*; +use russell_sparse::StrError; + +fn main() -> Result<(), StrError> { + // allocate a square matrix + 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(nrow, ncol); + coo.to_matrix(&mut a)?; + let correct = "┌ ┐\n\ + │ 2 3 0 0 0 │\n\ + │ 3 0 4 0 6 │\n\ + │ 0 -1 -3 2 0 │\n\ + │ 0 0 1 0 0 │\n\ + │ 0 4 2 0 1 │\n\ + └ ┘"; + assert_eq!(format!("{}", a), correct); + + // allocate x and rhs + 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, nrow, nnz, None)?; + solver.factorize(&coo)?; + solver.solve(&mut x, &rhs)?; + let correct = "┌ ┐\n\ + │ 1.000000 │\n\ + │ 2.000000 │\n\ + │ 3.000000 │\n\ + │ 4.000000 │\n\ + │ 5.000000 │\n\ + └ ┘"; + assert_eq!(format!("{:.6}", x), correct); + Ok(()) +} diff --git a/russell_sparse/examples/doc_solve_tiny_linear_system.rs b/russell_sparse/examples/doc_solve_tiny_linear_system.rs new file mode 100644 index 00000000..cb1317df --- /dev/null +++ b/russell_sparse/examples/doc_solve_tiny_linear_system.rs @@ -0,0 +1,51 @@ +use russell_lab::{Matrix, Vector}; +use russell_sparse::prelude::*; +use russell_sparse::StrError; + +fn main() -> Result<(), StrError> { + // allocate a square matrix + let nrow = 3; // number of equations + let ncol = nrow; // number of equations + let nnz = 5; // number of non-zeros + 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(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, &coo, &rhs1)?; + let correct1 = "┌ ┐\n\ + │ 3 │\n\ + │ 2 │\n\ + │ 4 │\n\ + └ ┘"; + assert_eq!(format!("{}", x1), correct1); + + // solve again + let mut x2 = Vector::new(nrow); + solver.solve(&mut x2, &rhs2)?; + let correct2 = "┌ ┐\n\ + │ 6 │\n\ + │ 4 │\n\ + │ 8 │\n\ + └ ┘"; + assert_eq!(format!("{}", x2), correct2); + Ok(()) +} diff --git a/russell_sparse/run_mem_check.bash b/russell_sparse/run_mem_check.bash deleted file mode 100755 index f446fc57..00000000 --- a/russell_sparse/run_mem_check.bash +++ /dev/null @@ -1,5 +0,0 @@ -#!/bin/bash - -cargo valgrind run --bin mem_check_build -cargo valgrind run --bin solve_mm_build -- data/matrix_market/bfwb62.mtx -cargo valgrind run --bin solve_mm_build -- data/matrix_market/bfwb62.mtx --mmp diff --git a/russell_sparse/run_solve_mm_build.bash b/russell_sparse/run_solve_mm_build.bash deleted file mode 100755 index ac6f3187..00000000 --- a/russell_sparse/run_solve_mm_build.bash +++ /dev/null @@ -1,4 +0,0 @@ -#!/bin/bash - -cargo run --release --bin solve_mm_build -- data/matrix_market/bfwb62.mtx - diff --git a/russell_sparse/src/bin/mem_check_build.rs b/russell_sparse/src/bin/mem_check_build.rs index 918725f0..2083b58a 100644 --- a/russell_sparse/src/bin/mem_check_build.rs +++ b/russell_sparse/src/bin/mem_check_build.rs @@ -1,5 +1,5 @@ use russell_lab::Vector; -use russell_sparse::{ConfigSolver, LinSolKind, Solver, SparseTriplet}; +use russell_sparse::prelude::*; fn test_solver(name: LinSolKind) { match name { @@ -7,33 +7,33 @@ fn test_solver(name: LinSolKind) { LinSolKind::Umf => println!("Testing UMF solver\n"), } - let (neq, nnz) = (5, 13); + let (nrow, ncol, nnz) = (5, 5, 13); - let mut trip = match SparseTriplet::new(neq, nnz) { + let mut coo = match CooMatrix::new(Layout::Full, nrow, ncol, nnz) { Ok(v) => v, Err(e) => { - println!("FAIL(new triplet): {}", e); + println!("FAIL(new CooMatrix): {}", e); return; } }; - trip.put(0, 0, 1.0).unwrap(); // << (0, 0, a00/2) - trip.put(0, 0, 1.0).unwrap(); // << (0, 0, a00/2) - trip.put(1, 0, 3.0).unwrap(); - trip.put(0, 1, 3.0).unwrap(); - trip.put(2, 1, -1.0).unwrap(); - trip.put(4, 1, 4.0).unwrap(); - trip.put(1, 2, 4.0).unwrap(); - trip.put(2, 2, -3.0).unwrap(); - trip.put(3, 2, 1.0).unwrap(); - trip.put(4, 2, 2.0).unwrap(); - trip.put(2, 3, 2.0).unwrap(); - trip.put(1, 4, 6.0).unwrap(); - trip.put(4, 4, 1.0).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(); + coo.put(0, 1, 3.0).unwrap(); + coo.put(2, 1, -1.0).unwrap(); + coo.put(4, 1, 4.0).unwrap(); + coo.put(1, 2, 4.0).unwrap(); + coo.put(2, 2, -3.0).unwrap(); + coo.put(3, 2, 1.0).unwrap(); + coo.put(4, 2, 2.0).unwrap(); + coo.put(2, 3, 2.0).unwrap(); + coo.put(1, 4, 6.0).unwrap(); + coo.put(4, 4, 1.0).unwrap(); let mut config = ConfigSolver::new(); config.lin_sol_kind(name); - let mut solver = match Solver::new(config, neq, nnz, None) { + let mut solver = match Solver::new(config, nrow, nnz, None) { Ok(v) => v, Err(e) => { println!("FAIL(new solver): {}", e); @@ -41,7 +41,7 @@ fn test_solver(name: LinSolKind) { } }; - match solver.factorize(&trip) { + match solver.factorize(&coo) { Err(e) => { println!("FAIL(factorize): {}", e); return; @@ -68,8 +68,6 @@ fn test_solver(name: LinSolKind) { _ => (), } - println!("{}", trip); - println!("{}", solver); println!("x =\n{}", x); } @@ -79,19 +77,19 @@ fn test_solver_singular(name: LinSolKind) { LinSolKind::Umf => println!("Testing UMF solver\n"), } - let (neq, nnz) = (2, 2); + let (nrow, ncol, nnz) = (2, 2, 2); - let trip_singular = match SparseTriplet::new(neq, nnz) { + let coo_singular = match CooMatrix::new(Layout::Full, nrow, ncol, nnz) { Ok(v) => v, Err(e) => { - println!("FAIL(new triplet): {}", e); + println!("FAIL(new CooMatrix): {}", e); return; } }; let mut config = ConfigSolver::new(); config.lin_sol_kind(name); - let mut solver = match Solver::new(config, neq, nnz, None) { + let mut solver = match Solver::new(config, nrow, nnz, None) { Ok(v) => v, Err(e) => { println!("FAIL(new solver): {}", e); @@ -99,7 +97,7 @@ fn test_solver_singular(name: LinSolKind) { } }; - match solver.factorize(&trip_singular) { + match solver.factorize(&coo_singular) { Err(e) => println!("\nOk(factorize singular matrix): {}\n", e), _ => (), }; diff --git a/russell_sparse/src/bin/solve_mm_build.rs b/russell_sparse/src/bin/solve_matrix_market_build.rs similarity index 55% rename from russell_sparse/src/bin/solve_mm_build.rs rename to russell_sparse/src/bin/solve_matrix_market_build.rs index a5566915..452337ac 100644 --- a/russell_sparse/src/bin/solve_mm_build.rs +++ b/russell_sparse/src/bin/solve_matrix_market_build.rs @@ -1,11 +1,44 @@ use russell_lab::{format_nanoseconds, Stopwatch, StrError, Vector}; -use russell_openblas::set_num_threads; -use russell_sparse::{ - enum_ordering, enum_scaling, read_matrix_market, ConfigSolver, LinSolKind, Solver, Symmetry, VerifyLinSys, -}; +use russell_openblas::{get_num_threads, set_num_threads}; +use russell_sparse::prelude::*; +use serde::{Deserialize, Serialize}; +use serde_json; use std::path::Path; use structopt::StructOpt; +/// Holds information about the solution of a linear system +#[derive(Clone, Debug, Deserialize, Serialize)] +pub struct SolutionInfo { + pub platform: String, + pub blas_lib: String, + pub solver_name: String, + pub matrix_name: String, + pub symmetry: String, + pub layout: String, + pub nrow: usize, + pub ncol: usize, + pub nnz: usize, + pub time_read_matrix_market_nanosecond: u128, + pub time_read_matrix_market_human: String, + pub time_factorize_nanosecond: u128, + pub time_factorize_human: String, + pub time_solve_nanosecond: u128, + pub time_solve_human: String, + pub time_factorize_and_solve_nanosecond: u128, + pub time_factorize_and_solve_human: String, + pub requested_ordering: String, + pub requested_scaling: String, + pub requested_openmp_num_threads: usize, + pub effective_ordering: String, + pub effective_scaling: String, + pub effective_openmp_num_threads: usize, + pub verify_max_abs_a: f64, + pub verify_max_abs_a_times_x: f64, + pub verify_relative_error: f64, + pub verify_time_nanosecond: u128, + pub verify_time_human: String, +} + /// Command line options #[derive(StructOpt, Debug)] #[structopt( @@ -49,25 +82,25 @@ fn main() -> Result<(), StrError> { // select linear solver let name = if opt.mmp { LinSolKind::Mmp } else { LinSolKind::Umf }; - // set the sym_mirror flag - let sym_mirror = match name { + // select the symmetric handling option + let handling = match name { LinSolKind::Mmp => { // MMP uses the lower-diagonal if symmetric. - false + SymmetricHandling::LeaveAsLower } LinSolKind::Umf => { // UMF uses the full matrix, if symmetric or not - true + SymmetricHandling::MakeItFull } }; // read matrix let mut sw = Stopwatch::new(""); - let (trip, symmetric) = read_matrix_market(&opt.matrix_market_file, sym_mirror)?; + let (coo, sym) = read_matrix_market(&opt.matrix_market_file, handling)?; let time_read = sw.stop(); // set the symmetry option - let symmetry = if symmetric { Some(Symmetry::General) } else { None }; + let symmetry = if sym { Some(Symmetry::General) } else { None }; // set configuration let mut config = ConfigSolver::new(); @@ -83,62 +116,70 @@ fn main() -> Result<(), StrError> { } // initialize and factorize - let (neq, nnz) = (trip.neq(), trip.nnz_current()); - let mut solver = Solver::new(config, neq, nnz, symmetry)?; - solver.factorize(&trip)?; + let (nrow, nnz) = (coo.nrow, coo.pos); + let mut solver = Solver::new(config, nrow, nnz, symmetry)?; + solver.factorize(&coo)?; // allocate vectors - let mut x = Vector::new(neq); - let rhs = Vector::filled(neq, 1.0); + let mut x = Vector::new(nrow); + let rhs = Vector::filled(nrow, 1.0); // solve linear system solver.solve(&mut x, &rhs)?; // verify solution - let triangular = symmetric && !sym_mirror; - let verify = VerifyLinSys::new(&trip, &x, &rhs, triangular)?; + let verify = VerifyLinSys::new(&coo, &x, &rhs)?; // matrix name let path = Path::new(&opt.matrix_market_file); - let matrix_name = path.file_stem().unwrap().to_str().unwrap(); + let matrix_name = match path.file_stem() { + Some(v) => match v.to_str() { + Some(w) => w.to_string(), + None => "Unknown".to_string(), + }, + None => "Unknown".to_string(), + }; // output - println!( - "{{\n\ - \x20\x20\"platform\": \"russell\",\n\ - \x20\x20\"blasLib\": \"OpenBLAS\",\n\ - \x20\x20\"matrixName\": \"{}\",\n\ - \x20\x20\"read\": {{\n\ - \x20\x20\x20\x20\"timeReadNs\": {},\n\ - \x20\x20\x20\x20\"timeReadStr\": \"{}\"\n\ - \x20\x20}},\n\ - \x20\x20\"triplet\": {{\n\ - {}\n\ - \x20\x20}},\n\ - \x20\x20\"symmetry\": \"{:?}\"\n\ - \x20\x20\"solver\": {{\n\ - {}\n\ - {}\n\ - \x20\x20}},\n\ - \x20\x20\"verify\": {{\n\ - {}\n\ - \x20\x20}}\n\ - }}", + let (time_fact, time_solve) = solver.get_elapsed_times(); + let info = SolutionInfo { + platform: "Russell".to_string(), + blas_lib: "OpenBLAS".to_string(), + solver_name: config.str_solver(), matrix_name, - time_read, - format_nanoseconds(time_read), - trip, - symmetry, - config, - solver, - verify - ); + symmetry: if sym { "General".to_string() } else { "None".to_string() }, + layout: format!("{:?}", coo.layout), + nrow: coo.nrow, + ncol: coo.ncol, + nnz: coo.pos, + time_read_matrix_market_nanosecond: time_read, + time_read_matrix_market_human: format_nanoseconds(time_read), + time_factorize_nanosecond: time_fact, + time_factorize_human: format_nanoseconds(time_fact), + time_solve_nanosecond: time_solve, + time_solve_human: format_nanoseconds(time_solve), + time_factorize_and_solve_nanosecond: time_fact + time_solve, + time_factorize_and_solve_human: format_nanoseconds(time_fact + time_solve), + requested_ordering: config.str_ordering(), + requested_scaling: config.str_scaling(), + requested_openmp_num_threads: opt.omp_nt as usize, + effective_ordering: solver.get_effective_ordering(), + effective_scaling: solver.get_effective_scaling(), + effective_openmp_num_threads: get_num_threads() as usize, + verify_max_abs_a: verify.max_abs_a, + verify_max_abs_a_times_x: verify.max_abs_ax, + verify_relative_error: verify.relative_error, + verify_time_nanosecond: verify.time_check, + verify_time_human: format_nanoseconds(verify.time_check), + }; + let info_json = serde_json::to_string_pretty(&info).unwrap(); + println!("{}", info_json); // check if path.ends_with("bfwb62.mtx") { let tolerance = if opt.mmp { 1e-10 } else { 1e-11 }; let correct_x = get_bfwb62_correct_x(); - for i in 0..neq { + for i in 0..nrow { let diff = f64::abs(x.get(i) - correct_x.get(i)); if diff > tolerance { println!("ERROR: diff({}) = {:.2e}", i, diff); diff --git a/russell_sparse/src/config_solver.rs b/russell_sparse/src/config_solver.rs index 82f1437f..754163e7 100644 --- a/russell_sparse/src/config_solver.rs +++ b/russell_sparse/src/config_solver.rs @@ -1,6 +1,5 @@ use super::{str_enum_ordering, str_enum_scaling, LinSolKind, Ordering, Scaling}; use russell_openblas::to_i32; -use std::fmt; /// Holds configuration options for the sparse Solver #[derive(Copy, Clone, Debug)] @@ -69,37 +68,29 @@ impl ConfigSolver { self.verbose = 1; self } -} -impl fmt::Display for ConfigSolver { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - let name = match self.lin_sol_kind { + /// Returns a string representation of the ordering option + pub fn str_ordering(&self) -> String { + str_enum_ordering(self.ordering).to_string() + } + + /// Returns a string representation of the scaling option + pub fn str_scaling(&self) -> String { + str_enum_scaling(self.scaling).to_string() + } + + /// Returns the name of the solver + pub fn str_solver(&self) -> String { + match self.lin_sol_kind { LinSolKind::Mmp => { if cfg!(local_mmp) { - "MMP-local" + "MUMPS-local".to_string() } else { - "MMP" + "MUMPS".to_string() } } - LinSolKind::Umf => "UMF", - }; - write!( - f, - "\x20\x20\x20\x20\"name\": \"{}\",\n\ - \x20\x20\x20\x20\"ordering\": \"{}\",\n\ - \x20\x20\x20\x20\"scaling\": \"{}\",\n\ - \x20\x20\x20\x20\"pctIncWorkspace\": {},\n\ - \x20\x20\x20\x20\"maxWorkMemory\": {},\n\ - \x20\x20\x20\x20\"openmpNumThreads\": {}", - name, - str_enum_ordering(self.ordering), - str_enum_scaling(self.scaling), - self.pct_inc_workspace, - self.max_work_memory, - self.openmp_num_threads, - ) - .unwrap(); - Ok(()) + LinSolKind::Umf => "UMFPACK".to_string(), + } } } @@ -186,32 +177,9 @@ mod tests { } #[test] - fn display_trait_works() { - let config1 = ConfigSolver::new(); - let correct1: &str = "\x20\x20\x20\x20\"name\": \"UMF\",\n\ - \x20\x20\x20\x20\"ordering\": \"Auto\",\n\ - \x20\x20\x20\x20\"scaling\": \"Auto\",\n\ - \x20\x20\x20\x20\"pctIncWorkspace\": 100,\n\ - \x20\x20\x20\x20\"maxWorkMemory\": 0,\n\ - \x20\x20\x20\x20\"openmpNumThreads\": 1"; - assert_eq!(format!("{}", config1), correct1); - let mut config2 = ConfigSolver::new(); - config2.lin_sol_kind(LinSolKind::Mmp); - let correct2: &str = if cfg!(local_mmp) { - "\x20\x20\x20\x20\"name\": \"MMP-local\",\n\ - \x20\x20\x20\x20\"ordering\": \"Auto\",\n\ - \x20\x20\x20\x20\"scaling\": \"Auto\",\n\ - \x20\x20\x20\x20\"pctIncWorkspace\": 100,\n\ - \x20\x20\x20\x20\"maxWorkMemory\": 0,\n\ - \x20\x20\x20\x20\"openmpNumThreads\": 1" - } else { - "\x20\x20\x20\x20\"name\": \"MMP\",\n\ - \x20\x20\x20\x20\"ordering\": \"Auto\",\n\ - \x20\x20\x20\x20\"scaling\": \"Auto\",\n\ - \x20\x20\x20\x20\"pctIncWorkspace\": 100,\n\ - \x20\x20\x20\x20\"maxWorkMemory\": 0,\n\ - \x20\x20\x20\x20\"openmpNumThreads\": 1" - }; - assert_eq!(format!("{}", config2), correct2); + fn string_methods_work() { + let config = ConfigSolver::new(); + assert_eq!(config.str_ordering(), "Auto"); + assert_eq!(config.str_scaling(), "Auto"); } } diff --git a/russell_sparse/src/coo_matrix.rs b/russell_sparse/src/coo_matrix.rs new file mode 100644 index 00000000..66fe32a3 --- /dev/null +++ b/russell_sparse/src/coo_matrix.rs @@ -0,0 +1,660 @@ +use super::Layout; +use crate::StrError; +use russell_lab::{Matrix, Vector}; +use russell_openblas::to_i32; + +/// Holds the row index, col index, and values of a matrix (also known as Triplet) +/// +/// # Remarks +/// +/// * Only the non-zero values are required +/// * Entries with repeated (i,j) indices are allowed +/// * Repeated (i,j) entries will have the aij values summed when solving a linear system +/// * The repeated (i,j) capability is of great convenience for Finite Element solvers +/// * A maximum number of entries must be decided prior to allocating a new Triplet +/// * The maximum number of entries includes possible entries with repeated indices +/// * See the `to_matrix` method for an example +pub struct CooMatrix { + /// Defines the stored layout: lower-triangular, upper-triangular, full-matrix + pub layout: Layout, + + /// Holds the number of rows (must fit i32) + pub nrow: usize, + + /// Holds the number of columns (must fit i32) + pub ncol: usize, + + /// Holds the current index (must fit i32) + /// + /// This will equal the number of non-zeros (nnz) after all items have been `put`. + pub pos: usize, + + /// Defines the maximum allowed number of entries (must fit i32) + /// + /// This may be greater than the number of non-zeros (nnz) + pub max: usize, + + /// Holds the row indices i + pub indices_i: Vec, + + /// Holds the column indices j + pub indices_j: Vec, + + /// Holds the values aij + pub values_aij: Vec, +} + +impl CooMatrix { + /// Creates a new CooMatrix representing a sparse matrix + /// + /// # Input + /// + /// * `layout` -- Defines the layout of the associated matrix + /// * `nrow` -- Is the number of rows of the sparse matrix (must be fit i32) + /// * `ncol` -- Is the number of columns of the sparse matrix (must be fit i32) + /// * `max` -- Maximum number of entries ≥ nnz (number of non-zeros), + /// including entries with repeated indices. (must be fit i32) + /// + /// # Example + /// + /// ``` + /// use russell_chk::vec_approx_eq; + /// use russell_sparse::{CooMatrix, Layout, StrError}; + /// + /// fn main() -> Result<(), StrError> { + /// let (nrow, ncol, nnz) = (3, 3, 4); + /// let coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?; + /// assert_eq!(coo.layout, Layout::Full); + /// assert_eq!(coo.pos, 0); + /// assert_eq!(coo.max, 4); + /// assert_eq!(coo.indices_i, &[0, 0, 0, 0]); + /// assert_eq!(coo.indices_j, &[0, 0, 0, 0]); + /// vec_approx_eq(&coo.values_aij, &[0.0, 0.0, 0.0, 0.0], 1e-15); + /// Ok(()) + /// } + /// ``` + pub fn new(layout: Layout, nrow: usize, ncol: usize, max: usize) -> Result { + if nrow < 1 { + return Err("nrow must be greater than zero"); + } + if ncol < 1 { + return Err("ncol must be greater than zero"); + } + if max < 1 { + return Err("max must be greater than zero"); + } + Ok(CooMatrix { + layout, + nrow, + ncol, + pos: 0, + max, + indices_i: vec![0; max], + indices_j: vec![0; max], + values_aij: vec![0.0; max], + }) + } + + /// Puts a new entry and updates pos (may be duplicate) + /// + /// # Input + /// + /// * `i` -- row index (indices start at zero; zero-based) + /// * `j` -- column index (indices start at zero; zero-based) + /// * `aij` -- the value A(i,j) + /// + /// # Example + /// + /// ``` + /// use russell_chk::vec_approx_eq; + /// use russell_sparse::{CooMatrix, Layout, StrError}; + /// + /// fn main() -> Result<(), StrError> { + /// let (nrow, ncol, nnz) = (3, 3, 4); + /// let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?; + /// coo.put(0, 0, 1.0)?; + /// coo.put(1, 1, 2.0)?; + /// coo.put(2, 2, 3.0)?; + /// coo.put(0, 1, 4.0)?; + /// assert_eq!(coo.layout, Layout::Full); + /// assert_eq!(coo.pos, 4); + /// assert_eq!(coo.max, 4); + /// assert_eq!(coo.indices_i, &[0, 1, 2, 0]); + /// assert_eq!(coo.indices_j, &[0, 1, 2, 1]); + /// vec_approx_eq(&coo.values_aij, &[1.0, 2.0, 3.0, 4.0], 1e-15); + /// Ok(()) + /// } + /// ``` + pub fn put(&mut self, i: usize, j: usize, aij: f64) -> Result<(), StrError> { + // check range + if i >= self.nrow { + return Err("index of row is outside range"); + } + if j >= self.ncol { + return Err("index of column is outside range"); + } + if self.pos >= self.max { + return Err("max number of items has been exceeded"); + } + if self.layout == Layout::Lower { + if j > i { + return Err("j > i is incorrect for lower triangular layout"); + } + } else if self.layout == Layout::Upper { + if j < i { + return Err("j < i is incorrect for upper triangular layout"); + } + } + + // insert a new entry + let i_i32 = to_i32(i); + let j_i32 = to_i32(j); + self.indices_i[self.pos] = i_i32; + self.indices_j[self.pos] = j_i32; + self.values_aij[self.pos] = aij; + self.pos += 1; + Ok(()) + } + + /// Resets the position of the current non-zero value + /// + /// This function allows using `put` all over again. + /// + /// # Example + /// + /// ``` + /// use russell_sparse::{CooMatrix, Layout, StrError}; + /// + /// fn main() -> Result<(), StrError> { + /// let (nrow, ncol, nnz) = (3, 3, 4); + /// let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?; + /// coo.put(0, 0, 1.0)?; + /// coo.put(1, 1, 2.0)?; + /// coo.put(2, 2, 3.0)?; + /// coo.put(0, 1, 4.0)?; + /// assert_eq!(coo.pos, 4); + /// coo.reset(); + /// assert_eq!(coo.pos, 0); + /// Ok(()) + /// } + /// ``` + pub fn reset(&mut self) { + self.pos = 0; + } + + /// Converts this CooMatrix to a dense matrix + /// + /// ``` + /// use russell_sparse::prelude::*; + /// use russell_sparse::StrError; + /// + /// fn main() -> Result<(), StrError> { + /// // define (4 x 4) sparse matrix with 6+1 non-zero values + /// // (with an extra ij-repeated entry) + /// let (nrow, ncol, nnz) = (4, 4, 7); + /// let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?; + /// coo.put(0, 0, 0.5)?; // (0, 0, a00/2) + /// coo.put(0, 0, 0.5)?; // (0, 0, a00/2) + /// coo.put(0, 1, 2.0)?; + /// coo.put(1, 0, 3.0)?; + /// coo.put(1, 1, 4.0)?; + /// coo.put(2, 2, 5.0)?; + /// coo.put(3, 3, 6.0)?; + /// + /// // convert to dense + /// let a = coo.as_matrix(); + /// let correct = "┌ ┐\n\ + /// │ 1 2 0 0 │\n\ + /// │ 3 4 0 0 │\n\ + /// │ 0 0 5 0 │\n\ + /// │ 0 0 0 6 │\n\ + /// └ ┘"; + /// assert_eq!(format!("{}", a), correct); + /// Ok(()) + /// } + /// ``` + pub fn as_matrix(&self) -> Matrix { + let mut a = Matrix::new(self.nrow, self.ncol); + self.to_matrix(&mut a).unwrap(); + a + } + + /// Converts this CooMatrix to a dense matrix + /// + /// # Input + /// + /// * `a` -- where to store the dense matrix; must be (nrow, ncol) + /// + /// # Examples + /// + /// ``` + /// use russell_lab::Matrix; + /// use russell_sparse::prelude::*; + /// use russell_sparse::StrError; + /// + /// fn main() -> Result<(), StrError> { + /// // define (4 x 4) sparse matrix with 6+1 non-zero values + /// // (with an extra ij-repeated entry) + /// let (nrow, ncol, nnz) = (4, 4, 7); + /// let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?; + /// coo.put(0, 0, 0.5)?; // (0, 0, a00/2) + /// coo.put(0, 0, 0.5)?; // (0, 0, a00/2) + /// coo.put(0, 1, 2.0)?; + /// coo.put(1, 0, 3.0)?; + /// coo.put(1, 1, 4.0)?; + /// coo.put(2, 2, 5.0)?; + /// coo.put(3, 3, 6.0)?; + /// + /// // convert to dense + /// let mut a = Matrix::new(4, 4); + /// coo.to_matrix(&mut a)?; + /// let correct = "┌ ┐\n\ + /// │ 1 2 0 0 │\n\ + /// │ 3 4 0 0 │\n\ + /// │ 0 0 5 0 │\n\ + /// │ 0 0 0 6 │\n\ + /// └ ┘"; + /// assert_eq!(format!("{}", a), correct); + /// Ok(()) + /// } + /// ``` + pub fn to_matrix(&self, a: &mut Matrix) -> Result<(), StrError> { + let (m, n) = a.dims(); + if m != self.nrow || n != self.ncol { + return Err("wrong matrix dimensions"); + } + a.fill(0.0); + for p in 0..self.pos { + let i = self.indices_i[p] as usize; + let j = self.indices_j[p] as usize; + a.add(i, j, self.values_aij[p]); + if self.layout != Layout::Full && i != j { + a.add(j, i, self.values_aij[p]); + } + } + Ok(()) + } + + /// Performs the matrix-vector multiplication + /// + /// ```text + /// v := a ⋅ u + /// (m) (m,n) (n) + /// ``` + /// + /// # Input + /// + /// * `u` -- Vector with dimension equal to the number of columns of the matrix + /// + /// # Output + /// + /// * `v` -- Vector with dimension equal to the number of rows of the matrix + /// + /// # Note + /// + /// This method is not highly efficient but should useful in verifications. + /// + /// # Example + /// + /// ``` + /// use russell_lab::{Matrix, Vector}; + /// use russell_sparse::{CooMatrix, Layout, StrError}; + /// + /// fn main() -> Result<(), StrError> { + /// // set sparse matrix (3 x 3) with 6 non-zeros + /// let (nrow, ncol, nnz) = (3, 3, 6); + /// let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?; + /// coo.put(0, 0, 1.0)?; + /// coo.put(1, 0, 2.0)?; + /// coo.put(1, 1, 3.0)?; + /// coo.put(2, 0, 4.0)?; + /// coo.put(2, 1, 5.0)?; + /// coo.put(2, 2, 6.0)?; + /// + /// // check matrix + /// let mut a = Matrix::new(nrow, ncol); + /// coo.to_matrix(&mut a)?; + /// let correct_a = "┌ ┐\n\ + /// │ 1 0 0 │\n\ + /// │ 2 3 0 │\n\ + /// │ 4 5 6 │\n\ + /// └ ┘"; + /// assert_eq!(format!("{}", a), correct_a); + /// + /// // perform mat-vec-mul + /// let u = Vector::from(&[1.0, 1.0, 1.0]); + /// let v = coo.mat_vec_mul(&u)?; + /// + /// // check vector + /// let correct_v = "┌ ┐\n\ + /// │ 1 │\n\ + /// │ 5 │\n\ + /// │ 15 │\n\ + /// └ ┘"; + /// assert_eq!(format!("{}", v), correct_v); + /// Ok(()) + /// } + /// ``` + pub fn mat_vec_mul(&self, u: &Vector) -> Result { + if u.dim() != self.ncol { + return Err("u.ndim must equal ncol"); + } + let mut v = Vector::new(self.nrow); + for p in 0..self.pos { + let i = self.indices_i[p] as usize; + let j = self.indices_j[p] as usize; + let aij = self.values_aij[p]; + v[i] += aij * u[j]; + if self.layout != Layout::Full && i != j { + v[j] += aij * u[i]; + } + } + Ok(v) + } +} + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +#[cfg(test)] +mod tests { + use super::CooMatrix; + use crate::Layout; + use russell_chk::vec_approx_eq; + use russell_lab::{Matrix, Vector}; + + #[test] + fn new_fails_on_wrong_input() { + assert_eq!( + CooMatrix::new(Layout::Full, 0, 1, 3).err(), + Some("nrow must be greater than zero") + ); + assert_eq!( + CooMatrix::new(Layout::Full, 1, 0, 3).err(), + Some("ncol must be greater than zero") + ); + assert_eq!( + CooMatrix::new(Layout::Full, 1, 1, 0).err(), + Some("max must be greater than zero") + ); + } + + #[test] + fn new_works() { + let coo = CooMatrix::new(Layout::Full, 1, 1, 3).unwrap(); + assert_eq!(coo.layout, Layout::Full); + assert_eq!(coo.nrow, 1); + assert_eq!(coo.ncol, 1); + assert_eq!(coo.pos, 0); + assert_eq!(coo.max, 3); + assert_eq!(coo.indices_i.len(), 3); + assert_eq!(coo.indices_j.len(), 3); + assert_eq!(coo.values_aij.len(), 3); + } + + #[test] + fn put_fails_on_wrong_values() { + let mut coo = CooMatrix::new(Layout::Full, 1, 1, 1).unwrap(); + assert_eq!(coo.put(1, 0, 0.0).err(), Some("index of row is outside range")); + assert_eq!(coo.put(0, 1, 0.0).err(), Some("index of column is outside range")); + assert_eq!(coo.put(0, 0, 0.0).err(), None); // << will take all spots + assert_eq!(coo.put(0, 0, 0.0).err(), Some("max number of items has been exceeded")); + let mut coo = CooMatrix::new(Layout::Lower, 2, 2, 4).unwrap(); + assert_eq!( + coo.put(0, 1, 0.0).err(), + Some("j > i is incorrect for lower triangular layout") + ); + let mut coo = CooMatrix::new(Layout::Upper, 2, 2, 4).unwrap(); + assert_eq!( + coo.put(1, 0, 0.0).err(), + Some("j < i is incorrect for upper triangular layout") + ); + } + + #[test] + fn put_works() { + let mut coo = CooMatrix::new(Layout::Full, 3, 3, 5).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + assert_eq!(coo.pos, 1); + coo.put(0, 1, 2.0).unwrap(); + assert_eq!(coo.pos, 2); + coo.put(1, 0, 3.0).unwrap(); + assert_eq!(coo.pos, 3); + coo.put(1, 1, 4.0).unwrap(); + assert_eq!(coo.pos, 4); + coo.put(2, 2, 5.0).unwrap(); + assert_eq!(coo.pos, 5); + } + + #[test] + fn reset_works() { + let mut coo = CooMatrix::new(Layout::Full, 2, 2, 4).unwrap(); + assert_eq!(coo.pos, 0); + coo.put(0, 0, 1.0).unwrap(); + coo.put(0, 1, 4.0).unwrap(); + coo.put(1, 0, 2.0).unwrap(); + coo.put(1, 1, 3.0).unwrap(); + assert_eq!(coo.pos, 4); + coo.reset(); + assert_eq!(coo.pos, 0); + } + + #[test] + fn to_matrix_fails_on_wrong_dims() { + let coo = CooMatrix::new(Layout::Full, 1, 1, 1).unwrap(); + let mut a_2x1 = Matrix::new(2, 1); + let mut a_1x2 = Matrix::new(1, 2); + assert_eq!(coo.to_matrix(&mut a_2x1), Err("wrong matrix dimensions")); + assert_eq!(coo.to_matrix(&mut a_1x2), Err("wrong matrix dimensions")); + } + + #[test] + fn to_matrix_works() { + let mut coo = CooMatrix::new(Layout::Full, 3, 3, 5).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(0, 1, 2.0).unwrap(); + coo.put(1, 0, 3.0).unwrap(); + coo.put(1, 1, 4.0).unwrap(); + coo.put(2, 2, 5.0).unwrap(); + let mut a = Matrix::new(3, 3); + coo.to_matrix(&mut a).unwrap(); + assert_eq!(a.get(0, 0), 1.0); + assert_eq!(a.get(0, 1), 2.0); + assert_eq!(a.get(1, 0), 3.0); + assert_eq!(a.get(1, 1), 4.0); + assert_eq!(a.get(2, 2), 5.0); + // using as_matrix + let bb = coo.as_matrix(); + assert_eq!(bb.get(0, 0), 1.0); + assert_eq!(bb.get(1, 0), 3.0); + } + + #[test] + fn to_matrix_with_duplicates_works() { + // allocate a square matrix + let (nrow, ncol, nnz) = (5, 5, 13); + let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz).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(); + coo.put(0, 1, 3.0).unwrap(); + coo.put(2, 1, -1.0).unwrap(); + coo.put(4, 1, 4.0).unwrap(); + coo.put(1, 2, 4.0).unwrap(); + coo.put(2, 2, -3.0).unwrap(); + coo.put(3, 2, 1.0).unwrap(); + coo.put(4, 2, 2.0).unwrap(); + coo.put(2, 3, 2.0).unwrap(); + coo.put(1, 4, 6.0).unwrap(); + coo.put(4, 4, 1.0).unwrap(); + + // print matrix + let mut a = Matrix::new(nrow as usize, ncol as usize); + coo.to_matrix(&mut a).unwrap(); + let correct = "┌ ┐\n\ + │ 2 3 0 0 0 │\n\ + │ 3 0 4 0 6 │\n\ + │ 0 -1 -3 2 0 │\n\ + │ 0 0 1 0 0 │\n\ + │ 0 4 2 0 1 │\n\ + └ ┘"; + assert_eq!(format!("{}", a), correct); + } + + #[test] + fn to_matrix_symmetric_lower_works() { + let mut coo = CooMatrix::new(Layout::Lower, 3, 3, 4).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(1, 0, 2.0).unwrap(); + coo.put(1, 1, 3.0).unwrap(); + coo.put(2, 1, 4.0).unwrap(); + let mut a = Matrix::new(3, 3); + coo.to_matrix(&mut a).unwrap(); + let correct = "┌ ┐\n\ + │ 1 2 0 │\n\ + │ 2 3 4 │\n\ + │ 0 4 0 │\n\ + └ ┘"; + assert_eq!(format!("{}", a), correct); + } + + #[test] + fn to_matrix_symmetric_upper_works() { + let mut coo = CooMatrix::new(Layout::Upper, 3, 3, 4).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(0, 1, 2.0).unwrap(); + coo.put(1, 1, 3.0).unwrap(); + coo.put(1, 2, 4.0).unwrap(); + let mut a = Matrix::new(3, 3); + coo.to_matrix(&mut a).unwrap(); + let correct = "┌ ┐\n\ + │ 1 2 0 │\n\ + │ 2 3 4 │\n\ + │ 0 4 0 │\n\ + └ ┘"; + assert_eq!(format!("{}", a), correct); + } + + #[test] + fn mat_vec_mul_fails_on_wrong_input() { + let coo = CooMatrix::new(Layout::Full, 2, 2, 1).unwrap(); + let u = Vector::new(3); + assert_eq!(coo.mat_vec_mul(&u).err(), Some("u.ndim must equal ncol")); + } + + #[test] + fn mat_vec_mul_works() { + // 1.0 2.0 3.0 + // 0.1 0.2 0.3 + // 10.0 20.0 30.0 + let mut coo = CooMatrix::new(Layout::Full, 3, 3, 9).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(0, 1, 2.0).unwrap(); + coo.put(0, 2, 3.0).unwrap(); + coo.put(1, 0, 0.1).unwrap(); + coo.put(1, 1, 0.2).unwrap(); + coo.put(1, 2, 0.3).unwrap(); + coo.put(2, 0, 10.0).unwrap(); + coo.put(2, 1, 20.0).unwrap(); + coo.put(2, 2, 30.0).unwrap(); + let u = Vector::from(&[0.1, 0.2, 0.3]); + let correct_v = &[1.4, 0.14, 14.0]; + let v = coo.mat_vec_mul(&u).unwrap(); + vec_approx_eq(v.as_data(), correct_v, 1e-15); + } + + #[test] + fn mat_vec_mul_symmetric_lower_works() { + // 2 + // 1 2 sym + // 1 2 9 + // 3 1 1 7 + // 2 1 5 1 8 + let (nrow, ncol, nnz) = (5, 5, 15); + let mut coo = CooMatrix::new(Layout::Lower, nrow, ncol, nnz).unwrap(); + coo.put(0, 0, 2.0).unwrap(); + coo.put(1, 1, 2.0).unwrap(); + coo.put(2, 2, 9.0).unwrap(); + coo.put(3, 3, 7.0).unwrap(); + coo.put(4, 4, 8.0).unwrap(); + + coo.put(1, 0, 1.0).unwrap(); + + coo.put(2, 0, 1.0).unwrap(); + coo.put(2, 1, 2.0).unwrap(); + + coo.put(3, 0, 3.0).unwrap(); + coo.put(3, 1, 1.0).unwrap(); + coo.put(3, 2, 1.0).unwrap(); + + coo.put(4, 0, 2.0).unwrap(); + coo.put(4, 1, 1.0).unwrap(); + coo.put(4, 2, 5.0).unwrap(); + coo.put(4, 3, 1.0).unwrap(); + let u = Vector::from(&[-629.0 / 98.0, 237.0 / 49.0, -53.0 / 49.0, 62.0 / 49.0, 23.0 / 14.0]); + let correct_v = &[-2.0, 4.0, 3.0, -5.0, 1.0]; + let v = coo.mat_vec_mul(&u).unwrap(); + vec_approx_eq(v.as_data(), correct_v, 1e-14); + } + + #[test] + fn mat_vec_mul_symmetric_full_works() { + // 2 1 1 3 2 + // 1 2 2 1 1 + // 1 2 9 1 5 + // 3 1 1 7 1 + // 2 1 5 1 8 + let (nrow, ncol, nnz) = (5, 5, 25); + let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz).unwrap(); + coo.put(0, 0, 2.0).unwrap(); + coo.put(1, 1, 2.0).unwrap(); + coo.put(2, 2, 9.0).unwrap(); + coo.put(3, 3, 7.0).unwrap(); + coo.put(4, 4, 8.0).unwrap(); + + coo.put(1, 0, 1.0).unwrap(); + coo.put(0, 1, 1.0).unwrap(); + + coo.put(2, 0, 1.0).unwrap(); + coo.put(0, 2, 1.0).unwrap(); + coo.put(2, 1, 2.0).unwrap(); + coo.put(1, 2, 2.0).unwrap(); + + coo.put(3, 0, 3.0).unwrap(); + coo.put(0, 3, 3.0).unwrap(); + coo.put(3, 1, 1.0).unwrap(); + coo.put(1, 3, 1.0).unwrap(); + coo.put(3, 2, 1.0).unwrap(); + coo.put(2, 3, 1.0).unwrap(); + + coo.put(4, 0, 2.0).unwrap(); + coo.put(0, 4, 2.0).unwrap(); + coo.put(4, 1, 1.0).unwrap(); + coo.put(1, 4, 1.0).unwrap(); + coo.put(4, 2, 5.0).unwrap(); + coo.put(2, 4, 5.0).unwrap(); + coo.put(4, 3, 1.0).unwrap(); + coo.put(3, 4, 1.0).unwrap(); + let u = Vector::from(&[-629.0 / 98.0, 237.0 / 49.0, -53.0 / 49.0, 62.0 / 49.0, 23.0 / 14.0]); + let correct_v = &[-2.0, 4.0, 3.0, -5.0, 1.0]; + let v = coo.mat_vec_mul(&u).unwrap(); + vec_approx_eq(v.as_data(), correct_v, 1e-14); + } + + #[test] + fn mat_vec_mul_positive_definite_works() { + // 2 -1 2 ... + // -1 2 -1 => -1 2 + // -1 2 -1 2 + let (nrow, ncol, nnz) = (3, 3, 5); + let mut coo = CooMatrix::new(Layout::Lower, nrow, ncol, nnz).unwrap(); + coo.put(0, 0, 2.0).unwrap(); + coo.put(1, 1, 2.0).unwrap(); + coo.put(2, 2, 2.0).unwrap(); + coo.put(1, 0, -1.0).unwrap(); + coo.put(2, 1, -1.0).unwrap(); + let u = Vector::from(&[5.0, 8.0, 7.0]); + let correct_v = &[2.0, 4.0, 6.0]; + let v = coo.mat_vec_mul(&u).unwrap(); + vec_approx_eq(v.as_data(), correct_v, 1e-15); + } +} diff --git a/russell_sparse/src/csr_matrix.rs b/russell_sparse/src/csr_matrix.rs new file mode 100644 index 00000000..e204b14d --- /dev/null +++ b/russell_sparse/src/csr_matrix.rs @@ -0,0 +1,733 @@ +use super::{CooMatrix, Layout}; +use crate::StrError; +use russell_lab::Matrix; +use russell_openblas::to_i32; + +/// Holds the arrays needed for a CSR (compressed sparse row) matrix +pub struct CsrMatrix { + /// Defines the stored layout: lower-triangular, upper-triangular, full-matrix + pub layout: Layout, + + /// Holds the number of rows (must fit i32) + pub nrow: usize, + + /// Holds the number of columns (must fit i32) + pub ncol: usize, + + /// Defines the row pointers array with size = nrow + 1 + pub row_pointers: Vec, + + /// Defines the column indices array with size = nnz (number of non-zeros) + pub col_indices: Vec, + + /// Defines the values array with size = nnz (number of non-zeros) + pub values: Vec, +} + +impl CsrMatrix { + /// Allocates an empty CsrMatrix + /// + /// This function simply allocates the following arrays: + /// + /// * row_pointers: `vec![0; nrow + 1]` + /// * col_indices: `vec![0; nnz]` + /// * values: `vec![0.0; nnz]` + /// + /// # Examples + /// + /// ``` + /// use russell_sparse::prelude::*; + /// use russell_sparse::StrError; + /// + /// fn main() { + /// let csr = CsrMatrix::new(Layout::Full, 3, 3, 4); + /// assert_eq!(csr.layout, Layout::Full); + /// assert_eq!(csr.nrow, 3); + /// assert_eq!(csr.ncol, 3); + /// assert_eq!(csr.row_pointers, &[0, 0, 0, 0]); + /// assert_eq!(csr.col_indices, &[0, 0, 0, 0]); + /// assert_eq!(csr.values, &[0.0, 0.0, 0.0, 0.0]); + /// } + /// ``` + pub fn new(layout: Layout, nrow: usize, ncol: usize, nnz: usize) -> Self { + CsrMatrix { + layout, + nrow, + ncol, + row_pointers: vec![0; nrow + 1], + col_indices: vec![0; nnz], + values: vec![0.0; nnz], + } + } + + /// Creates a new CsrMatrix from a CooMatrix + /// + /// # Examples + /// + /// ``` + /// use russell_sparse::prelude::*; + /// use russell_sparse::StrError; + /// + /// fn main() -> Result<(), StrError> { + /// // allocate a square matrix and store as COO matrix + /// // ┌ ┐ + /// // │ 1 0 2 │ + /// // │ 0 0 3 │ << the diagonal 0 entry is optional, + /// // │ 4 5 6 │ but should be saved for Intel DSS + /// // └ ┘ + /// let (nrow, ncol, nnz) = (3, 3, 6); + /// let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?; + /// coo.put(0, 0, 1.0)?; + /// coo.put(0, 2, 2.0)?; + /// coo.put(1, 2, 3.0)?; + /// coo.put(2, 0, 4.0)?; + /// coo.put(2, 1, 5.0)?; + /// coo.put(2, 2, 6.0)?; + /// + /// // convert to CSR matrix + /// let csr = CsrMatrix::from(&coo); + /// let correct_v = &[ + /// // p + /// 1.0, 2.0, // i = 0, count = 0, 1 + /// 3.0, // i = 1, count = 2 + /// 4.0, 5.0, 6.0, // i = 2, count = 3, 4, 5 + /// // count = 6 + /// ]; + /// let correct_j = &[ + /// // p + /// 0, 2, // i = 0, count = 0, 1 + /// 2, // i = 1, count = 2 + /// 0, 1, 2, // i = 2, count = 3, 4, 5 + /// // count = 6 + /// ]; + /// let correct_p = &[0, 2, 3, 6]; + /// + /// // check + /// assert_eq!(&csr.row_pointers, correct_p); + /// assert_eq!(&csr.col_indices, correct_j); + /// assert_eq!(&csr.values, correct_v); + /// Ok(()) + /// } + /// ``` + pub fn from(coo: &CooMatrix) -> Self { + // Based on the SciPy code from here: + // + // https://github.com/scipy/scipy/blob/main/scipy/sparse/sparsetools/coo.h + // + // Notes: + // + // * The row and column indices may be unordered + // * Linear complexity: O(nnz(A) + max(nrow,ncol)) + + // access triplet data + let ai = &coo.indices_i; + let aj = &coo.indices_j; + let ax = &coo.values_aij; + + // allocate vectors + let nrow = coo.nrow; + let nnz = coo.pos; + let mut csr = CsrMatrix { + layout: coo.layout, + nrow: coo.nrow, + ncol: coo.ncol, + row_pointers: vec![0; nrow + 1], + col_indices: vec![0; nnz], + values: vec![0.0; nnz], + }; + let bp = &mut csr.row_pointers; + let bj = &mut csr.col_indices; + let bx = &mut csr.values; + + // compute number of non-zero entries per row of A + for k in 0..nnz { + bp[ai[k] as usize] += 1; + } + + // perform the cumulative sum of the nnz per row to get bp[] + let mut cum_sum = 0; + for i in 0..nrow { + let temp = bp[i]; + bp[i] = cum_sum; + cum_sum += temp; + } + bp[nrow] = to_i32(nnz); + + // write aj and ax into bj and bx (will use bp as workspace) + for k in 0..nnz { + let row = ai[k]; + let dest = bp[row as usize]; + bj[dest as usize] = aj[k]; + bx[dest as usize] = ax[k]; + bp[row as usize] += 1; + } + + // fix bp + let mut last = 0; + for i in 0..nrow { + let temp = bp[i]; + bp[i] = last; + last = temp; + } + + // sort rows + let mut temp: Vec<(i32, f64)> = Vec::new(); + for i in 0..nrow { + let row_start = bp[i]; + let row_end = bp[i + 1]; + temp.resize((row_end - row_start) as usize, (0, 0.0)); + let mut n = 0; + for jj in row_start..row_end { + temp[n].0 = bj[jj as usize]; + temp[n].1 = bx[jj as usize]; + n += 1; + } + temp.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap()); + n = 0; + for jj in row_start..row_end { + bj[jj as usize] = temp[n].0; + bx[jj as usize] = temp[n].1; + n += 1; + } + } + + // sum duplicates + let final_nnz = csr_sum_duplicates(nrow, bp, bj, bx); + bj.resize(final_nnz, 0); + bx.resize(final_nnz, 0.0); + + // results + csr + } + + /// Converts this CsrMatrix to a dense matrix + /// + /// # Examples + /// + /// ``` + /// use russell_sparse::prelude::*; + /// use russell_sparse::StrError; + /// + /// fn main() -> Result<(), StrError> { + /// // allocate a square matrix and store as CSR matrix + /// // ┌ ┐ + /// // │ 2 3 0 0 0 │ + /// // │ 3 0 4 0 6 │ + /// // │ 0 -1 -3 2 0 │ + /// // │ 0 0 1 0 0 │ + /// // │ 0 4 2 0 1 │ + /// // └ ┘ + /// let csr = CsrMatrix { + /// layout: Layout::Full, + /// nrow: 5, + /// ncol: 5, + /// row_pointers: vec![0, 2, 5, 8, 9, 12], + /// col_indices: vec![ + /// // p + /// 0, 1, // i = 0, count = 0, 1 + /// 0, 2, 4, // i = 1, count = 2, 3, 4 + /// 1, 2, 3, // i = 2, count = 5, 6, 7 + /// 2, // i = 3, count = 8 + /// 1, 2, 4, // i = 4, count = 9, 10, 11 + /// // count = 12 + /// ], + /// values: vec![ + /// // p + /// 2.0, 3.0, // i = 0, count = 0, 1 + /// 3.0, 4.0, 6.0, // i = 1, count = 2, 3, 4 + /// -1.0, -3.0, 2.0, // i = 2, count = 5, 6, 7 + /// 1.0, // i = 3, count = 8 + /// 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 + /// // count = 12 + /// ], + /// }; + /// + /// // covert to dense + /// let a = csr.as_matrix()?; + /// let correct = "┌ ┐\n\ + /// │ 2 3 0 0 0 │\n\ + /// │ 3 0 4 0 6 │\n\ + /// │ 0 -1 -3 2 0 │\n\ + /// │ 0 0 1 0 0 │\n\ + /// │ 0 4 2 0 1 │\n\ + /// └ ┘"; + /// assert_eq!(format!("{}", a), correct); + /// Ok(()) + /// } + /// ``` + pub fn as_matrix(&self) -> Result { + let mut a = Matrix::new(self.nrow, self.ncol); + self.to_matrix(&mut a).unwrap(); + Ok(a) + } + + /// Converts this CsrMatrix to a dense matrix + /// + /// # Input + /// + /// * `a` -- where to store the dense matrix; must be (nrow, ncol) + /// + /// # Examples + /// + /// ``` + /// use russell_lab::Matrix; + /// use russell_sparse::prelude::*; + /// use russell_sparse::StrError; + /// + /// fn main() -> Result<(), StrError> { + /// // allocate a square matrix and store as CSR matrix + /// // ┌ ┐ + /// // │ 2 3 0 0 0 │ + /// // │ 3 0 4 0 6 │ + /// // │ 0 -1 -3 2 0 │ + /// // │ 0 0 1 0 0 │ + /// // │ 0 4 2 0 1 │ + /// // └ ┘ + /// let csr = CsrMatrix { + /// layout: Layout::Full, + /// nrow: 5, + /// ncol: 5, + /// row_pointers: vec![0, 2, 5, 8, 9, 12], + /// col_indices: vec![ + /// // p + /// 0, 1, // i = 0, count = 0, 1 + /// 0, 2, 4, // i = 1, count = 2, 3, 4 + /// 1, 2, 3, // i = 2, count = 5, 6, 7 + /// 2, // i = 3, count = 8 + /// 1, 2, 4, // i = 4, count = 9, 10, 11 + /// // count = 12 + /// ], + /// values: vec![ + /// // p + /// 2.0, 3.0, // i = 0, count = 0, 1 + /// 3.0, 4.0, 6.0, // i = 1, count = 2, 3, 4 + /// -1.0, -3.0, 2.0, // i = 2, count = 5, 6, 7 + /// 1.0, // i = 3, count = 8 + /// 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 + /// // count = 12 + /// ], + /// }; + /// + /// // covert to dense + /// let mut a = Matrix::new(5, 5); + /// csr.to_matrix(&mut a)?; + /// let correct = "┌ ┐\n\ + /// │ 2 3 0 0 0 │\n\ + /// │ 3 0 4 0 6 │\n\ + /// │ 0 -1 -3 2 0 │\n\ + /// │ 0 0 1 0 0 │\n\ + /// │ 0 4 2 0 1 │\n\ + /// └ ┘"; + /// assert_eq!(format!("{}", a), correct); + /// Ok(()) + /// } + /// ``` + pub fn to_matrix(&self, a: &mut Matrix) -> Result<(), StrError> { + let (m, n) = a.dims(); + if m != self.nrow || n != self.ncol { + return Err("wrong matrix dimensions"); + } + a.fill(0.0); + for i in 0..self.nrow { + for p in self.row_pointers[i]..self.row_pointers[i + 1] { + let j = self.col_indices[p as usize] as usize; + a.add(i, j, self.values[p as usize]); + if self.layout != Layout::Full && i != j { + a.add(j, i, self.values[p as usize]); + } + } + } + Ok(()) + } +} + +/// brief Sums duplicate column entries in each row of a CSR matrix +/// +/// Returns The final number of non-zeros (nnz) after duplicates have been handled +fn csr_sum_duplicates(nrow: usize, ap: &mut [i32], aj: &mut [i32], ax: &mut [f64]) -> usize { + // Based on the SciPy code from here: + // + // https://github.com/scipy/scipy/blob/main/scipy/sparse/sparsetools/csr.h + // + // * ap[n_row+1] -- row pointer + // * aj[nnz(A)] -- column indices + // * ax[nnz(A)] -- non-zeros + // * The column indices within each row must be sorted + // * Explicit zeros are retained + // * ap, aj, and ax will be modified in place + + let mut nnz: i32 = 0; + let mut row_end = 0; + for i in 0..nrow { + let mut jj = row_end; + row_end = ap[i + 1]; + while jj < row_end { + let j = aj[jj as usize]; + let mut x = ax[jj as usize]; + jj += 1; + while jj < row_end && aj[jj as usize] == j { + x += ax[jj as usize]; + jj += 1; + } + aj[nnz as usize] = j; + ax[nnz as usize] = x; + nnz += 1; + } + ap[i + 1] = nnz; + } + nnz as usize +} + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +#[cfg(test)] +mod tests { + use super::CsrMatrix; + use crate::{CooMatrix, Layout}; + use russell_chk::vec_approx_eq; + use russell_lab::Matrix; + + #[test] + fn new_works() { + let csr = CsrMatrix::new(Layout::Full, 3, 3, 4); + assert_eq!(csr.layout, Layout::Full); + assert_eq!(csr.nrow, 3); + assert_eq!(csr.ncol, 3); + assert_eq!(csr.row_pointers, &[0, 0, 0, 0]); + assert_eq!(csr.col_indices, &[0, 0, 0, 0]); + assert_eq!(csr.values, &[0.0, 0.0, 0.0, 0.0]); + } + + #[test] + fn csr_matrix_first_triplet_with_shuffled_entries() { + // 1 -1 . -3 . + // -2 5 . . . + // . . 4 6 4 + // -4 . 2 7 . + // . 8 . . -5 + // first triplet with shuffled entries + let mut coo = CooMatrix::new(Layout::Full, 5, 5, 13).unwrap(); + coo.put(2, 4, 4.0).unwrap(); + coo.put(4, 1, 8.0).unwrap(); + coo.put(0, 1, -1.0).unwrap(); + coo.put(2, 2, 4.0).unwrap(); + coo.put(4, 4, -5.0).unwrap(); + coo.put(3, 0, -4.0).unwrap(); + coo.put(0, 3, -3.0).unwrap(); + coo.put(2, 3, 6.0).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(1, 1, 5.0).unwrap(); + coo.put(3, 2, 2.0).unwrap(); + coo.put(1, 0, -2.0).unwrap(); + coo.put(3, 3, 7.0).unwrap(); + let csr = CsrMatrix::from(&coo); + // solution + let correct_p = vec![0, 3, 5, 8, 11, 13]; + let correct_j = vec![ + 0, 1, 3, /**/ 0, 1, /**/ 2, 3, 4, /**/ 0, 2, 3, /**/ 1, 4, + ]; + let correct_x = vec![1.0, -1.0, -3.0, -2.0, 5.0, 4.0, 6.0, 4.0, -4.0, 2.0, 7.0, 8.0, -5.0]; + assert_eq!(&csr.row_pointers, &correct_p); + assert_eq!(&csr.col_indices, &correct_j); + vec_approx_eq(&csr.values, &correct_x, 1e-15); + } + + #[test] + fn csr_matrix_small_triplet_with_shuffled_entries() { + // 1 2 . . . + // 3 4 . . . + // . . 5 6 . + // . . 7 8 . + // . . . . 9 + // small triplet with shuffled entries + let mut coo = CooMatrix::new(Layout::Full, 5, 5, 9).unwrap(); + coo.put(4, 4, 9.0).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(1, 0, 3.0).unwrap(); + coo.put(2, 2, 5.0).unwrap(); + coo.put(2, 3, 6.0).unwrap(); + coo.put(0, 1, 2.0).unwrap(); + coo.put(3, 2, 7.0).unwrap(); + coo.put(1, 1, 4.0).unwrap(); + coo.put(3, 3, 8.0).unwrap(); + let csr = CsrMatrix::from(&coo); + // solution + let correct_p = vec![0, 2, 4, 6, 8, 9]; + let correct_j = vec![0, 1, 0, 1, 2, 3, 2, 3, 4]; + let correct_x = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]; + assert_eq!(&csr.row_pointers, &correct_p); + assert_eq!(&csr.col_indices, &correct_j); + vec_approx_eq(&csr.values, &correct_x, 1e-15); + } + + #[test] + fn csr_matrix_small_triplet_with_duplicates() { + // 1 2 . . . + // 3 4 . . . + // . . 5 6 . + // . . 7 8 . + // . . . . 9 + // with duplicates + let mut coo = CooMatrix::new(Layout::Full, 5, 5, 11).unwrap(); + coo.put(4, 4, 9.0).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(1, 0, 3.0).unwrap(); + coo.put(2, 2, 5.0).unwrap(); + coo.put(2, 3, 3.0).unwrap(); // << + coo.put(0, 1, 2.0).unwrap(); + coo.put(3, 2, 7.0).unwrap(); + coo.put(1, 1, 2.0).unwrap(); // << + coo.put(3, 3, 8.0).unwrap(); + coo.put(2, 3, 3.0).unwrap(); // << duplicate + coo.put(1, 1, 2.0).unwrap(); // << duplicate + let csr = CsrMatrix::from(&coo); + // solution + let correct_p = vec![0, 2, 4, 6, 8, 9]; + let correct_j = vec![0, 1, 0, 1, 2, 3, 2, 3, 4]; + let correct_x = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]; + assert_eq!(&csr.row_pointers, &correct_p); + assert_eq!(&csr.col_indices, &correct_j); + vec_approx_eq(&csr.values, &correct_x, 1e-15); + } + + #[test] + fn csr_matrix_upper_triangular_with_ordered_entries() { + // 9.00 1.5 6.0 0.750 3.0 + // 1.50 0.5 0.0 0.000 0.0 + // 6.00 0.0 12.0 0.000 0.0 + // 0.75 0.0 0.0 0.625 0.0 + // 3.00 0.0 0.0 0.000 16.0 + // upper triangular with ordered entries + let mut coo = CooMatrix::new(Layout::Upper, 5, 5, 9).unwrap(); + coo.put(0, 0, 9.0).unwrap(); + coo.put(0, 1, 1.5).unwrap(); + coo.put(1, 1, 0.5).unwrap(); + coo.put(0, 2, 6.0).unwrap(); + coo.put(2, 2, 12.0).unwrap(); + coo.put(0, 3, 0.75).unwrap(); + coo.put(3, 3, 0.625).unwrap(); + coo.put(0, 4, 3.0).unwrap(); + coo.put(4, 4, 16.0).unwrap(); + let csr = CsrMatrix::from(&coo); + // solution + let correct_p = vec![0, 5, 6, 7, 8, 9]; + let correct_j = vec![0, 1, 2, 3, 4, 1, 2, 3, 4]; + let correct_x = vec![9.0, 1.5, 6.0, 0.75, 3.0, 0.5, 12.0, 0.625, 16.0]; + assert_eq!(&csr.row_pointers, &correct_p); + assert_eq!(&csr.col_indices, &correct_j); + vec_approx_eq(&csr.values, &correct_x, 1e-15); + } + + #[test] + fn csr_matrix_upper_triangular_with_shuffled_entries() { + // 9.00 1.5 6.0 0.750 3.0 + // 1.50 0.5 0.0 0.000 0.0 + // 6.00 0.0 12.0 0.000 0.0 + // 0.75 0.0 0.0 0.625 0.0 + // 3.00 0.0 0.0 0.000 16.0 + // upper triangular with shuffled entries + let mut coo = CooMatrix::new(Layout::Upper, 5, 5, 9).unwrap(); + coo.put(2, 2, 12.0).unwrap(); + coo.put(0, 0, 9.0).unwrap(); + coo.put(3, 3, 0.625).unwrap(); + coo.put(0, 1, 1.5).unwrap(); + coo.put(0, 2, 6.0).unwrap(); + coo.put(4, 4, 16.0).unwrap(); + coo.put(0, 3, 0.75).unwrap(); + coo.put(1, 1, 0.5).unwrap(); + coo.put(0, 4, 3.0).unwrap(); + let csr = CsrMatrix::from(&coo); + // solution + let correct_p = vec![0, 5, 6, 7, 8, 9]; + let correct_j = vec![0, 1, 2, 3, 4, 1, 2, 3, 4]; + let correct_x = vec![9.0, 1.5, 6.0, 0.75, 3.0, 0.5, 12.0, 0.625, 16.0]; + assert_eq!(&csr.row_pointers, &correct_p); + assert_eq!(&csr.col_indices, &correct_j); + vec_approx_eq(&csr.values, &correct_x, 1e-15); + } + + #[test] + fn csr_matrix_upper_triangular_with_diagonal_entries_first() { + // 9.00 1.5 6.0 0.750 3.0 + // 1.50 0.5 0.0 0.000 0.0 + // 6.00 0.0 12.0 0.000 0.0 + // 0.75 0.0 0.0 0.625 0.0 + // 3.00 0.0 0.0 0.000 16.0 + // upper triangular with diagonal entries being set first + let mut coo = CooMatrix::new(Layout::Upper, 5, 5, 9).unwrap(); + // diagonal + coo.put(0, 0, 9.0).unwrap(); + coo.put(1, 1, 0.5).unwrap(); + coo.put(2, 2, 12.0).unwrap(); + coo.put(3, 3, 0.625).unwrap(); + coo.put(4, 4, 16.0).unwrap(); + // upper diagonal + coo.put(0, 1, 1.5).unwrap(); + coo.put(0, 2, 6.0).unwrap(); + coo.put(0, 3, 0.75).unwrap(); + coo.put(0, 4, 3.0).unwrap(); + let csr = CsrMatrix::from(&coo); + // solution + let correct_p = vec![0, 5, 6, 7, 8, 9]; + let correct_j = vec![0, 1, 2, 3, 4, 1, 2, 3, 4]; + let correct_x = vec![9.0, 1.5, 6.0, 0.75, 3.0, 0.5, 12.0, 0.625, 16.0]; + assert_eq!(&csr.row_pointers, &correct_p); + assert_eq!(&csr.col_indices, &correct_j); + vec_approx_eq(&csr.values, &correct_x, 1e-15); + } + + #[test] + fn csr_matrix_lower_triangular_with_ordered_entries() { + // 9.00 1.5 6.0 0.750 3.0 + // 1.50 0.5 0.0 0.000 0.0 + // 6.00 0.0 12.0 0.000 0.0 + // 0.75 0.0 0.0 0.625 0.0 + // 3.00 0.0 0.0 0.000 16.0 + // lower diagonal with ordered entries + let mut coo = CooMatrix::new(Layout::Lower, 5, 5, 9).unwrap(); + coo.put(0, 0, 9.0).unwrap(); + coo.put(1, 0, 1.5).unwrap(); + coo.put(1, 1, 0.5).unwrap(); + coo.put(2, 0, 6.0).unwrap(); + coo.put(2, 2, 12.0).unwrap(); + coo.put(3, 0, 0.75).unwrap(); + coo.put(3, 3, 0.625).unwrap(); + coo.put(4, 0, 3.0).unwrap(); + coo.put(4, 4, 16.0).unwrap(); + let csr = CsrMatrix::from(&coo); + // solution + let correct_p = vec![0, 1, 3, 5, 7, 9]; + let correct_j = vec![0, 0, 1, 0, 2, 0, 3, 0, 4]; + let correct_x = vec![9.0, 1.5, 0.5, 6.0, 12.0, 0.75, 0.625, 3.0, 16.0]; + assert_eq!(&csr.row_pointers, &correct_p); + assert_eq!(&csr.col_indices, &correct_j); + vec_approx_eq(&csr.values, &correct_x, 1e-15); + } + + #[test] + fn csr_matrix_lower_triangular_with_diagonal_first() { + // 9.00 1.5 6.0 0.750 3.0 + // 1.50 0.5 0.0 0.000 0.0 + // 6.00 0.0 12.0 0.000 0.0 + // 0.75 0.0 0.0 0.625 0.0 + // 3.00 0.0 0.0 0.000 16.0 + // lower triangular with diagonal entries being set first + let mut coo = CooMatrix::new(Layout::Lower, 5, 5, 9).unwrap(); + // diagonal + coo.put(0, 0, 9.0).unwrap(); + coo.put(1, 1, 0.5).unwrap(); + coo.put(2, 2, 12.0).unwrap(); + coo.put(3, 3, 0.625).unwrap(); + coo.put(4, 4, 16.0).unwrap(); + // lower diagonal + coo.put(1, 0, 1.5).unwrap(); + coo.put(2, 0, 6.0).unwrap(); + coo.put(3, 0, 0.75).unwrap(); + coo.put(4, 0, 3.0).unwrap(); + let csr = CsrMatrix::from(&coo); + // solution + let correct_p = vec![0, 1, 3, 5, 7, 9]; + let correct_j = vec![0, 0, 1, 0, 2, 0, 3, 0, 4]; + let correct_x = vec![9.0, 1.5, 0.5, 6.0, 12.0, 0.75, 0.625, 3.0, 16.0]; + assert_eq!(&csr.row_pointers, &correct_p); + assert_eq!(&csr.col_indices, &correct_j); + vec_approx_eq(&csr.values, &correct_x, 1e-15); + } + + #[test] + fn to_matrix_fails_on_wrong_dims() { + let csr = CsrMatrix::new(Layout::Upper, 2, 2, 3); + let mut a_2x1 = Matrix::new(3, 1); + let mut a_1x2 = Matrix::new(1, 3); + assert_eq!(csr.to_matrix(&mut a_2x1), Err("wrong matrix dimensions")); + assert_eq!(csr.to_matrix(&mut a_1x2), Err("wrong matrix dimensions")); + } + + #[test] + fn to_matrix_and_as_matrix_work() { + let csr = CsrMatrix { + layout: Layout::Full, + nrow: 5, + ncol: 5, + row_pointers: vec![0, 2, 5, 8, 9, 12], + col_indices: vec![ + // p + 0, 1, // i = 0, count = 0, 1 + 0, 2, 4, // i = 1, count = 2, 3, 4 + 1, 2, 3, // i = 2, count = 5, 6, 7 + 2, // i = 3, count = 8 + 1, 2, 4, // i = 4, count = 9, 10, 11 + // count = 12 + ], + values: vec![ + // p + 2.0, 3.0, // i = 0, count = 0, 1 + 3.0, 4.0, 6.0, // i = 1, count = 2, 3, 4 + -1.0, -3.0, 2.0, // i = 2, count = 5, 6, 7 + 1.0, // i = 3, count = 8 + 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 + // count = 12 + ], + }; + + // covert to dense + let mut a = Matrix::new(5, 5); + csr.to_matrix(&mut a).unwrap(); + let correct = "┌ ┐\n\ + │ 2 3 0 0 0 │\n\ + │ 3 0 4 0 6 │\n\ + │ 0 -1 -3 2 0 │\n\ + │ 0 0 1 0 0 │\n\ + │ 0 4 2 0 1 │\n\ + └ ┘"; + assert_eq!(format!("{}", a), correct); + + // use as_matrix + let b = csr.as_matrix().unwrap(); + assert_eq!(format!("{}", b), correct); + } + + #[test] + fn to_matrix_upper_works() { + // upper triangular + let csr = CsrMatrix { + layout: Layout::Upper, + nrow: 5, + ncol: 5, + row_pointers: vec![0, 5, 6, 7, 8, 9], + col_indices: vec![0, 1, 2, 3, 4, 1, 2, 3, 4], + values: vec![9.0, 1.5, 6.0, 0.75, 3.0, 0.5, 12.0, 0.625, 16.0], + }; + let a = csr.as_matrix().unwrap(); + let correct = "┌ ┐\n\ + │ 9 1.5 6 0.75 3 │\n\ + │ 1.5 0.5 0 0 0 │\n\ + │ 6 0 12 0 0 │\n\ + │ 0.75 0 0 0.625 0 │\n\ + │ 3 0 0 0 16 │\n\ + └ ┘"; + assert_eq!(format!("{}", a), correct); + } + + #[test] + fn to_matrix_lower_works() { + // upper triangular + let csr = CsrMatrix { + layout: Layout::Lower, + nrow: 5, + ncol: 5, + row_pointers: vec![0, 1, 3, 5, 7, 9], + col_indices: vec![0, 0, 1, 0, 2, 0, 3, 0, 4], + values: vec![9.0, 1.5, 0.5, 6.0, 12.0, 0.75, 0.625, 3.0, 16.0], + }; + let a = csr.as_matrix().unwrap(); + let correct = "┌ ┐\n\ + │ 9 1.5 6 0.75 3 │\n\ + │ 1.5 0.5 0 0 0 │\n\ + │ 6 0 12 0 0 │\n\ + │ 0.75 0 0 0.625 0 │\n\ + │ 3 0 0 0 16 │\n\ + └ ┘"; + assert_eq!(format!("{}", a), correct); + } +} diff --git a/russell_sparse/src/enums.rs b/russell_sparse/src/enums.rs index a268cb56..5f415e98 100644 --- a/russell_sparse/src/enums.rs +++ b/russell_sparse/src/enums.rs @@ -1,5 +1,48 @@ use crate::StrError; +/// Defines how the CooMatrix (triplet) represents a matrix +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub enum Layout { + /// Lower triangular (e.g., for MUMPS) + Lower, + + /// Upper triangular (e.g., for Intel DSS) + Upper, + + /// Full matrix (e.g., for UMFPACK) + Full, +} + +/// 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. +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub enum SymmetricHandling { + /// Leave the layout as lower triangular (if symmetric) + /// + /// **Note:** Lower triangular is the standard MatrixMarket format. + /// Thus, this option will do nothing. + /// + /// This option is useful for the MUMPS solver. + LeaveAsLower, + + /// Convert the layout to upper triangular (if symmetric) + /// + /// **Note:** Since lower triangular is standard in MatrixMarket, + /// this option will swap the lower triangle to the upper triangle. + /// + /// This option is useful for the Intel DSS solver. + SwapToUpper, + + /// Make the matrix full (if symmetric) + /// + /// **Note:: Mirror the lower triangle to the upper triangle (duplicate data). + /// The number of non-zeros will be slightly larger than just duplicating the lower triangle. + /// + /// This option is useful for the UMFPACK solver. + MakeItFull, +} + /// Matrix symmetry option #[derive(Clone, Copy, Debug)] pub enum Symmetry { @@ -229,11 +272,26 @@ pub(crate) fn str_umf_scaling(umf_code: i32) -> &'static str { mod tests { use super::{ code_symmetry_mmp, code_symmetry_umf, enum_ordering, enum_scaling, str_enum_ordering, str_enum_scaling, - str_mmp_ordering, str_mmp_scaling, str_umf_ordering, str_umf_scaling, LinSolKind, Ordering, Scaling, Symmetry, + str_mmp_ordering, str_mmp_scaling, str_umf_ordering, str_umf_scaling, Layout, LinSolKind, Ordering, Scaling, + SymmetricHandling, Symmetry, }; #[test] fn clone_copy_and_debug_work() { + let layout = Layout::Full; + let copy = layout; + let clone = layout.clone(); + assert_eq!(format!("{:?}", layout), "Full"); + assert_eq!(copy, Layout::Full); + assert_eq!(clone, Layout::Full); + + let handling = SymmetricHandling::LeaveAsLower; + let copy = handling; + let clone = handling.clone(); + assert_eq!(format!("{:?}", handling), "LeaveAsLower"); + assert_eq!(copy, SymmetricHandling::LeaveAsLower); + assert_eq!(clone, SymmetricHandling::LeaveAsLower); + let symmetry = Symmetry::General; let copy = symmetry; let clone = symmetry.clone(); diff --git a/russell_sparse/src/lib.rs b/russell_sparse/src/lib.rs index b2f88b1a..cdce7cda 100644 --- a/russell_sparse/src/lib.rs +++ b/russell_sparse/src/lib.rs @@ -6,30 +6,30 @@ //! //! ``` //! 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, nnz) = (5, 13); -//! 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\ @@ -40,13 +40,13 @@ //! 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\ @@ -64,18 +64,22 @@ pub type StrError = &'static str; mod config_solver; +mod coo_matrix; +mod csr_matrix; mod enums; pub mod prelude; mod read_matrix_market; mod solver; -mod sparse_triplet; mod verify_lin_sys; +mod write_matrix_market; pub use crate::config_solver::*; +pub use crate::coo_matrix::*; +pub use crate::csr_matrix::*; pub use crate::enums::*; pub use crate::read_matrix_market::*; pub use crate::solver::*; -pub use crate::sparse_triplet::*; pub use crate::verify_lin_sys::*; +pub use crate::write_matrix_market::*; // run code from README file #[cfg(doctest)] diff --git a/russell_sparse/src/prelude.rs b/russell_sparse/src/prelude.rs index 4bbe7a76..ba9f0f5a 100644 --- a/russell_sparse/src/prelude.rs +++ b/russell_sparse/src/prelude.rs @@ -4,5 +4,9 @@ //! access to commonly used functionality. pub use crate::config_solver::ConfigSolver; +pub use crate::coo_matrix::CooMatrix; +pub use crate::csr_matrix::CsrMatrix; +pub use crate::enums::*; +pub use crate::read_matrix_market; pub use crate::solver::Solver; -pub use crate::sparse_triplet::SparseTriplet; +pub use crate::verify_lin_sys::VerifyLinSys; diff --git a/russell_sparse/src/read_matrix_market.rs b/russell_sparse/src/read_matrix_market.rs index 8e4f1062..ca8f9304 100644 --- a/russell_sparse/src/read_matrix_market.rs +++ b/russell_sparse/src/read_matrix_market.rs @@ -1,7 +1,9 @@ -use super::SparseTriplet; +use super::{CooMatrix, Layout, SymmetricHandling}; use crate::StrError; +use std::ffi::OsStr; use std::fs::File; use std::io::{BufRead, BufReader}; +use std::path::Path; struct MatrixMarketData { // header @@ -34,7 +36,7 @@ impl MatrixMarketData { } #[inline] - fn parse_header(&mut self, line: &String) -> Result<(), StrError> { + fn parse_header(&mut self, line: &str) -> Result<(), StrError> { let mut data = line.trim_start().trim_end_matches("\n").split_whitespace(); match data.next() { @@ -86,7 +88,7 @@ impl MatrixMarketData { } #[inline] - fn parse_dimensions(&mut self, line: &String) -> Result { + fn parse_dimensions(&mut self, line: &str) -> Result { let maybe_data = line.trim_start().trim_end_matches("\n"); if maybe_data.starts_with("%") || maybe_data == "" { return Ok(false); // ignore comments or empty lines; returns false == not parsed @@ -118,7 +120,7 @@ impl MatrixMarketData { } #[inline] - fn parse_triple(&mut self, line: &String) -> Result { + fn parse_triple(&mut self, line: &str) -> Result { let maybe_data = line.trim_start().trim_end_matches("\n"); if maybe_data.starts_with("%") || maybe_data == "" { return Ok(false); // ignore comments or empty lines @@ -146,7 +148,7 @@ impl MatrixMarketData { None => return Err("cannot read value aij"), }; - self.i -= 1; // MatrixMarket is one-based + self.i -= 1; // MatrixMarket is one-based, so make it zero-based here self.j -= 1; if self.i < 0 || self.i >= self.m || self.j < 0 || self.j >= self.n { @@ -159,42 +161,31 @@ impl MatrixMarketData { } } -/// Reads a MatrixMarket file into a SparseTriplet -/// -/// **Note:** This function works only with square matrices. +/// Reads a MatrixMarket file into a CooMatrix /// /// # Input /// -/// * `filepath` -- The full file path with filename -/// * `sym_mirror` -- Tells the reader to mirror the **off diagonal** entries, -/// if the symmetric option is found in the header. -/// -/// ## Remarks on sym_mirror +/// * `full_path` -- may be a String, &str, or Path +/// * `symmetric_handling` -- Options to handle symmetric matrices /// -/// ```text -/// if i != j, read line and set a(i,j) = a(j,i) = line data -/// ``` +/// # Output /// -/// If the matrix is symmetric, only entries in the **lower triangular** portion -/// are present in the MatrixMarket file (see reference). However, some solvers -/// (e.g., UMFPACK) require the complete sparse dataset (both off-diagonals), -/// even if the matrix is symmetric. Other solvers (e.g. Mu-M-P-S) must **not** -/// receive both off-diagonal sides when working with symmetric matrices. -/// Therefore, the user has to decide when to use the `sym_mirror` flag. +/// * `coo` -- the CooMatrix +/// * `symmetric` -- whether MatrixMarket is flagged as symmetric or not /// -/// If `sym_mirror` is true, the reader will set `nnz` (number of non-zero values) -/// with twice the specified `nnz` value because we cannot know how many entries -/// are on the diagonal until the whole file is read. Nonetheless, the `SparseTriplet` -/// can be used normally by the user, since this information is internal to `SparseTriplet`. +/// ## Remarks on symmetric matrices /// -/// # Output +/// If the matrix is symmetric, only entries in the **lower triangular** portion +/// are present in the MatrixMarket file (see reference). Thus, the `symmetric_handling` +/// may be used to: /// -/// * A SparseTriplet or an error message -/// * Returns true if the `symmetric` keyword is present in the header +/// 1. Leave the data as it is, i.e., return a Lower Triangular CooMatrix (for MUMPS solver) +/// 2. Swap the lower triangle with the upper triangle, i.e., return an Upper Triangular CooMatrix (for Intel DSS solver) +/// 3. Duplicate the data to make a full matrix, i.e., return a Full CooMatrix (for UMFPACK solver) /// /// # Panics /// -/// This function may panic but should not panic (please contact us if it panics). +/// This function may panic but should not panic (please contact us if it panics :-). /// /// # Example of MatrixMarket file /// @@ -253,7 +244,7 @@ impl MatrixMarketData { /// /// ## Example 1 - General matrix /// -/// Given the following `simple_gen.mtx` file: +/// Given the following `ok_simple_general.mtx` file: /// /// ```text /// %%MatrixMarket matrix coordinate real general @@ -269,20 +260,20 @@ impl MatrixMarketData { /// /// ``` /// use russell_lab::Matrix; -/// use russell_sparse::{read_matrix_market, StrError}; +/// use russell_sparse::{read_matrix_market, Layout, SymmetricHandling, StrError}; /// /// fn main() -> Result<(), StrError> { -/// let filepath = "./data/matrix_market/simple_gen.mtx".to_string(); -/// let (trip, symmetric) = read_matrix_market(&filepath, false)?; -/// let neq = trip.neq(); -/// let mut a = Matrix::new(neq, neq); -/// trip.to_matrix(&mut a)?; +/// let filepath = "./data/matrix_market/ok_simple_general.mtx".to_string(); +/// let (coo, sym) = read_matrix_market(&filepath, SymmetricHandling::LeaveAsLower)?; +/// let mut a = Matrix::new(coo.nrow, coo.ncol); +/// coo.to_matrix(&mut a)?; /// let correct = "┌ ┐\n\ /// │ 1 2 0 │\n\ /// │ 3 4 0 │\n\ /// │ 0 0 5 │\n\ /// └ ┘"; -/// assert_eq!(symmetric, false); +/// assert!(!sym); +/// assert_eq!(coo.layout, Layout::Full); /// assert_eq!(format!("{}", a), correct); /// Ok(()) /// } @@ -290,7 +281,7 @@ impl MatrixMarketData { /// /// ## Example 2 - Symmetric matrix /// -/// Given the following `simple_sym.mtx` file: +/// Given the following `ok_simple_symmetric.mtx` file: /// /// ```text /// %%MatrixMarket matrix coordinate real symmetric @@ -305,26 +296,33 @@ impl MatrixMarketData { /// /// ``` /// use russell_lab::Matrix; -/// use russell_sparse::{read_matrix_market, StrError}; +/// use russell_sparse::{read_matrix_market, Layout, SymmetricHandling, StrError}; /// /// fn main() -> Result<(), StrError> { -/// let filepath = "./data/matrix_market/simple_sym.mtx".to_string(); -/// let (trip, symmetric) = read_matrix_market(&filepath, true)?; -/// let neq = trip.neq(); -/// let mut a = Matrix::new(neq, neq); -/// trip.to_matrix(&mut a)?; +/// let filepath = "./data/matrix_market/ok_simple_symmetric.mtx".to_string(); +/// let (coo, sym) = read_matrix_market(&filepath, SymmetricHandling::LeaveAsLower)?; +/// let mut a = Matrix::new(coo.nrow, coo.ncol); +/// coo.to_matrix(&mut a)?; /// let correct = "┌ ┐\n\ /// │ 1 2 0 │\n\ /// │ 2 3 4 │\n\ /// │ 0 4 0 │\n\ /// └ ┘"; -/// assert_eq!(symmetric, true); +/// assert!(sym); +/// assert_eq!(coo.layout, Layout::Lower); /// assert_eq!(format!("{}", a), correct); /// Ok(()) /// } /// ``` -pub fn read_matrix_market(filepath: &String, sym_mirror: bool) -> Result<(SparseTriplet, bool), StrError> { - let input = File::open(filepath).map_err(|_| "cannot open file")?; +pub fn read_matrix_market

( + full_path: &P, + symmetric_handling: SymmetricHandling, +) -> Result<(CooMatrix, bool), StrError> +where + P: AsRef + ?Sized, +{ + let path = Path::new(full_path).to_path_buf(); + let input = File::open(path).map_err(|_| "cannot open file")?; let buffered = BufReader::new(input); let mut lines_iter = buffered.lines(); @@ -348,19 +346,28 @@ pub fn read_matrix_market(filepath: &String, sym_mirror: bool) -> Result<(Sparse } } - // check dimensions - if data.m != data.n { - return Err("cannot read non-square matrix"); - } + // specify the CooMatrix layout + let layout = 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"); + } + match symmetric_handling { + SymmetricHandling::LeaveAsLower => Layout::Lower, + SymmetricHandling::SwapToUpper => Layout::Upper, + SymmetricHandling::MakeItFull => Layout::Full, + } + } else { + Layout::Full + }; // set max number of entries let mut max = data.nnz; - if data.symmetric && sym_mirror { + if data.symmetric && symmetric_handling == SymmetricHandling::MakeItFull { max = 2 * data.nnz; } // allocate triplet - let mut trip = SparseTriplet::new(data.m as usize, max as usize)?; + let mut coo = CooMatrix::new(layout, data.m as usize, data.n as usize, max as usize).unwrap(); // read and parse triples loop { @@ -368,10 +375,24 @@ pub fn read_matrix_market(filepath: &String, sym_mirror: bool) -> Result<(Sparse Some(v) => { let line = v.unwrap(); // must panic because no error expected here if data.parse_triple(&line)? { - trip.put(data.i as usize, data.j as usize, data.aij)?; - if data.symmetric && sym_mirror && data.i != data.j { - trip.put(data.j as usize, data.i as usize, data.aij)?; - } + if data.symmetric { + match symmetric_handling { + SymmetricHandling::LeaveAsLower => { + coo.put(data.i as usize, data.j as usize, data.aij).unwrap(); + } + SymmetricHandling::SwapToUpper => { + coo.put(data.j as usize, data.i as usize, data.aij).unwrap(); + } + SymmetricHandling::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(); + } + } + } + } else { + coo.put(data.i as usize, data.j as usize, data.aij).unwrap(); + }; } } None => break, @@ -383,7 +404,7 @@ pub fn read_matrix_market(filepath: &String, sym_mirror: bool) -> Result<(Sparse return Err("not all triples (i,j,aij) have been found"); } - Ok((trip, data.symmetric)) + Ok((coo, data.symmetric)) } //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -391,6 +412,7 @@ pub fn read_matrix_market(filepath: &String, sym_mirror: bool) -> Result<(Sparse #[cfg(test)] mod tests { use super::{read_matrix_market, MatrixMarketData}; + use crate::{Layout, SymmetricHandling}; use russell_lab::Matrix; #[test] @@ -398,47 +420,47 @@ mod tests { let mut data = MatrixMarketData::new(); assert_eq!( - data.parse_header(&String::from(" \n")), + data.parse_header(" \n"), Err("cannot find the keyword %%MatrixMarket on the first line") ); assert_eq!( - data.parse_header(&String::from("MatrixMarket ")), + data.parse_header("MatrixMarket "), Err("the header (first line) must start with %%MatrixMarket"), ); assert_eq!( - data.parse_header(&String::from(" %%MatrixMarket")), + data.parse_header(" %%MatrixMarket"), Err("cannot find the first option in the header line"), ); assert_eq!( - data.parse_header(&String::from("%%MatrixMarket wrong")), + data.parse_header("%%MatrixMarket wrong"), Err("after %%MatrixMarket, the first option must be \"matrix\""), ); assert_eq!( - data.parse_header(&String::from("%%MatrixMarket matrix ")), + data.parse_header("%%MatrixMarket matrix "), Err("cannot find the second option in the header line"), ); assert_eq!( - data.parse_header(&String::from("%%MatrixMarket matrix wrong")), + data.parse_header("%%MatrixMarket matrix wrong"), Err("after %%MatrixMarket, the second option must be \"coordinate\""), ); assert_eq!( - data.parse_header(&String::from("%%MatrixMarket matrix coordinate")), + data.parse_header("%%MatrixMarket matrix coordinate"), Err("cannot find the third option in the header line"), ); assert_eq!( - data.parse_header(&String::from("%%MatrixMarket matrix coordinate wrong")), + data.parse_header("%%MatrixMarket matrix coordinate wrong"), Err("after %%MatrixMarket, the third option must be \"real\""), ); assert_eq!( - data.parse_header(&String::from("%%MatrixMarket matrix coordinate real")), + data.parse_header("%%MatrixMarket matrix coordinate real"), Err("cannot find the fourth option in the header line"), ); assert_eq!( - data.parse_header(&String::from(" %%MatrixMarket matrix coordinate real wrong")), + data.parse_header(" %%MatrixMarket matrix coordinate real wrong"), Err("after %%MatrixMarket, the fourth option must be either \"general\" or \"symmetric\""), ); } @@ -448,38 +470,38 @@ mod tests { let mut data = MatrixMarketData::new(); assert_eq!( - data.parse_dimensions(&String::from(" wrong \n")).err(), + data.parse_dimensions(" wrong \n").err(), Some("cannot parse number of rows") ); assert_eq!( - data.parse_dimensions(&String::from(" 1 \n")).err(), + data.parse_dimensions(" 1 \n").err(), Some("cannot read number of columns") ); assert_eq!( - data.parse_dimensions(&String::from(" 1 wrong")).err(), + data.parse_dimensions(" 1 wrong").err(), Some("cannot parse number of columns") ); assert_eq!( - data.parse_dimensions(&String::from(" 1 1 \n")).err(), + data.parse_dimensions(" 1 1 \n").err(), Some("cannot read number of non-zeros") ); assert_eq!( - data.parse_dimensions(&String::from(" 1 1 wrong")).err(), + data.parse_dimensions(" 1 1 wrong").err(), Some("cannot parse number of non-zeros") ); assert_eq!( - data.parse_dimensions(&String::from(" 0 1 1")).err(), + data.parse_dimensions(" 0 1 1").err(), Some("found invalid (zero or negative) dimensions") ); assert_eq!( - data.parse_dimensions(&String::from(" 1 0 1")).err(), + data.parse_dimensions(" 1 0 1").err(), Some("found invalid (zero or negative) dimensions") ); assert_eq!( - data.parse_dimensions(&String::from(" 1 1 0")).err(), + data.parse_dimensions(" 1 1 0").err(), Some("found invalid (zero or negative) dimensions") ); } @@ -491,121 +513,114 @@ mod tests { data.n = 2; data.nnz = 1; - assert_eq!( - data.parse_triple(&String::from(" wrong \n")).err(), - Some("cannot parse index i") - ); + assert_eq!(data.parse_triple(" wrong \n").err(), Some("cannot parse index i")); - assert_eq!( - data.parse_triple(&String::from(" 1 \n")).err(), - Some("cannot read index j") - ); - assert_eq!( - data.parse_triple(&String::from(" 1 wrong")).err(), - Some("cannot parse index j") - ); + assert_eq!(data.parse_triple(" 1 \n").err(), Some("cannot read index j")); + assert_eq!(data.parse_triple(" 1 wrong").err(), Some("cannot parse index j")); - assert_eq!( - data.parse_triple(&String::from(" 1 1 \n")).err(), - Some("cannot read value aij") - ); - assert_eq!( - data.parse_triple(&String::from(" 1 1 wrong")).err(), - Some("cannot parse value aij") - ); + assert_eq!(data.parse_triple(" 1 1 \n").err(), Some("cannot read value aij")); + assert_eq!(data.parse_triple(" 1 1 wrong").err(), Some("cannot parse value aij")); - assert_eq!( - data.parse_triple(&String::from(" 0 1 1")).err(), - Some("found invalid indices") - ); - assert_eq!( - data.parse_triple(&String::from(" 3 1 1")).err(), - Some("found invalid indices") - ); - assert_eq!( - data.parse_triple(&String::from(" 1 0 1")).err(), - Some("found invalid indices") - ); - assert_eq!( - data.parse_triple(&String::from(" 1 3 1")).err(), - Some("found invalid indices") - ); + assert_eq!(data.parse_triple(" 0 1 1").err(), Some("found invalid indices")); + assert_eq!(data.parse_triple(" 3 1 1").err(), Some("found invalid indices")); + assert_eq!(data.parse_triple(" 1 0 1").err(), Some("found invalid indices")); + assert_eq!(data.parse_triple(" 1 3 1").err(), Some("found invalid indices")); } #[test] fn read_matrix_market_handle_wrong_files() { + let h = SymmetricHandling::LeaveAsLower; + assert_eq!(read_matrix_market("__wrong__", h).err(), Some("cannot open file")); assert_eq!( - read_matrix_market(&String::from("__wrong__"), false).err(), - Some("cannot open file") - ); - assert_eq!( - read_matrix_market(&String::from("./data/matrix_market/bad_empty_file.mtx"), false).err(), + read_matrix_market("./data/matrix_market/bad_empty_file.mtx", h).err(), Some("file is empty") ); assert_eq!( - read_matrix_market(&String::from("./data/matrix_market/bad_wrong_header.mtx"), false).err(), + read_matrix_market("./data/matrix_market/bad_wrong_header.mtx", h).err(), Some("after %%MatrixMarket, the first option must be \"matrix\"") ); assert_eq!( - read_matrix_market(&String::from("./data/matrix_market/bad_wrong_dims.mtx"), false).err(), + read_matrix_market("./data/matrix_market/bad_wrong_dims.mtx", h).err(), Some("found invalid (zero or negative) dimensions") ); assert_eq!( - read_matrix_market(&String::from("./data/matrix_market/bad_rectangular.mtx"), false).err(), - Some("cannot read non-square matrix") - ); - assert_eq!( - read_matrix_market(&String::from("./data/matrix_market/bad_missing_data.mtx"), false).err(), + read_matrix_market("./data/matrix_market/bad_missing_data.mtx", h).err(), Some("not all triples (i,j,aij) have been found") ); assert_eq!( - read_matrix_market(&String::from("./data/matrix_market/bad_many_lines.mtx"), false).err(), + read_matrix_market("./data/matrix_market/bad_many_lines.mtx", h).err(), Some("there are more (i,j,aij) triples than specified") ); + 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") + ); } #[test] fn read_matrix_market_works() { - let filepath = "./data/matrix_market/ok1.mtx".to_string(); - let (trip, sym) = read_matrix_market(&filepath, false).unwrap(); - assert_eq!(sym, false); - assert_eq!((trip.neq, trip.pos, trip.max), (5, 12, 12)); - assert_eq!(trip.indices_i, &[0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4]); - assert_eq!(trip.indices_j, &[0, 0, 1, 1, 1, 2, 2, 2, 2, 3, 4, 4]); - assert_eq!( - trip.values_aij, + let h = SymmetricHandling::LeaveAsLower; + let filepath = "./data/matrix_market/ok_general.mtx".to_string(); + let (coo, sym) = read_matrix_market(&filepath, h).unwrap(); + assert!(!sym); + assert_eq!(coo.layout, Layout::Full); + assert_eq!((coo.nrow, coo.ncol, coo.pos, coo.max), (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]); + assert_eq!( + coo.values_aij, &[2.0, 3.0, 3.0, -1.0, 4.0, 4.0, -3.0, 1.0, 2.0, 2.0, 6.0, 1.0] ); } #[test] - fn read_matrix_market_sym_triangle_works() { - let filepath = "./data/matrix_market/ok2.mtx".to_string(); - let (trip, sym) = read_matrix_market(&filepath, false).unwrap(); - assert_eq!(sym, true); - assert_eq!((trip.neq, trip.pos, trip.max), (5, 15, 15)); - assert_eq!(trip.indices_i, &[0, 1, 2, 3, 4, 0, 0, 0, 0, 1, 1, 1, 2, 2, 3]); - assert_eq!(trip.indices_j, &[0, 1, 2, 3, 4, 1, 2, 3, 4, 2, 3, 4, 3, 4, 4]); - assert_eq!( - trip.values_aij, + fn read_matrix_market_symmetric_lower_works() { + let h = SymmetricHandling::LeaveAsLower; + let filepath = "./data/matrix_market/ok_symmetric.mtx".to_string(); + let (coo, sym) = read_matrix_market(&filepath, h).unwrap(); + assert!(sym); + assert_eq!(coo.layout, Layout::Lower); + assert_eq!((coo.nrow, coo.ncol, coo.pos, coo.max), (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]); + assert_eq!( + coo.values_aij, &[2.0, 2.0, 9.0, 7.0, 8.0, 1.0, 1.0, 3.0, 2.0, 2.0, 1.0, 1.0, 1.0, 5.0, 1.0], ); } #[test] - fn read_matrix_market_sym_mirror_works() { - let filepath = "./data/matrix_market/ok3.mtx".to_string(); - let (trip, sym) = read_matrix_market(&filepath, true).unwrap(); - assert_eq!(sym, true); - assert_eq!((trip.neq, trip.pos, trip.max), (5, 11, 14)); - assert_eq!(trip.indices_i, &[0, 1, 0, 2, 1, 3, 2, 3, 4, 1, 4, 0, 0, 0]); - assert_eq!(trip.indices_j, &[0, 0, 1, 1, 2, 2, 3, 3, 1, 4, 4, 0, 0, 0]); - assert_eq!( - trip.values_aij, + fn read_matrix_market_symmetric_upper_works() { + let h = SymmetricHandling::SwapToUpper; + let filepath = "./data/matrix_market/ok_symmetric.mtx".to_string(); + let (coo, sym) = read_matrix_market(&filepath, h).unwrap(); + assert!(sym); + assert_eq!(coo.layout, Layout::Upper); + assert_eq!((coo.nrow, coo.ncol, coo.pos, coo.max), (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]); + assert_eq!( + coo.values_aij, + &[2.0, 2.0, 9.0, 7.0, 8.0, 1.0, 1.0, 3.0, 2.0, 2.0, 1.0, 1.0, 1.0, 5.0, 1.0], + ); + } + + #[test] + fn read_matrix_market_symmetric_to_full_works() { + let h = SymmetricHandling::MakeItFull; + let filepath = "./data/matrix_market/ok_symmetric_small.mtx".to_string(); + let (coo, sym) = read_matrix_market(&filepath, h).unwrap(); + assert!(sym); + assert_eq!(coo.layout, Layout::Full); + assert_eq!((coo.nrow, coo.ncol, coo.pos, coo.max), (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]); + assert_eq!( + coo.values_aij, &[2.0, 3.0, 3.0, -1.0, -1.0, 2.0, 2.0, 3.0, 6.0, 6.0, 1.0, 0.0, 0.0, 0.0] ); let mut a = Matrix::new(5, 5); - trip.to_matrix(&mut a).unwrap(); + coo.to_matrix(&mut a).unwrap(); let correct = "┌ ┐\n\ │ 2 3 0 0 0 │\n\ │ 3 0 -1 0 6 │\n\ diff --git a/russell_sparse/src/solver.rs b/russell_sparse/src/solver.rs index 8cf9b37d..e2cd58d9 100644 --- a/russell_sparse/src/solver.rs +++ b/russell_sparse/src/solver.rs @@ -1,11 +1,10 @@ use super::{ code_symmetry_mmp, code_symmetry_umf, str_enum_ordering, str_enum_scaling, str_mmp_ordering, str_mmp_scaling, - str_umf_ordering, str_umf_scaling, ConfigSolver, LinSolKind, SparseTriplet, + str_umf_ordering, str_umf_scaling, ConfigSolver, CooMatrix, LinSolKind, }; use crate::{StrError, Symmetry}; -use russell_lab::{format_nanoseconds, vec_copy, Stopwatch, Vector}; +use russell_lab::{vec_copy, Stopwatch, Vector}; use russell_openblas::to_i32; -use std::fmt; #[repr(C)] pub(crate) struct ExtSolver { @@ -76,7 +75,7 @@ pub struct Solver { kind: LinSolKind, // solver kind verbose: i32, // verbose mode done_factorize: bool, // factorization completed - neq: usize, // number of equations == nrow(a) where a*x=rhs + nrow: usize, // number of equations == nrow(a) where a*x=rhs solver: *mut ExtSolver, // data allocated by the c-code stopwatch: Stopwatch, // stopwatch to measure elapsed time time_fact: u128, // elapsed time during factorize @@ -87,8 +86,8 @@ pub struct Solver { impl Solver { /// Creates a new solver - pub fn new(config: ConfigSolver, neq: usize, nnz: usize, symmetry: Option) -> Result { - let n = to_i32(neq); + pub fn new(config: ConfigSolver, nrow: usize, nnz: usize, symmetry: Option) -> Result { + let n = to_i32(nrow); let nnz = to_i32(nnz); unsafe { let solver = match config.lin_sol_kind { @@ -136,7 +135,7 @@ impl Solver { kind: config.lin_sol_kind, verbose: config.verbose, done_factorize: false, - neq, + nrow, solver, stopwatch: Stopwatch::new(""), time_fact: 0, @@ -148,9 +147,12 @@ impl Solver { } /// Performs the factorization - pub fn factorize(&mut self, trip: &SparseTriplet) -> Result<(), StrError> { - if trip.neq != self.neq { - return Err("cannot factorize because the triplet has incompatible number of equations"); + pub fn factorize(&mut self, coo: &CooMatrix) -> Result<(), StrError> { + if coo.nrow != coo.ncol { + return Err("cannot factorize because the CooMatrix is not square"); + } + if coo.nrow != self.nrow { + return Err("cannot factorize because the CooMatrix has incompatible number of rows"); } self.stopwatch.reset(); unsafe { @@ -158,9 +160,9 @@ impl Solver { LinSolKind::Mmp => { let res = solver_mmp_factorize( self.solver, - trip.indices_i.as_ptr(), - trip.indices_j.as_ptr(), - trip.values_aij.as_ptr(), + coo.indices_i.as_ptr(), + coo.indices_j.as_ptr(), + coo.values_aij.as_ptr(), self.verbose, ); if res != 0 { @@ -174,9 +176,9 @@ impl Solver { LinSolKind::Umf => { let res = solver_umf_factorize( self.solver, - trip.indices_i.as_ptr(), - trip.indices_j.as_ptr(), - trip.values_aij.as_ptr(), + coo.indices_i.as_ptr(), + coo.indices_j.as_ptr(), + coo.values_aij.as_ptr(), self.verbose, ); if res != 0 { @@ -200,29 +202,29 @@ impl Solver { /// /// ``` /// use russell_lab::{Matrix, Vector}; - /// use russell_sparse::{ConfigSolver, SparseTriplet, Solver, StrError}; + /// use russell_sparse::{ConfigSolver, CooMatrix, Layout, Solver, StrError}; /// /// fn main() -> Result<(), StrError> { /// // allocate a square matrix - /// let (neq, nnz) = (5, 13); - /// 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, nrow); + /// coo.to_matrix(&mut a)?; /// let correct = "┌ ┐\n\ /// │ 2 3 0 0 0 │\n\ /// │ 3 0 4 0 6 │\n\ @@ -233,13 +235,13 @@ impl Solver { /// 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\ @@ -256,7 +258,7 @@ impl Solver { if !self.done_factorize { return Err("factorization must be done before calling solve"); } - if x.dim() != self.neq || rhs.dim() != self.neq { + if x.dim() != self.nrow || rhs.dim() != self.nrow { return Err("x.ndim() and rhs.ndim() must equal the number of equations"); } self.stopwatch.reset(); @@ -309,21 +311,21 @@ impl Solver { /// /// ``` /// use russell_lab::{Matrix, Vector}; - /// use russell_sparse::{ConfigSolver, Solver, SparseTriplet, StrError}; + /// use russell_sparse::{ConfigSolver, CooMatrix, Layout, Solver, StrError}; /// /// fn main() -> Result<(), StrError> { /// // allocate a square matrix - /// let (neq, nnz) = (3, 5); - /// 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 (nrow, ncol, nnz) = (3, 3, 5); + /// 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, nrow); + /// coo.to_matrix(&mut a)?; /// let correct = "┌ ┐\n\ /// │ 0.2 0.2 0 │\n\ /// │ 0.5 -0.25 0 │\n\ @@ -337,7 +339,7 @@ impl Solver { /// /// // 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\ @@ -346,7 +348,7 @@ impl Solver { /// 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\ @@ -357,10 +359,13 @@ impl Solver { /// Ok(()) /// } /// ``` - pub fn compute(config: ConfigSolver, trip: &SparseTriplet, rhs: &Vector) -> Result<(Self, Vector), StrError> { - let mut solver = Solver::new(config, trip.neq, trip.pos, None)?; - let mut x = Vector::new(trip.neq()); - solver.factorize(&trip)?; + pub fn compute(config: ConfigSolver, coo: &CooMatrix, rhs: &Vector) -> Result<(Self, Vector), StrError> { + if coo.nrow != coo.ncol { + return Err("CooMatrix must be symmetric"); + } + let mut solver = Solver::new(config, coo.nrow, coo.pos, None)?; + let mut x = Vector::new(coo.nrow); + solver.factorize(&coo)?; solver.solve(&mut x, &rhs)?; Ok((solver, x)) } @@ -374,6 +379,16 @@ impl Solver { (self.time_fact, self.time_solve) } + /// Returns the ordering effectively used by the solver + pub fn get_effective_ordering(&self) -> String { + self.used_ordering.to_string() + } + + /// Returns the scaling effectively used by the solver + pub fn get_effective_scaling(&self) -> String { + self.used_scaling.to_string() + } + /// Handles error code fn handle_mmp_error_code(err: i32) -> StrError { match err { @@ -492,96 +507,67 @@ impl Drop for Solver { } } -impl fmt::Display for Solver { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - let time_total = self.time_fact + self.time_solve; - write!( - f, - "\x20\x20\x20\x20\"usedOrdering\": \"{}\",\n\ - \x20\x20\x20\x20\"usedScaling\": \"{}\",\n\ - \x20\x20\x20\x20\"doneFactorize\": {},\n\ - \x20\x20\x20\x20\"neq\": {},\n\ - \x20\x20\x20\x20\"timeFactNs\": {},\n\ - \x20\x20\x20\x20\"timeSolveNs\": {},\n\ - \x20\x20\x20\x20\"timeTotalNs\": {},\n\ - \x20\x20\x20\x20\"timeFactStr\": \"{}\",\n\ - \x20\x20\x20\x20\"timeSolveStr\": \"{}\",\n\ - \x20\x20\x20\x20\"timeTotalStr\": \"{}\"", - self.used_ordering, - self.used_scaling, - self.done_factorize, - self.neq, - self.time_fact, - self.time_solve, - time_total, - format_nanoseconds(self.time_fact), - format_nanoseconds(self.time_solve), - format_nanoseconds(time_total) - )?; - Ok(()) - } -} - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #[cfg(test)] mod tests { - use super::{ConfigSolver, LinSolKind, Solver, SparseTriplet}; + use super::{ConfigSolver, CooMatrix, LinSolKind, Solver}; + use crate::Layout; use russell_chk::vec_approx_eq; use russell_lab::Vector; #[test] fn new_works() { let config = ConfigSolver::new(); - let (neq, nnz) = (2, 2); - let solver = Solver::new(config, neq, nnz, None).unwrap(); + let (nrow, nnz) = (2, 2); + let solver = Solver::new(config, nrow, nnz, None).unwrap(); assert_eq!(solver.done_factorize, false); - assert_eq!(solver.neq, 2); + assert_eq!(solver.nrow, 2); } #[test] fn factorize_fails_on_incompatible_triplet() { let config = ConfigSolver::new(); let mut solver = Solver::new(config, 1, 1, None).unwrap(); - let trip = SparseTriplet::new(2, 2).unwrap(); + let coo = CooMatrix::new(Layout::Full, 2, 2, 2).unwrap(); assert_eq!( - solver.factorize(&trip).err(), - Some("cannot factorize because the triplet has incompatible number of equations") + solver.factorize(&coo).err(), + Some("cannot factorize because the CooMatrix has incompatible number of rows") ); } #[test] fn factorize_fails_on_singular_matrix() { let config = ConfigSolver::new(); - let (neq, nnz) = (2, 2); - let mut solver = Solver::new(config, neq, nnz, None).unwrap(); - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 1.0).unwrap(); - trip.put(1, 1, 0.0).unwrap(); - assert_eq!(solver.factorize(&trip), Err("Error(1): Matrix is singular")); + let (nrow, ncol, nnz) = (2, 2, 2); + let mut solver = Solver::new(config, nrow, nnz, None).unwrap(); + let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(1, 1, 0.0).unwrap(); + assert_eq!(solver.factorize(&coo), Err("Error(1): Matrix is singular")); } #[test] fn factorize_works() { let config = ConfigSolver::new(); - let (neq, nnz) = (2, 2); - let mut solver = Solver::new(config, neq, nnz, None).unwrap(); - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 1.0).unwrap(); - trip.put(1, 1, 1.0).unwrap(); - solver.factorize(&trip).unwrap(); + let (nrow, ncol, nnz) = (2, 2, 2); + let mut solver = Solver::new(config, nrow, nnz, None).unwrap(); + let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(1, 1, 1.0).unwrap(); + solver.factorize(&coo).unwrap(); assert!(solver.done_factorize); } #[test] fn solve_fails_on_non_factorized() { let config = ConfigSolver::new(); - let (neq, nnz) = (2, 2); - let mut solver = Solver::new(config, neq, nnz, None).unwrap(); - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 1.0).unwrap(); - trip.put(1, 1, 1.0).unwrap(); - let mut x = Vector::new(neq); + let (nrow, ncol, nnz) = (2, 2, 2); + let mut solver = Solver::new(config, nrow, nnz, None).unwrap(); + let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(1, 1, 1.0).unwrap(); + let mut x = Vector::new(nrow); let rhs = Vector::from(&[1.0, 1.0]); assert_eq!( solver.solve(&mut x, &rhs), @@ -592,12 +578,12 @@ mod tests { #[test] fn solve_fails_on_wrong_vectors() { let config = ConfigSolver::new(); - let (neq, nnz) = (2, 2); - let mut solver = Solver::new(config, neq, nnz, None).unwrap(); - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 1.0).unwrap(); - trip.put(1, 1, 1.0).unwrap(); - solver.factorize(&trip).unwrap(); + let (nrow, ncol, nnz) = (2, 2, 2); + let mut solver = Solver::new(config, nrow, nnz, None).unwrap(); + let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(1, 1, 1.0).unwrap(); + solver.factorize(&coo).unwrap(); let mut x = Vector::new(2); let rhs = Vector::from(&[1.0, 1.0]); let mut x_wrong = Vector::new(1); @@ -615,24 +601,24 @@ mod tests { #[test] fn solve_works() { let config = ConfigSolver::new(); - let (neq, nnz) = (5, 13); - let mut solver = Solver::new(config, neq, nnz, None).unwrap(); + let (nrow, ncol, nnz) = (5, 5, 13); + let mut solver = Solver::new(config, nrow, nnz, None).unwrap(); // allocate a square matrix - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 1.0).unwrap(); // << (0, 0, a00/2) - trip.put(0, 0, 1.0).unwrap(); // << (0, 0, a00/2) - trip.put(1, 0, 3.0).unwrap(); - trip.put(0, 1, 3.0).unwrap(); - trip.put(2, 1, -1.0).unwrap(); - trip.put(4, 1, 4.0).unwrap(); - trip.put(1, 2, 4.0).unwrap(); - trip.put(2, 2, -3.0).unwrap(); - trip.put(3, 2, 1.0).unwrap(); - trip.put(4, 2, 2.0).unwrap(); - trip.put(2, 3, 2.0).unwrap(); - trip.put(1, 4, 6.0).unwrap(); - trip.put(4, 4, 1.0).unwrap(); + let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz).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(); + coo.put(0, 1, 3.0).unwrap(); + coo.put(2, 1, -1.0).unwrap(); + coo.put(4, 1, 4.0).unwrap(); + coo.put(1, 2, 4.0).unwrap(); + coo.put(2, 2, -3.0).unwrap(); + coo.put(3, 2, 1.0).unwrap(); + coo.put(4, 2, 2.0).unwrap(); + coo.put(2, 3, 2.0).unwrap(); + coo.put(1, 4, 6.0).unwrap(); + coo.put(4, 4, 1.0).unwrap(); // allocate x and rhs let mut x = Vector::new(5); @@ -640,7 +626,7 @@ mod tests { let x_correct = &[1.0, 2.0, 3.0, 4.0, 5.0]; // initialize, factorize, and solve - solver.factorize(&trip).unwrap(); + solver.factorize(&coo).unwrap(); solver.solve(&mut x, &rhs).unwrap(); // check @@ -654,33 +640,33 @@ mod tests { fn solver_mmp_behaves_as_expected() { // allocate a new solver let mut config = ConfigSolver::new(); - let (neq, nnz) = (5, 13); + let (nrow, ncol, nnz) = (5, 5, 13); config.lin_sol_kind(LinSolKind::Mmp); - let mut solver = Solver::new(config, neq, nnz, None).unwrap(); + let mut solver = Solver::new(config, nrow, nnz, None).unwrap(); - // factorize fails on incompatible triplet - let mut trip_wrong = SparseTriplet::new(1, 1).unwrap(); - trip_wrong.put(0, 0, 1.0).unwrap(); + // factorize fails on incompatible CooMatrix + let mut coo_wrong = CooMatrix::new(Layout::Full, 1, 1, 1).unwrap(); + coo_wrong.put(0, 0, 1.0).unwrap(); assert_eq!( - solver.factorize(&trip_wrong).err(), - Some("cannot factorize because the triplet has incompatible number of equations") + solver.factorize(&coo_wrong).err(), + Some("cannot factorize because the CooMatrix has incompatible number of rows") ); // allocate a square matrix - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 1.0).unwrap(); // << (0, 0, a00/2) - trip.put(0, 0, 1.0).unwrap(); // << (0, 0, a00/2) - trip.put(1, 0, 3.0).unwrap(); - trip.put(0, 1, 3.0).unwrap(); - trip.put(2, 1, -1.0).unwrap(); - trip.put(4, 1, 4.0).unwrap(); - trip.put(1, 2, 4.0).unwrap(); - trip.put(2, 2, -3.0).unwrap(); - trip.put(3, 2, 1.0).unwrap(); - trip.put(4, 2, 2.0).unwrap(); - trip.put(2, 3, 2.0).unwrap(); - trip.put(1, 4, 6.0).unwrap(); - trip.put(4, 4, 1.0).unwrap(); + let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz).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(); + coo.put(0, 1, 3.0).unwrap(); + coo.put(2, 1, -1.0).unwrap(); + coo.put(4, 1, 4.0).unwrap(); + coo.put(1, 2, 4.0).unwrap(); + coo.put(2, 2, -3.0).unwrap(); + coo.put(3, 2, 1.0).unwrap(); + coo.put(4, 2, 2.0).unwrap(); + coo.put(2, 3, 2.0).unwrap(); + coo.put(1, 4, 6.0).unwrap(); + coo.put(4, 4, 1.0).unwrap(); // allocate x and rhs let mut x = Vector::new(5); @@ -694,7 +680,7 @@ mod tests { ); // factorize works - solver.factorize(&trip).unwrap(); + solver.factorize(&coo).unwrap(); assert!(solver.done_factorize); // solve fails on wrong x vector @@ -721,40 +707,40 @@ mod tests { vec_approx_eq(x_again.as_data(), x_correct, 1e-14); // factorize fails on singular matrix - let mut trip_singular = SparseTriplet::new(5, 2).unwrap(); - trip_singular.put(0, 0, 1.0).unwrap(); - trip_singular.put(4, 4, 1.0).unwrap(); + let mut coo_singular = CooMatrix::new(Layout::Full, 5, 5, 2).unwrap(); + coo_singular.put(0, 0, 1.0).unwrap(); + coo_singular.put(4, 4, 1.0).unwrap(); let mut solver = Solver::new(config, 5, 2, None).unwrap(); assert_eq!( - solver.factorize(&trip_singular), + solver.factorize(&coo_singular), Err("Error(-10): numerically singular matrix") ); } #[test] fn compute_works() { - let (neq, nnz) = (3, 6); - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 1.0).unwrap(); - trip.put(0, 1, 1.0).unwrap(); - trip.put(1, 0, 2.0).unwrap(); - trip.put(1, 1, 1.0).unwrap(); - trip.put(1, 2, 1.0).unwrap(); - trip.put(2, 2, 1.0).unwrap(); + let (nrow, ncol, nnz) = (3, 3, 6); + let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(0, 1, 1.0).unwrap(); + coo.put(1, 0, 2.0).unwrap(); + coo.put(1, 1, 1.0).unwrap(); + coo.put(1, 2, 1.0).unwrap(); + coo.put(2, 2, 1.0).unwrap(); let rhs1 = Vector::from(&[1.0, 2.0, 3.0]); let rhs2 = Vector::from(&[2.0, 4.0, 6.0]); let config = ConfigSolver::new(); - let (mut solver, x1) = Solver::compute(config, &trip, &rhs1).unwrap(); + let (mut solver, x1) = Solver::compute(config, &coo, &rhs1).unwrap(); vec_approx_eq(x1.as_data(), &[-2.0, 3.0, 3.0], 1e-15); - let mut x2 = Vector::new(neq); + let mut x2 = Vector::new(nrow); solver.solve(&mut x2, &rhs2).unwrap(); } #[test] fn get_elapsed_times_works() { let config = ConfigSolver::new(); - let (neq, nnz) = (2, 2); - let solver = Solver::new(config, neq, nnz, None).unwrap(); + let (nrow, nnz) = (2, 2); + let solver = Solver::new(config, nrow, nnz, None).unwrap(); let times = solver.get_elapsed_times(); assert_eq!(times, (0, 0)); } @@ -808,22 +794,4 @@ mod tests { ); assert_eq!(Solver::handle_umf_error_code(123), default); } - - #[test] - fn display_trait_works() { - let config = ConfigSolver::new(); - let (neq, nnz) = (2, 2); - let solver = Solver::new(config, neq, nnz, None).unwrap(); - let b: &str = "\x20\x20\x20\x20\"usedOrdering\": \"Auto\",\n\ - \x20\x20\x20\x20\"usedScaling\": \"Auto\",\n\ - \x20\x20\x20\x20\"doneFactorize\": false,\n\ - \x20\x20\x20\x20\"neq\": 2,\n\ - \x20\x20\x20\x20\"timeFactNs\": 0,\n\ - \x20\x20\x20\x20\"timeSolveNs\": 0,\n\ - \x20\x20\x20\x20\"timeTotalNs\": 0,\n\ - \x20\x20\x20\x20\"timeFactStr\": \"0ns\",\n\ - \x20\x20\x20\x20\"timeSolveStr\": \"0ns\",\n\ - \x20\x20\x20\x20\"timeTotalStr\": \"0ns\""; - assert_eq!(format!("{}", solver), b); - } } diff --git a/russell_sparse/src/sparse_triplet.rs b/russell_sparse/src/sparse_triplet.rs deleted file mode 100644 index 73ba81c3..00000000 --- a/russell_sparse/src/sparse_triplet.rs +++ /dev/null @@ -1,652 +0,0 @@ -use crate::StrError; -use russell_lab::{Matrix, Vector}; -use russell_openblas::to_i32; -use std::fmt; - -/// Holds triples (i,j,aij) representing a sparse matrix -/// -/// # Remarks -/// -/// - Only the non-zero values are required -/// - Entries with repeated (i,j) indices are allowed -/// - Repeated (i,j) entries will have the aij values summed when solving a linear system -/// - The repeated (i,j) capability is of great convenience for Finite Element solvers -/// - A maximum number of entries must be decided prior to allocating a new Triplet -/// - The maximum number of entries includes possible entries with repeated indices -/// - See the `to_matrix` method for an example -pub struct SparseTriplet { - pub(crate) neq: usize, // [i32] number of rows = number of columns = n_equation - pub(crate) pos: usize, // [i32] current index => nnz in the end - pub(crate) max: usize, // [i32] max allowed number of entries (may be > nnz) - pub(crate) indices_i: Vec, // [nnz] indices i - pub(crate) indices_j: Vec, // [nnz] indices j - pub(crate) values_aij: Vec, // [nnz] values aij -} - -impl SparseTriplet { - /// Creates a new SparseTriplet representing a sparse matrix - /// - /// ```text - /// trip := sparse(a) - /// (max) (nrow,ncol) - /// ``` - /// - /// # Input - /// - /// * `neq` -- The number of rows (= ncol) of the sparse matrix - /// * `max` -- The maximum number fo non-zero (nnz) values in the sparse matrix, - /// including entries with repeated indices - /// **Note:** This value must be greater than or equal to the actual nnz. - /// - /// # Example - /// - /// ``` - /// use russell_sparse::{SparseTriplet, StrError}; - /// - /// fn main() -> Result<(), StrError> { - /// let (neq, nnz) = (3, 4); - /// let trip = SparseTriplet::new(neq, nnz)?; - /// Ok(()) - /// } - /// ``` - pub fn new(neq: usize, max: usize) -> Result { - if neq == 0 || max == 0 { - return Err("neq and max must be greater than zero"); - } - Ok(SparseTriplet { - neq, - pos: 0, - max, - indices_i: vec![0; max], - indices_j: vec![0; max], - values_aij: vec![0.0; max], - }) - } - - /// Puts the next triple (i,j,aij) into the Triplet - /// - /// # Example - /// - /// ``` - /// use russell_sparse::{SparseTriplet, StrError}; - /// - /// fn main() -> Result<(), StrError> { - /// let (neq, nnz) = (3, 4); - /// let mut trip = SparseTriplet::new(neq, nnz)?; - /// trip.put(0, 0, 1.0)?; - /// trip.put(1, 1, 2.0)?; - /// trip.put(2, 2, 3.0)?; - /// trip.put(0, 1, 4.0)?; - /// Ok(()) - /// } - /// ``` - pub fn put(&mut self, i: usize, j: usize, aij: f64) -> Result<(), StrError> { - if i >= self.neq { - return Err("sparse matrix row index is out of bounds"); - } - if j >= self.neq { - return Err("sparse matrix column index is out of bounds"); - } - if self.pos >= self.max { - return Err("current nnz (number of non-zeros) reached maximum limit"); - } - let i_i32 = to_i32(i); - let j_i32 = to_i32(j); - self.indices_i[self.pos] = i_i32; - self.indices_j[self.pos] = j_i32; - self.values_aij[self.pos] = aij; - self.pos += 1; - Ok(()) - } - - /// Returns the (nrow = ncol) dimensions of the matrix represented by this Triplet - /// - /// # Example - /// - /// ``` - /// use russell_sparse::{SparseTriplet, StrError}; - /// - /// fn main() -> Result<(), StrError> { - /// let (neq, nnz) = (2, 1); - /// let trip = SparseTriplet::new(neq, nnz)?; - /// assert_eq!(trip.neq(), 2); - /// Ok(()) - /// } - /// ``` - pub fn neq(&self) -> usize { - self.neq - } - - /// Returns the (current) number of non-zero values (nnz) - /// - /// # Example - /// - /// ``` - /// use russell_sparse::{SparseTriplet, StrError}; - /// - /// fn main() -> Result<(), StrError> { - /// let (neq, nnz) = (2, 1); - /// let mut trip = SparseTriplet::new(neq, neq)?; - /// assert_eq!(trip.nnz_current(), 0); - /// trip.put(0, 0, 1.0); - /// assert_eq!(trip.nnz_current(), 1); - /// Ok(()) - /// } - /// ``` - pub fn nnz_current(&self) -> usize { - self.pos - } - - /// Returns the maximum allowed number of non-zero values (max) - /// - /// # Example - /// - /// ``` - /// use russell_sparse::{SparseTriplet, StrError}; - /// - /// fn main() -> Result<(), StrError> { - /// let (neq, nnz) = (2, 1); - /// let trip = SparseTriplet::new(neq, nnz)?; - /// assert_eq!(trip.nnz_maximum(), 1); - /// Ok(()) - /// } - /// ``` - pub fn nnz_maximum(&self) -> usize { - self.max - } - - /// Resets the position of the current non-zero value, allowing using "put" from scratch - /// - /// # Example - /// - /// ``` - /// use russell_sparse::{SparseTriplet, StrError}; - /// - /// fn main() -> Result<(), StrError> { - /// let (neq, nnz) = (3, 4); - /// let mut trip = SparseTriplet::new(neq, nnz)?; - /// trip.put(0, 0, 1.0)?; - /// trip.put(1, 1, 2.0)?; - /// trip.put(2, 2, 3.0)?; - /// trip.put(0, 1, 4.0)?; - /// assert_eq!(trip.nnz_current(), 4); - /// trip.reset(); - /// assert_eq!(trip.nnz_current(), 0); - /// Ok(()) - /// } - /// ``` - pub fn reset(&mut self) { - self.pos = 0; - } - - /// Returns the Matrix corresponding to this Triplet - /// - /// Note: this function calls [SparseTriplet::to_matrix]. - /// - /// ``` - /// use russell_sparse::{SparseTriplet, StrError}; - /// - /// fn main() -> Result<(), StrError> { - /// // define (4 x 4) sparse matrix with 6+1 non-zero values - /// // (with an extra ij-repeated entry) - /// let (neq, nnz) = (4, 7); - /// let mut trip = SparseTriplet::new(neq, nnz)?; - /// trip.put(0, 0, 0.5)?; // (0, 0, a00/2) - /// trip.put(0, 0, 0.5)?; // (0, 0, a00/2) - /// trip.put(0, 1, 2.0)?; - /// trip.put(1, 0, 3.0)?; - /// trip.put(1, 1, 4.0)?; - /// trip.put(2, 2, 5.0)?; - /// trip.put(3, 3, 6.0)?; - /// - /// // convert to matrix - /// let a = trip.as_matrix(); - /// let correct = "┌ ┐\n\ - /// │ 1 2 0 0 │\n\ - /// │ 3 4 0 0 │\n\ - /// │ 0 0 5 0 │\n\ - /// │ 0 0 0 6 │\n\ - /// └ ┘"; - /// assert_eq!(format!("{}", a), correct); - /// Ok(()) - /// } - /// ``` - pub fn as_matrix(&self) -> Matrix { - let mut a = Matrix::new(self.neq, self.neq); - self.to_matrix(&mut a).unwrap(); - a - } - - /// Converts the triplet data to a matrix, up to a limit - /// - /// Note: see the function [SparseTriplet::as_matrix] that returns the Matrix already. - /// - /// # Input - /// - /// `a` -- (nrow_max, ncol_max) matrix to hold the triplet data. - /// The output matrix may have fewer rows or fewer columns than the triplet data. - /// - /// # Example - /// - /// ``` - /// use russell_lab::{Matrix}; - /// use russell_sparse::{SparseTriplet, StrError}; - /// - /// fn main() -> Result<(), StrError> { - /// // define (4 x 4) sparse matrix with 6+1 non-zero values - /// // (with an extra ij-repeated entry) - /// let (neq, nnz) = (4, 7); - /// let mut trip = SparseTriplet::new(neq, nnz)?; - /// trip.put(0, 0, 0.5)?; // (0, 0, a00/2) - /// trip.put(0, 0, 0.5)?; // (0, 0, a00/2) - /// trip.put(0, 1, 2.0)?; - /// trip.put(1, 0, 3.0)?; - /// trip.put(1, 1, 4.0)?; - /// trip.put(2, 2, 5.0)?; - /// trip.put(3, 3, 6.0)?; - /// - /// // convert the first (3 x 3) values - /// let mut a = Matrix::new(3, 3); - /// trip.to_matrix(&mut a)?; - /// let correct = "┌ ┐\n\ - /// │ 1 2 0 │\n\ - /// │ 3 4 0 │\n\ - /// │ 0 0 5 │\n\ - /// └ ┘"; - /// assert_eq!(format!("{}", a), correct); - /// - /// // convert the first (4 x 4) values - /// let mut b = Matrix::new(4, 4); - /// trip.to_matrix(&mut b)?; - /// let correct = "┌ ┐\n\ - /// │ 1 2 0 0 │\n\ - /// │ 3 4 0 0 │\n\ - /// │ 0 0 5 0 │\n\ - /// │ 0 0 0 6 │\n\ - /// └ ┘"; - /// assert_eq!(format!("{}", b), correct); - /// Ok(()) - /// } - /// ``` - pub fn to_matrix(&self, a: &mut Matrix) -> Result<(), StrError> { - let (m, n) = a.dims(); - if m > self.neq || n > self.neq { - return Err("wrong matrix dimensions"); - } - let m_i32 = to_i32(m); - let n_i32 = to_i32(n); - a.fill(0.0); - for p in 0..self.pos { - if self.indices_i[p] < m_i32 && self.indices_j[p] < n_i32 { - let (i, j) = (self.indices_i[p] as usize, self.indices_j[p] as usize); - a.add(i, j, self.values_aij[p]); - } - } - Ok(()) - } - - /// Performs the matrix-vector multiplication - /// - /// ```text - /// v := a ⋅ u - /// (m) (m,n) (n) - /// ``` - /// - /// # Input - /// - /// * `triangular` -- must be set to true if the triplet stores - /// the components of the matrix in triangular format, i.e., - /// only the upper/lower diagonal values and the diagonal are stored - /// - /// # Note - /// - /// This method is not highly efficient but should useful in verifications. - /// - /// # Example - /// - /// ``` - /// use russell_lab::{Matrix, Vector}; - /// use russell_sparse::{SparseTriplet, StrError}; - /// - /// fn main() -> Result<(), StrError> { - /// // set sparse matrix (3 x 3) with 6 non-zeros - /// let (neq, nnz) = (3, 6); - /// let mut trip = SparseTriplet::new(neq, nnz)?; - /// trip.put(0, 0, 1.0)?; - /// trip.put(1, 0, 2.0)?; - /// trip.put(1, 1, 3.0)?; - /// trip.put(2, 0, 4.0)?; - /// trip.put(2, 1, 5.0)?; - /// trip.put(2, 2, 6.0)?; - /// - /// // check matrix - /// let mut a = Matrix::new(neq, neq); - /// trip.to_matrix(&mut a)?; - /// let correct_a = "┌ ┐\n\ - /// │ 1 0 0 │\n\ - /// │ 2 3 0 │\n\ - /// │ 4 5 6 │\n\ - /// └ ┘"; - /// assert_eq!(format!("{}", a), correct_a); - /// - /// // perform mat-vec-mul - /// let u = Vector::from(&[1.0, 1.0, 1.0]); - /// let v = trip.mat_vec_mul(&u, false)?; - /// - /// // check vector - /// let correct_v = "┌ ┐\n\ - /// │ 1 │\n\ - /// │ 5 │\n\ - /// │ 15 │\n\ - /// └ ┘"; - /// assert_eq!(format!("{}", v), correct_v); - /// Ok(()) - /// } - /// ``` - pub fn mat_vec_mul(&self, u: &Vector, triangular: bool) -> Result { - if u.dim() != self.neq { - return Err("u.ndim must equal neq"); - } - let mut v = Vector::new(self.neq); - for p in 0..self.pos { - let i = self.indices_i[p] as usize; - let j = self.indices_j[p] as usize; - let aij = self.values_aij[p]; - v[i] += aij * u[j]; - if triangular && i != j { - v[j] += aij * u[i]; - } - } - Ok(v) - } -} - -impl fmt::Display for SparseTriplet { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!( - f, - "\x20\x20\x20\x20\"neq\": {},\n\ - \x20\x20\x20\x20\"nnz_current\": {},\n\ - \x20\x20\x20\x20\"nnz_maximum\": {},\n", - self.neq, self.pos, self.max, - ) - .unwrap(); - Ok(()) - } -} - -//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -#[cfg(test)] -mod tests { - use super::SparseTriplet; - use russell_chk::vec_approx_eq; - use russell_lab::{Matrix, Vector}; - - #[test] - fn new_fails_on_wrong_input() { - assert_eq!( - SparseTriplet::new(0, 3).err(), - Some("neq and max must be greater than zero") - ); - assert_eq!( - SparseTriplet::new(3, 0).err(), - Some("neq and max must be greater than zero") - ); - } - - #[test] - fn new_works() { - let trip = SparseTriplet::new(3, 5).unwrap(); - assert_eq!(trip.neq, 3); - assert_eq!(trip.pos, 0); - assert_eq!(trip.max, 5); - } - - #[test] - fn put_fails_on_wrong_values() { - let mut trip = SparseTriplet::new(1, 1).unwrap(); - assert_eq!( - trip.put(1, 0, 0.0).err(), - Some("sparse matrix row index is out of bounds") - ); - assert_eq!( - trip.put(0, 1, 0.0).err(), - Some("sparse matrix column index is out of bounds") - ); - assert_eq!(trip.put(0, 0, 0.0).err(), None); // << will tak all spots - assert_eq!( - trip.put(0, 0, 0.0).err(), - Some("current nnz (number of non-zeros) reached maximum limit") - ); - } - - #[test] - fn put_works() { - let mut trip = SparseTriplet::new(3, 5).unwrap(); - trip.put(0, 0, 1.0).unwrap(); - assert_eq!(trip.pos, 1); - trip.put(0, 1, 2.0).unwrap(); - assert_eq!(trip.pos, 2); - trip.put(1, 0, 3.0).unwrap(); - assert_eq!(trip.pos, 3); - trip.put(1, 1, 4.0).unwrap(); - assert_eq!(trip.pos, 4); - trip.put(2, 2, 5.0).unwrap(); - assert_eq!(trip.pos, 5); - } - - #[test] - fn getters_and_reset_work() { - let mut trip = SparseTriplet::new(2, 4).unwrap(); - assert_eq!(trip.nnz_current(), 0); - trip.put(0, 0, 1.0).unwrap(); - trip.put(0, 1, 4.0).unwrap(); - trip.put(1, 0, 2.0).unwrap(); - trip.put(1, 1, 3.0).unwrap(); - assert_eq!(trip.neq(), 2); - assert_eq!(trip.nnz_current(), 4); - assert_eq!(trip.nnz_maximum(), 4); - trip.reset(); - assert_eq!(trip.nnz_current(), 0); - } - - #[test] - fn to_matrix_fails_on_wrong_dims() { - let trip = SparseTriplet::new(1, 1).unwrap(); - let mut a_2x1 = Matrix::new(2, 1); - let mut a_1x2 = Matrix::new(1, 2); - assert_eq!(trip.to_matrix(&mut a_2x1), Err("wrong matrix dimensions")); - assert_eq!(trip.to_matrix(&mut a_1x2), Err("wrong matrix dimensions")); - } - - #[test] - fn to_matrix_works() { - let mut trip = SparseTriplet::new(3, 5).unwrap(); - trip.put(0, 0, 1.0).unwrap(); - trip.put(0, 1, 2.0).unwrap(); - trip.put(1, 0, 3.0).unwrap(); - trip.put(1, 1, 4.0).unwrap(); - trip.put(2, 2, 5.0).unwrap(); - let mut a = Matrix::new(3, 3); - trip.to_matrix(&mut a).unwrap(); - assert_eq!(a.get(0, 0), 1.0); - assert_eq!(a.get(0, 1), 2.0); - assert_eq!(a.get(1, 0), 3.0); - assert_eq!(a.get(1, 1), 4.0); - assert_eq!(a.get(2, 2), 5.0); - let mut b = Matrix::new(2, 1); - trip.to_matrix(&mut b).unwrap(); - assert_eq!(b.get(0, 0), 1.0); - assert_eq!(b.get(1, 0), 3.0); - // using as_matrix - let bb = trip.as_matrix(); - assert_eq!(bb.get(0, 0), 1.0); - assert_eq!(bb.get(1, 0), 3.0); - } - - #[test] - fn to_matrix_with_duplicates_works() { - // allocate a square matrix - let (neq, nnz) = (5, 13); - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 1.0).unwrap(); // << (0, 0, a00/2) - trip.put(0, 0, 1.0).unwrap(); // << (0, 0, a00/2) - trip.put(1, 0, 3.0).unwrap(); - trip.put(0, 1, 3.0).unwrap(); - trip.put(2, 1, -1.0).unwrap(); - trip.put(4, 1, 4.0).unwrap(); - trip.put(1, 2, 4.0).unwrap(); - trip.put(2, 2, -3.0).unwrap(); - trip.put(3, 2, 1.0).unwrap(); - trip.put(4, 2, 2.0).unwrap(); - trip.put(2, 3, 2.0).unwrap(); - trip.put(1, 4, 6.0).unwrap(); - trip.put(4, 4, 1.0).unwrap(); - - // print matrix - let mut a = Matrix::new(neq, neq); - trip.to_matrix(&mut a).unwrap(); - let correct = "┌ ┐\n\ - │ 2 3 0 0 0 │\n\ - │ 3 0 4 0 6 │\n\ - │ 0 -1 -3 2 0 │\n\ - │ 0 0 1 0 0 │\n\ - │ 0 4 2 0 1 │\n\ - └ ┘"; - assert_eq!(format!("{}", a), correct); - } - - #[test] - fn mat_vec_mul_fails_on_wrong_input() { - let trip = SparseTriplet::new(2, 1).unwrap(); - let u = Vector::new(3); - assert_eq!(trip.mat_vec_mul(&u, false).err(), Some("u.ndim must equal neq")); - } - - #[test] - fn mat_vec_mul_works() { - // 1.0 2.0 3.0 - // 0.1 0.2 0.3 - // 10.0 20.0 30.0 - let mut trip = SparseTriplet::new(3, 9).unwrap(); - trip.put(0, 0, 1.0).unwrap(); - trip.put(0, 1, 2.0).unwrap(); - trip.put(0, 2, 3.0).unwrap(); - trip.put(1, 0, 0.1).unwrap(); - trip.put(1, 1, 0.2).unwrap(); - trip.put(1, 2, 0.3).unwrap(); - trip.put(2, 0, 10.0).unwrap(); - trip.put(2, 1, 20.0).unwrap(); - trip.put(2, 2, 30.0).unwrap(); - let u = Vector::from(&[0.1, 0.2, 0.3]); - let correct_v = &[1.4, 0.14, 14.0]; - let v = trip.mat_vec_mul(&u, false).unwrap(); - vec_approx_eq(v.as_data(), correct_v, 1e-15); - } - - #[test] - fn mat_vec_mul_sym_part_works() { - // 2 - // 1 2 sym - // 1 2 9 - // 3 1 1 7 - // 2 1 5 1 8 - let (neq, nnz) = (5, 15); - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 2.0).unwrap(); - trip.put(1, 1, 2.0).unwrap(); - trip.put(2, 2, 9.0).unwrap(); - trip.put(3, 3, 7.0).unwrap(); - trip.put(4, 4, 8.0).unwrap(); - - trip.put(1, 0, 1.0).unwrap(); - - trip.put(2, 0, 1.0).unwrap(); - trip.put(2, 1, 2.0).unwrap(); - - trip.put(3, 0, 3.0).unwrap(); - trip.put(3, 1, 1.0).unwrap(); - trip.put(3, 2, 1.0).unwrap(); - - trip.put(4, 0, 2.0).unwrap(); - trip.put(4, 1, 1.0).unwrap(); - trip.put(4, 2, 5.0).unwrap(); - trip.put(4, 3, 1.0).unwrap(); - let u = Vector::from(&[-629.0 / 98.0, 237.0 / 49.0, -53.0 / 49.0, 62.0 / 49.0, 23.0 / 14.0]); - let correct_v = &[-2.0, 4.0, 3.0, -5.0, 1.0]; - let v = trip.mat_vec_mul(&u, true).unwrap(); - vec_approx_eq(v.as_data(), correct_v, 1e-14); - } - - #[test] - fn mat_vec_mul_sym_full_works() { - // 2 1 1 3 2 - // 1 2 2 1 1 - // 1 2 9 1 5 - // 3 1 1 7 1 - // 2 1 5 1 8 - let (neq, nnz) = (5, 25); - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 2.0).unwrap(); - trip.put(1, 1, 2.0).unwrap(); - trip.put(2, 2, 9.0).unwrap(); - trip.put(3, 3, 7.0).unwrap(); - trip.put(4, 4, 8.0).unwrap(); - - trip.put(1, 0, 1.0).unwrap(); - trip.put(0, 1, 1.0).unwrap(); - - trip.put(2, 0, 1.0).unwrap(); - trip.put(0, 2, 1.0).unwrap(); - trip.put(2, 1, 2.0).unwrap(); - trip.put(1, 2, 2.0).unwrap(); - - trip.put(3, 0, 3.0).unwrap(); - trip.put(0, 3, 3.0).unwrap(); - trip.put(3, 1, 1.0).unwrap(); - trip.put(1, 3, 1.0).unwrap(); - trip.put(3, 2, 1.0).unwrap(); - trip.put(2, 3, 1.0).unwrap(); - - trip.put(4, 0, 2.0).unwrap(); - trip.put(0, 4, 2.0).unwrap(); - trip.put(4, 1, 1.0).unwrap(); - trip.put(1, 4, 1.0).unwrap(); - trip.put(4, 2, 5.0).unwrap(); - trip.put(2, 4, 5.0).unwrap(); - trip.put(4, 3, 1.0).unwrap(); - trip.put(3, 4, 1.0).unwrap(); - let u = Vector::from(&[-629.0 / 98.0, 237.0 / 49.0, -53.0 / 49.0, 62.0 / 49.0, 23.0 / 14.0]); - let correct_v = &[-2.0, 4.0, 3.0, -5.0, 1.0]; - let v = trip.mat_vec_mul(&u, false).unwrap(); - vec_approx_eq(v.as_data(), correct_v, 1e-14); - } - - #[test] - fn mat_vec_mul_pos_def_works() { - // 2 -1 2 ... - // -1 2 -1 => -1 2 - // -1 2 -1 2 - let (neq, nnz) = (3, 5); - let mut trip = SparseTriplet::new(neq, nnz).unwrap(); - trip.put(0, 0, 2.0).unwrap(); - trip.put(1, 1, 2.0).unwrap(); - trip.put(2, 2, 2.0).unwrap(); - trip.put(1, 0, -1.0).unwrap(); - trip.put(2, 1, -1.0).unwrap(); - let u = Vector::from(&[5.0, 8.0, 7.0]); - let correct_v = &[2.0, 4.0, 6.0]; - let v = trip.mat_vec_mul(&u, true).unwrap(); - vec_approx_eq(v.as_data(), correct_v, 1e-15); - } - - #[test] - fn display_trait_works() { - let trip = SparseTriplet::new(3, 1).unwrap(); - let correct: &str = "\x20\x20\x20\x20\"neq\": 3,\n\ - \x20\x20\x20\x20\"nnz_current\": 0,\n\ - \x20\x20\x20\x20\"nnz_maximum\": 1,\n"; - assert_eq!(format!("{}", trip), correct); - } -} diff --git a/russell_sparse/src/verify_lin_sys.rs b/russell_sparse/src/verify_lin_sys.rs index 1f892d7a..b52d5f10 100644 --- a/russell_sparse/src/verify_lin_sys.rs +++ b/russell_sparse/src/verify_lin_sys.rs @@ -1,8 +1,7 @@ -use super::SparseTriplet; +use super::CooMatrix; use crate::StrError; -use russell_lab::{format_nanoseconds, vec_norm, vec_update, Norm, Stopwatch, Vector}; +use russell_lab::{vec_norm, vec_update, Norm, Stopwatch, Vector}; use russell_openblas::{idamax, to_i32}; -use std::fmt; /// Verifies the linear system a ⋅ x = rhs pub struct VerifyLinSys { @@ -24,20 +23,20 @@ impl VerifyLinSys { /// /// ``` /// use russell_lab::{Matrix, Vector}; - /// use russell_sparse::{SparseTriplet, VerifyLinSys, StrError}; + /// use russell_sparse::{CooMatrix, Layout, StrError, VerifyLinSys}; /// /// fn main() -> Result<(), StrError> { /// // set sparse matrix (3 x 3) with 4 non-zeros - /// let (neq, nnz) = (3, 4); - /// let mut trip = SparseTriplet::new(neq, nnz)?; - /// trip.put(0, 0, 1.0)?; - /// trip.put(0, 2, 4.0)?; - /// trip.put(1, 1, 2.0)?; - /// trip.put(2, 2, 3.0)?; + /// let (nrow, ncol, nnz) = (3, 3, 4); + /// let mut coo = CooMatrix::new(Layout::Full, nrow, ncol, nnz)?; + /// coo.put(0, 0, 1.0)?; + /// coo.put(0, 2, 4.0)?; + /// coo.put(1, 1, 2.0)?; + /// coo.put(2, 2, 3.0)?; /// /// // check matrix - /// let mut a = Matrix::new(neq, neq); - /// trip.to_matrix(&mut a)?; + /// let mut a = Matrix::new(nrow, nrow); + /// coo.to_matrix(&mut a)?; /// let correct_a = "┌ ┐\n\ /// │ 1 0 4 │\n\ /// │ 0 2 0 │\n\ @@ -48,7 +47,7 @@ impl VerifyLinSys { /// // verify lin-sys /// let x = Vector::from(&[1.0, 1.0, 1.0]); /// let rhs = Vector::from(&[5.0, 2.0, 3.0]); - /// let verify = VerifyLinSys::new(&trip, &x, &rhs, false)?; + /// let verify = VerifyLinSys::new(&coo, &x, &rhs)?; /// assert_eq!(verify.max_abs_a, 4.0); /// assert_eq!(verify.max_abs_ax, 5.0); /// assert_eq!(verify.max_abs_diff, 0.0); @@ -57,20 +56,23 @@ impl VerifyLinSys { /// Ok(()) /// } /// ``` - pub fn new(trip: &SparseTriplet, x: &Vector, rhs: &Vector, triangular: bool) -> Result { - if x.dim() != trip.neq || rhs.dim() != trip.neq { + pub fn new(coo: &CooMatrix, x: &Vector, rhs: &Vector) -> Result { + if coo.nrow != coo.ncol { + return Err("CooMatrix mut be square"); + } + if x.dim() != coo.nrow || rhs.dim() != coo.nrow { return Err("vector dimensions are incompatible"); } // start stopwatch let mut sw = Stopwatch::new(""); // compute max_abs_a - let nnz = to_i32(trip.pos); - let idx = idamax(nnz, &trip.values_aij, 1); - let max_abs_a = f64::abs(trip.values_aij[idx as usize]); + let nnz = to_i32(coo.pos); + let idx = idamax(nnz, &coo.values_aij, 1); + let max_abs_a = f64::abs(coo.values_aij[idx as usize]); // compute max_abs_ax - let mut ax = trip.mat_vec_mul(&x, triangular).unwrap(); // already checked + let mut ax = coo.mat_vec_mul(&x).unwrap(); // already checked let max_abs_ax = vec_norm(&ax, Norm::Max); // compute max_abs_diff @@ -94,48 +96,27 @@ impl VerifyLinSys { } } -impl fmt::Display for VerifyLinSys { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!( - f, - "\x20\x20\x20\x20\"maxAbsA\": {},\n\ - \x20\x20\x20\x20\"maxAbsAx\": {},\n\ - \x20\x20\x20\x20\"maxAbsDiff\": {:e},\n\ - \x20\x20\x20\x20\"relativeError\": {:e},\n\ - \x20\x20\x20\x20\"timeCheckNs\": {},\n\ - \x20\x20\x20\x20\"timeCheckStr\": \"{}\"", - self.max_abs_a, - self.max_abs_ax, - self.max_abs_diff, - self.relative_error, - self.time_check, - format_nanoseconds(self.time_check), - ) - .unwrap(); - Ok(()) - } -} - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #[cfg(test)] mod tests { - use super::{SparseTriplet, VerifyLinSys}; + use super::{CooMatrix, VerifyLinSys}; + use crate::Layout; use russell_lab::Vector; #[test] fn new_fails_on_wrong_vectors() { - let trip = SparseTriplet::new(1, 1).unwrap(); + let coo = CooMatrix::new(Layout::Full, 1, 1, 1).unwrap(); let x = Vector::new(2); let rhs = Vector::new(3); let x_wrong = Vector::new(3); let rhs_wrong = Vector::new(2); assert_eq!( - VerifyLinSys::new(&trip, &x_wrong, &rhs, false).err(), + VerifyLinSys::new(&coo, &x_wrong, &rhs).err(), Some("vector dimensions are incompatible") ); assert_eq!( - VerifyLinSys::new(&trip, &x, &rhs_wrong, false).err(), + VerifyLinSys::new(&coo, &x, &rhs_wrong).err(), Some("vector dimensions are incompatible") ); } @@ -145,41 +126,23 @@ mod tests { // | 1 3 -2 | // | 3 5 6 | // | 2 4 3 | - let mut trip = SparseTriplet::new(3, 9).unwrap(); - trip.put(0, 0, 1.0).unwrap(); - trip.put(0, 1, 3.0).unwrap(); - trip.put(0, 2, -2.0).unwrap(); - trip.put(1, 0, 3.0).unwrap(); - trip.put(1, 1, 5.0).unwrap(); - trip.put(1, 2, 6.0).unwrap(); - trip.put(2, 0, 2.0).unwrap(); - trip.put(2, 1, 4.0).unwrap(); - trip.put(2, 2, 3.0).unwrap(); + let mut coo = CooMatrix::new(Layout::Full, 3, 3, 9).unwrap(); + coo.put(0, 0, 1.0).unwrap(); + coo.put(0, 1, 3.0).unwrap(); + coo.put(0, 2, -2.0).unwrap(); + coo.put(1, 0, 3.0).unwrap(); + coo.put(1, 1, 5.0).unwrap(); + coo.put(1, 2, 6.0).unwrap(); + coo.put(2, 0, 2.0).unwrap(); + coo.put(2, 1, 4.0).unwrap(); + coo.put(2, 2, 3.0).unwrap(); let x = Vector::from(&[-15.0, 8.0, 2.0]); let rhs = Vector::from(&[5.0, 7.0, 8.0]); - let verify = VerifyLinSys::new(&trip, &x, &rhs, false).unwrap(); + let verify = VerifyLinSys::new(&coo, &x, &rhs).unwrap(); assert_eq!(verify.max_abs_a, 6.0); assert_eq!(verify.max_abs_ax, 8.0); assert_eq!(verify.max_abs_diff, 0.0); assert_eq!(verify.relative_error, 0.0); assert!(verify.time_check > 0); } - - #[test] - fn display_trait_works() { - let mut trip = SparseTriplet::new(2, 2).unwrap(); - trip.put(0, 0, 1.0).unwrap(); - trip.put(1, 1, 1.0).unwrap(); - let x = Vector::from(&[1.0, 1.0]); - let rhs = Vector::from(&[1.0, 1.0]); - let mut verify = VerifyLinSys::new(&trip, &x, &rhs, false).unwrap(); - verify.time_check = 0; - let correct: &str = "\x20\x20\x20\x20\"maxAbsA\": 1,\n\ - \x20\x20\x20\x20\"maxAbsAx\": 1,\n\ - \x20\x20\x20\x20\"maxAbsDiff\": 0e0,\n\ - \x20\x20\x20\x20\"relativeError\": 0e0,\n\ - \x20\x20\x20\x20\"timeCheckNs\": 0,\n\ - \x20\x20\x20\x20\"timeCheckStr\": \"0ns\""; - assert_eq!(format!("{}", verify), correct); - } } diff --git a/russell_sparse/src/write_matrix_market.rs b/russell_sparse/src/write_matrix_market.rs new file mode 100644 index 00000000..333b7d92 --- /dev/null +++ b/russell_sparse/src/write_matrix_market.rs @@ -0,0 +1,238 @@ +use super::{CsrMatrix, Layout}; +use crate::StrError; +use std::ffi::OsStr; +use std::fmt::Write; +use std::fs::{self, File}; +use std::io::Write as IoWrite; +use std::path::Path; + +impl CsrMatrix { + /// Writes a MatrixMarket file from a CooMatrix + /// + /// # Input + /// + /// * `full_path` -- may be a String, &str, or Path + /// * `vismatrix` -- generate a SMAT file for Vismatrix instead of a MatrixMarket + /// + /// **Note:** The vismatrix format is is like the MatrixMarket format + /// without the header, and the indices start at zero. + /// + /// # References + /// + /// * MatrixMarket: + /// * Vismatrix: + /// + /// # Examples + /// + /// ``` + /// use russell_sparse::prelude::*; + /// use russell_sparse::StrError; + /// + /// const SAVE_FILE: bool = false; + /// + /// fn main() -> Result<(), StrError> { + /// // allocate a square matrix and store as CSR matrix + /// // ┌ ┐ + /// // │ 2 3 0 0 0 │ + /// // │ 3 0 4 0 6 │ + /// // │ 0 -1 -3 2 0 │ + /// // │ 0 0 1 0 0 │ + /// // │ 0 4 2 0 1 │ + /// // └ ┘ + /// let csr = CsrMatrix { + /// layout: Layout::Full, + /// nrow: 5, + /// ncol: 5, + /// row_pointers: vec![0, 2, 5, 8, 9, 12], + /// col_indices: vec![ + /// // p + /// 0, 1, // i = 0, count = 0, 1 + /// 0, 2, 4, // i = 1, count = 2, 3, 4 + /// 1, 2, 3, // i = 2, count = 5, 6, 7 + /// 2, // i = 3, count = 8 + /// 1, 2, 4, // i = 4, count = 9, 10, 11 + /// // count = 12 + /// ], + /// values: vec![ + /// // p + /// 2.0, 3.0, // i = 0, count = 0, 1 + /// 3.0, 4.0, 6.0, // i = 1, count = 2, 3, 4 + /// -1.0, -3.0, 2.0, // i = 2, count = 5, 6, 7 + /// 1.0, // i = 3, count = 8 + /// 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 + /// // count = 12 + /// ], + /// }; + /// if SAVE_FILE { + /// let full_path = "/tmp/russell_sparse/doc-example-vismatrix.smat"; + /// csr.write_matrix_market(full_path, true)?; + /// } + /// Ok(()) + /// } + /// ``` + /// + /// By running `vismatrix doc-example-vismatrix.smat` you get the following screen: + /// + /// ![doc-example-vismatrix](https://raw.githubusercontent.com/cpmech/russell/main/russell_sparse/data/figures/doc-example-vismatrix.svg) + pub fn write_matrix_market

(&self, full_path: &P, vismatrix: bool) -> Result<(), StrError> + where + P: AsRef + ?Sized, + { + // output buffer + let mut buffer = String::new(); + + // write to buffer + let d = if vismatrix { 0 } else { 1 }; + let nnz = self.values.len(); + if !vismatrix { + if self.layout == Layout::Full { + write!(&mut buffer, "%%MatrixMarket matrix coordinate real general\n").unwrap() + } else { + write!(&mut buffer, "%%MatrixMarket matrix coordinate real symmetric\n").unwrap() + } + } + write!(&mut buffer, "{} {} {}\n", self.nrow, self.ncol, nnz).unwrap(); + for i in 0..self.nrow { + for p in self.row_pointers[i]..self.row_pointers[i + 1] { + let j = self.col_indices[p as usize] as usize; + let aij = self.values[p as usize]; + write!(&mut buffer, "{} {} {:?}\n", i + d, j + d, aij).unwrap() + } + } + + // create directory + let path = Path::new(full_path); + if let Some(p) = path.parent() { + fs::create_dir_all(p).map_err(|_| "cannot create directory")?; + } + + // write file + let mut file = File::create(path).map_err(|_| "cannot create file")?; + file.write_all(buffer.as_bytes()).map_err(|_| "cannot write file")?; + + // force sync + file.sync_all().map_err(|_| "cannot sync file")?; + Ok(()) + } +} + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +#[cfg(test)] +mod tests { + use crate::{CsrMatrix, Layout}; + use std::fs; + + #[test] + fn write_matrix_market_works() { + // allocate a square matrix and store as CSR matrix + // ┌ ┐ + // │ 2 3 0 0 0 │ + // │ 3 0 4 0 6 │ + // │ 0 -1 -3 2 0 │ + // │ 0 0 1 0 0 │ + // │ 0 4 2 0 1 │ + // └ ┘ + let csr = CsrMatrix { + layout: Layout::Full, + nrow: 5, + ncol: 5, + row_pointers: vec![0, 2, 5, 8, 9, 12], + col_indices: vec![ + // p + 0, 1, // i = 0, count = 0, 1 + 0, 2, 4, // i = 1, count = 2, 3, 4 + 1, 2, 3, // i = 2, count = 5, 6, 7 + 2, // i = 3, count = 8 + 1, 2, 4, // i = 4, count = 9, 10, 11 + // count = 12 + ], + values: vec![ + // p + 2.0, 3.0, // i = 0, count = 0, 1 + 3.0, 4.0, 6.0, // i = 1, count = 2, 3, 4 + -1.0, -3.0, 2.0, // i = 2, count = 5, 6, 7 + 1.0, // i = 3, count = 8 + 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 + // count = 12 + ], + }; + + let full_path = "/tmp/russell_sparse/test_write_matrix_market.mtx"; + csr.write_matrix_market(full_path, false).unwrap(); + let contents = fs::read_to_string(full_path).map_err(|_| "cannot open file").unwrap(); + assert_eq!( + contents, + "%%MatrixMarket matrix coordinate real general\n\ + 5 5 12\n\ + 1 1 2.0\n\ + 1 2 3.0\n\ + 2 1 3.0\n\ + 2 3 4.0\n\ + 2 5 6.0\n\ + 3 2 -1.0\n\ + 3 3 -3.0\n\ + 3 4 2.0\n\ + 4 3 1.0\n\ + 5 2 4.0\n\ + 5 3 2.0\n\ + 5 5 1.0\n" + ); + } + + #[test] + fn write_matrix_market_vismatrix_works() { + // allocate a square matrix and store as CSR matrix + // ┌ ┐ + // │ 2 3 0 0 0 │ + // │ 3 0 4 0 6 │ + // │ 0 -1 -3 2 0 │ + // │ 0 0 1 0 0 │ + // │ 0 4 2 0 1 │ + // └ ┘ + let csr = CsrMatrix { + layout: Layout::Full, + nrow: 5, + ncol: 5, + row_pointers: vec![0, 2, 5, 8, 9, 12], + col_indices: vec![ + // p + 0, 1, // i = 0, count = 0, 1 + 0, 2, 4, // i = 1, count = 2, 3, 4 + 1, 2, 3, // i = 2, count = 5, 6, 7 + 2, // i = 3, count = 8 + 1, 2, 4, // i = 4, count = 9, 10, 11 + // count = 12 + ], + values: vec![ + // p + 2.0, 3.0, // i = 0, count = 0, 1 + 3.0, 4.0, 6.0, // i = 1, count = 2, 3, 4 + -1.0, -3.0, 2.0, // i = 2, count = 5, 6, 7 + 1.0, // i = 3, count = 8 + 4.0, 2.0, 1.0, // i = 4, count = 9, 10, 11 + // count = 12 + ], + }; + + let full_path = "/tmp/russell_sparse/test_write_matrix_market_vismatrix.smat"; + csr.write_matrix_market(full_path, true).unwrap(); + let contents = fs::read_to_string(full_path).map_err(|_| "cannot open file").unwrap(); + assert_eq!( + contents, + "5 5 12\n\ + 0 0 2.0\n\ + 0 1 3.0\n\ + 1 0 3.0\n\ + 1 2 4.0\n\ + 1 4 6.0\n\ + 2 1 -1.0\n\ + 2 2 -3.0\n\ + 2 3 2.0\n\ + 3 2 1.0\n\ + 4 1 4.0\n\ + 4 2 2.0\n\ + 4 4 1.0\n" + ); + } +} diff --git a/russell_sparse/tests/test_nonlinear_system.rs b/russell_sparse/tests/test_nonlinear_system.rs index 466c5f3a..a2beb8f9 100644 --- a/russell_sparse/tests/test_nonlinear_system.rs +++ b/russell_sparse/tests/test_nonlinear_system.rs @@ -1,6 +1,6 @@ use russell_chk::{deriv_central5, vec_approx_eq}; use russell_lab::{mat_approx_eq, vec_norm, vec_update, Matrix, Norm, Vector}; -use russell_sparse::{ConfigSolver, LinSolKind, Solver, SparseTriplet, StrError}; +use russell_sparse::{ConfigSolver, CooMatrix, Layout, LinSolKind, Solver, StrError}; fn calc_residual(rr: &mut Vector, uu: &Vector) { let (d1, d2, d3, d4) = (uu[0], uu[1], uu[2], uu[3]); @@ -10,7 +10,7 @@ fn calc_residual(rr: &mut Vector, uu: &Vector) { rr[3] = -9.0 * d1 + 4.0 * d1 * d4 * d4 * d4 + 7.0 * d2 + 2.0 * d3 + 5.0 * d4 - 0.5; } -fn calc_jacobian(jj: &mut SparseTriplet, uu: &Vector) -> Result<(), StrError> { +fn calc_jacobian(jj: &mut CooMatrix, uu: &Vector) -> Result<(), StrError> { let (d1, d2, d3, d4) = (uu[0], uu[1], uu[2], uu[3]); jj.reset(); @@ -67,7 +67,7 @@ fn check_jacobian() { } } let nnz = neq * neq; - let mut jj_tri = SparseTriplet::new(neq, nnz).unwrap(); + let mut jj_tri = CooMatrix::new(Layout::Full, neq, neq, nnz).unwrap(); calc_jacobian(&mut jj_tri, &uu).unwrap(); let mut jj_ana = Matrix::new(neq, neq); jj_tri.to_matrix(&mut jj_ana).unwrap(); @@ -80,7 +80,7 @@ fn solve_nonlinear_system(kind: LinSolKind) -> Result<(), StrError> { // config.verbose(); let (neq, nnz) = (4, 16); let mut solver = Solver::new(config, neq, nnz, None)?; - let mut jj = SparseTriplet::new(neq, nnz).unwrap(); + let mut jj = CooMatrix::new(Layout::Full, neq, neq, nnz).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/zscripts/memcheck.bash b/russell_sparse/zscripts/memcheck.bash new file mode 100755 index 00000000..dd58b04e --- /dev/null +++ b/russell_sparse/zscripts/memcheck.bash @@ -0,0 +1,5 @@ +#!/bin/bash + +cargo valgrind run --bin mem_check_build +cargo valgrind run --bin solve_matrix_market_build -- data/matrix_market/bfwb62.mtx +cargo valgrind run --bin solve_matrix_market_build -- data/matrix_market/bfwb62.mtx --mmp diff --git a/russell_sparse/zscripts/run-examples.bash b/russell_sparse/zscripts/run-examples.bash new file mode 100755 index 00000000..2bba7c66 --- /dev/null +++ b/russell_sparse/zscripts/run-examples.bash @@ -0,0 +1,7 @@ +#!/bin/bash + +for example in examples/*.rs; do + filename="$(basename "$example")" + filekey="${filename%%.*}" + cargo run --example $filekey +done diff --git a/russell_sparse/zscripts/run-solve-matrix-markert.bash b/russell_sparse/zscripts/run-solve-matrix-markert.bash new file mode 100755 index 00000000..f3e75e83 --- /dev/null +++ b/russell_sparse/zscripts/run-solve-matrix-markert.bash @@ -0,0 +1,4 @@ +#!/bin/bash + +cargo run --release --bin solve_matrix_market_build -- data/matrix_market/bfwb62.mtx +cargo run --release --bin solve_matrix_market_build -- data/matrix_market/bfwb62.mtx --mmp diff --git a/zscripts/cargo_workspaces.bash b/zscripts/cargo-workspaces.bash similarity index 100% rename from zscripts/cargo_workspaces.bash rename to zscripts/cargo-workspaces.bash