Impl compute determinant

.vscode/settings.json

@@ -15,6 +15,7 @@
     "dscal",
     "dsyev",
     "dsyrk",
+    "dtype",
     "Flannery",
     "IIIₛ",
     "ᵢⱼₖₗ",
@@ -25,12 +26,14 @@
     "ldvl",
     "ldvr",
     "ldvt",
+    "linalg",
     "lredrhs",
     "lrhs",
     "lsol",
     "nelt",
     "odyad",
     "PERMUT",
+    "RINFOG",
     "Schur",
     "Teukolsky",
     "tgamma",

russell_sparse/c_code/solver_mmp.h

@@ -9,9 +9,10 @@
 #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 RINFOG(i) rinfog[(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) {
@@ -82,7 +83,8 @@ int32_t solver_mmp_initialize(struct SolverMMP *solver,
                               int32_t scaling,
                               int32_t pct_inc_workspace,
                               int32_t max_work_memory,
-                              int32_t openmp_num_threads) {
+                              int32_t openmp_num_threads,
+                              int32_t compute_determinant) {
     if (solver == NULL) {
         return NULL_POINTER_ERROR;
     }
@@ -141,13 +143,25 @@ int32_t solver_mmp_initialize(struct SolverMMP *solver,
     solver->data.ICNTL(28) = MUMPS_ICNTL28_SEQUENTIAL;
     solver->data.ICNTL(29) = MUMPS_IGNORED;
 
+    if (compute_determinant == C_TRUE) {
+        // The determinant is obtained by computing
+        // (a + ib) * 2^c where a = RINFOG(12), b = RINFOG(13) and c = INFOG(34).
+        // In real arithmetic b = RINFOG(13) is equal to 0.
+        solver->data.ICNTL(33) = 1;
+        solver->data.ICNTL(8) = 0; // it's recommended to disable scaling when computing the determinant
+    } else {
+        solver->data.ICNTL(33) = 0;
+    }
+
     return 0; // success
 }
 
 int32_t solver_mmp_factorize(struct SolverMMP *solver,
                              int32_t const *indices_i,
                              int32_t const *indices_j,
                              double const *values_aij,
+                             double *determinant_coefficient_a,
+                             double *determinant_exponent_c,
                              int32_t verbose) {
     if (solver == NULL) {
         return NULL_POINTER_ERROR;
@@ -177,6 +191,16 @@ int32_t solver_mmp_factorize(struct SolverMMP *solver,
     solver->data.job = MUMPS_JOB_FACTORIZE;
     dmumps_c(&solver->data);
 
+    // read determinant
+
+    if (solver->data.ICNTL(33) == 1) {
+        *determinant_coefficient_a = solver->data.RINFOG(12);
+        *determinant_exponent_c = solver->data.INFOG(34);
+    } else {
+        *determinant_coefficient_a = 0.0;
+        *determinant_exponent_c = 0.0;
+    }
+
     return solver->data.INFOG(1);
 }
 

russell_sparse/data/python-validations/determinant-small-sparse-matrix.py

@@ -0,0 +1,14 @@
+import numpy as np
+import scipy.linalg as la
+
+a = np.array([
+    [ 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]
+], dtype=float)
+
+print(a)
+print(f"determinant(a) = {la.det(a)}")
+

russell_sparse/src/bin/solve_matrix_market_build.rs

@@ -37,6 +37,8 @@ pub struct SolutionInfo {
     pub verify_relative_error: f64,
     pub verify_time_nanosecond: u128,
     pub verify_time_human: String,
+    pub determinant_coefficient_a: f64, // det = a * 2^c
+    pub determinant_exponent_c: f64,    // det = a * 2^c
 }
 
 /// Command line options
@@ -68,6 +70,10 @@ struct Options {
     /// Activate verbose mode
     #[structopt(short = "v", long)]
     verbose: bool,
+
+    /// Computes determinant
+    #[structopt(short = "d", long)]
+    determinant: bool,
 }
 
 fn main() -> Result<(), StrError> {
@@ -114,6 +120,9 @@ fn main() -> Result<(), StrError> {
     if opt.verbose {
         config.verbose();
     }
+    if opt.determinant {
+        config.set_compute_determinant();
+    }
 
     // initialize and factorize
     let (nrow, nnz) = (coo.nrow, coo.pos);
@@ -142,6 +151,7 @@ fn main() -> Result<(), StrError> {
 
     // output
     let (time_fact, time_solve) = solver.get_elapsed_times();
+    let (det_a, det_c) = solver.get_determinant();
     let info = SolutionInfo {
         platform: "Russell".to_string(),
         blas_lib: "OpenBLAS".to_string(),
@@ -171,6 +181,8 @@ fn main() -> Result<(), StrError> {
         verify_relative_error: verify.relative_error,
         verify_time_nanosecond: verify.time_check,
         verify_time_human: format_nanoseconds(verify.time_check),
+        determinant_coefficient_a: det_a,
+        determinant_exponent_c: det_c,
     };
     let info_json = serde_json::to_string_pretty(&info).unwrap();
     println!("{}", info_json);

russell_sparse/src/config_solver.rs

@@ -11,6 +11,7 @@ pub struct ConfigSolver {
     pub(crate) max_work_memory: i32,     // max size of the working memory in mega bytes (MMP-only)
     pub(crate) openmp_num_threads: i32,  // number of OpenMP threads (MMP-only)
     pub(crate) verbose: i32,             // show lower-level messages
+    pub(crate) compute_determinant: i32, // compute determinant
 }
 
 impl ConfigSolver {
@@ -24,6 +25,7 @@ impl ConfigSolver {
             max_work_memory: 0,     // (MMP-only) 0 => Auto
             openmp_num_threads: 1,  // (MMP-only)
             verbose: 0,
+            compute_determinant: 0,
         }
     }
 
@@ -69,6 +71,12 @@ impl ConfigSolver {
         self
     }
 
+    /// Sets option to compute the determinant
+    pub fn set_compute_determinant(&mut self) -> &mut Self {
+        self.compute_determinant = 1;
+        self
+    }
+
     /// Returns a string representation of the ordering option
     pub fn str_ordering(&self) -> String {
         str_enum_ordering(self.ordering).to_string()
@@ -102,7 +110,7 @@ mod tests {
 
     #[test]
     fn clone_copy_and_debug_work() {
-        let correct = "ConfigSolver { lin_sol_kind: Umf, ordering: 2, scaling: 0, pct_inc_workspace: 100, max_work_memory: 0, openmp_num_threads: 1, verbose: 0 }";
+        let correct = "ConfigSolver { lin_sol_kind: Umf, ordering: 2, scaling: 0, pct_inc_workspace: 100, max_work_memory: 0, openmp_num_threads: 1, verbose: 0, compute_determinant: 0 }";
         let config = ConfigSolver::new();
         let copy = config;
         let clone = config.clone();

russell_sparse/src/solver.rs

@@ -26,12 +26,15 @@ extern "C" {
         pct_inc_workspace: i32,
         max_work_memory: i32,
         openmp_num_threads: i32,
+        compute_determinant: i32,
     ) -> i32;
     fn solver_mmp_factorize(
         solver: *mut ExtSolver,
         indices_i: *const i32,
         indices_j: *const i32,
         values_aij: *const f64,
+        determinant_coefficient_a: *mut f64,
+        determinant_exponent_c: *mut f64,
         verbose: i32,
     ) -> i32;
     fn solver_mmp_solve(solver: *mut ExtSolver, rhs: *mut f64, verbose: i32) -> i32;
@@ -72,16 +75,18 @@ extern "C" {
 /// (m,m)   (m)    (m)
 /// ```
 pub struct Solver {
-    kind: LinSolKind,            // solver kind
-    verbose: i32,                // verbose mode
-    done_factorize: bool,        // factorization completed
-    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
-    time_solve: u128,            // elapsed time during solve
-    used_ordering: &'static str, // used ordering strategy
-    used_scaling: &'static str,  // used scaling strategy
+    kind: LinSolKind,               // solver kind
+    verbose: i32,                   // verbose mode
+    done_factorize: bool,           // factorization completed
+    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
+    time_solve: u128,               // elapsed time during solve
+    used_ordering: &'static str,    // used ordering strategy
+    used_scaling: &'static str,     // used scaling strategy
+    determinant_coefficient_a: f64, // coefficient a of determinant = a * 2^c
+    determinant_exponent_c: f64,    // exponent c of determinant = a * 2^c
 }
 
 impl Solver {
@@ -109,6 +114,7 @@ impl Solver {
                         config.pct_inc_workspace,
                         config.max_work_memory,
                         config.openmp_num_threads,
+                        config.compute_determinant,
                     );
                     if res != 0 {
                         drop_solver_mmp(solver);
@@ -142,6 +148,8 @@ impl Solver {
                 time_solve: 0,
                 used_ordering: str_enum_ordering(config.ordering),
                 used_scaling: str_enum_scaling(config.scaling),
+                determinant_coefficient_a: 0.0,
+                determinant_exponent_c: 0.0,
             })
         }
     }
@@ -163,6 +171,8 @@ impl Solver {
                         coo.indices_i.as_ptr(),
                         coo.indices_j.as_ptr(),
                         coo.values_aij.as_ptr(),
+                        &mut self.determinant_coefficient_a,
+                        &mut self.determinant_exponent_c,
                         self.verbose,
                     );
                     if res != 0 {
@@ -389,6 +399,17 @@ impl Solver {
         self.used_scaling.to_string()
     }
 
+    /// Returns the coefficient and exponent of the determinant, if this option is requested
+    ///
+    /// Returns `(a, c)`, such that
+    ///
+    /// ```text
+    /// det(a) = a * 2^c
+    /// ```
+    pub fn get_determinant(&self) -> (f64, f64) {
+        (self.determinant_coefficient_a, self.determinant_exponent_c)
+    }
+
     /// Handles error code
     fn handle_mmp_error_code(err: i32) -> StrError {
         match err {
@@ -513,8 +534,8 @@ impl Drop for Solver {
 mod tests {
     use super::{ConfigSolver, CooMatrix, LinSolKind, Solver};
     use crate::Layout;
-    use russell_chk::vec_approx_eq;
-    use russell_lab::Vector;
+    use russell_chk::{approx_eq, vec_approx_eq};
+    use russell_lab::{mat_inverse, Matrix, Vector};
 
     #[test]
     fn new_works() {
@@ -559,6 +580,41 @@ mod tests {
         assert!(solver.done_factorize);
     }
 
+    #[test]
+    fn factorize_computes_the_determinant() {
+        let mut config = ConfigSolver::new();
+        config.lin_sol_kind(LinSolKind::Mmp);
+        config.set_compute_determinant();
+
+        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();
+
+        let mat = coo.as_matrix();
+        let mut inv = Matrix::new(nrow, ncol);
+        let det = mat_inverse(&mut inv, &mat).unwrap();
+        approx_eq(det, 114.0, 1e-15);
+
+        let mut solver = Solver::new(config, nrow, nnz, None).unwrap();
+        solver.factorize(&coo).unwrap();
+        let d = solver.determinant_coefficient_a * f64::powf(2.0, solver.determinant_exponent_c);
+        approx_eq(solver.determinant_coefficient_a, 57.0 / 64.0, 1e-15);
+        approx_eq(solver.determinant_exponent_c, 7.0, 1e-15);
+        approx_eq(d, 114.0, 1e-15);
+    }
+
     #[test]
     fn solve_fails_on_non_factorized() {
         let config = ConfigSolver::new();