[wip] Move SolverMUMPS to own file and structure

russell_sparse/c_code/interface_mumps.c

@@ -5,11 +5,6 @@
 #include "constants.h"
 #include "interface_mumps.h"
 
-#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_mumps_verbose(DMUMPS_STRUC_C *data, int32_t verbose) {
     if (verbose == C_TRUE) {
         data->ICNTL(1) = 6; // standard output stream
@@ -24,13 +19,8 @@ static inline void set_mumps_verbose(DMUMPS_STRUC_C *data, int32_t verbose) {
     }
 }
 
-struct SolverMUMPS {
-    DMUMPS_STRUC_C data;   // data structure
-    int32_t done_job_init; // job init successfully
-};
-
-struct SolverMUMPS *new_solver_mumps() {
-    struct SolverMUMPS *solver = (struct SolverMUMPS *)malloc(sizeof(struct SolverMUMPS));
+struct InterfaceMUMPS *solver_mumps_new() {
+    struct InterfaceMUMPS *solver = (struct InterfaceMUMPS *)malloc(sizeof(struct InterfaceMUMPS));
 
     if (solver == NULL) {
         return NULL;
@@ -44,7 +34,7 @@ struct SolverMUMPS *new_solver_mumps() {
     return solver;
 }
 
-void drop_solver_mumps(struct SolverMUMPS *solver) {
+void solver_mumps_drop(struct InterfaceMUMPS *solver) {
     if (solver == NULL) {
         return;
     }
@@ -71,23 +61,23 @@ void drop_solver_mumps(struct SolverMUMPS *solver) {
     free(solver);
 }
 
-int32_t solver_mumps_initialize(struct SolverMUMPS *solver,
-                              int32_t n,
-                              int32_t nnz,
-                              int32_t symmetry,
-                              int32_t ordering,
-                              int32_t scaling,
-                              int32_t pct_inc_workspace,
-                              int32_t max_work_memory,
-                              int32_t openmp_num_threads,
-                              int32_t compute_determinant) {
+int32_t solver_mumps_initialize(struct InterfaceMUMPS *solver,
+                                int32_t n,
+                                int32_t nnz,
+                                int32_t symmetry,
+                                int32_t ordering,
+                                int32_t scaling,
+                                int32_t pct_inc_workspace,
+                                int32_t max_work_memory,
+                                int32_t openmp_num_threads,
+                                int32_t compute_determinant) {
     if (solver == NULL) {
         return NULL_POINTER_ERROR;
     }
 
     solver->data.comm_fortran = MUMPS_IGNORED;
     solver->data.par = MUMPS_PAR_HOST_ALSO_WORKS;
-    solver->data.sym = MUMPS_SYMMETRY[symmetry];
+    solver->data.sym = symmetry;
 
     set_mumps_verbose(&solver->data, C_FALSE);
     solver->data.job = MUMPS_JOB_INITIALIZE;
@@ -130,8 +120,8 @@ int32_t solver_mumps_initialize(struct SolverMUMPS *solver,
 
     solver->data.ICNTL(5) = MUMPS_ICNTL5_ASSEMBLED_MATRIX;
     solver->data.ICNTL(6) = MUMPS_ICNTL6_PERMUT_AUTO;
-    solver->data.ICNTL(7) = MUMPS_ORDERING[ordering];
-    solver->data.ICNTL(8) = MUMPS_SCALING[scaling];
+    solver->data.ICNTL(7) = ordering;
+    solver->data.ICNTL(8) = scaling;
     solver->data.ICNTL(14) = pct_inc_workspace;
     solver->data.ICNTL(16) = openmp_num_threads;
     solver->data.ICNTL(18) = MUMPS_ICNTL18_CENTRALIZED;
@@ -152,13 +142,11 @@ int32_t solver_mumps_initialize(struct SolverMUMPS *solver,
     return 0; // success
 }
 
-int32_t solver_mumps_factorize(struct SolverMUMPS *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) {
+int32_t solver_mumps_factorize(struct InterfaceMUMPS *solver,
+                               int32_t const *indices_i,
+                               int32_t const *indices_j,
+                               double const *values_aij,
+                               int32_t verbose) {
     if (solver == NULL) {
         return NULL_POINTER_ERROR;
     }
@@ -190,17 +178,17 @@ int32_t solver_mumps_factorize(struct SolverMUMPS *solver,
     // read determinant
 
     if (solver->data.ICNTL(33) == 1) {
-        *determinant_coefficient_a = solver->data.RINFOG(12);
-        *determinant_exponent_c = solver->data.INFOG(34);
+        solver->determinant_coefficient_a = solver->data.RINFOG(12);
+        solver->determinant_exponent_c = solver->data.INFOG(34);
     } else {
-        *determinant_coefficient_a = 0.0;
-        *determinant_exponent_c = 0.0;
+        solver->determinant_coefficient_a = 0.0;
+        solver->determinant_exponent_c = 0.0;
     }
 
     return solver->data.INFOG(1);
 }
 
-int32_t solver_mumps_solve(struct SolverMUMPS *solver, double *rhs, int32_t verbose) {
+int32_t solver_mumps_solve(struct InterfaceMUMPS *solver, double *rhs, int32_t verbose) {
     if (solver == NULL) {
         return NULL_POINTER_ERROR;
     }
@@ -214,10 +202,30 @@ int32_t solver_mumps_solve(struct SolverMUMPS *solver, double *rhs, int32_t verb
     return solver->data.INFOG(1);
 }
 
-int32_t solver_mumps_used_ordering(struct SolverMUMPS *solver) {
+int32_t solver_mumps_get_ordering(const struct InterfaceMUMPS *solver) {
+    if (solver == NULL) {
+        return 0;
+    }
     return solver->data.INFOG(7);
 }
 
-int32_t solver_mumps_used_scaling(struct SolverMUMPS *solver) {
+int32_t solver_mumps_get_scaling(const struct InterfaceMUMPS *solver) {
+    if (solver == NULL) {
+        return 0;
+    }
     return solver->data.INFOG(33);
 }
+
+double solver_mumps_get_det_coef_a(const struct InterfaceMUMPS *solver) {
+    if (solver == NULL) {
+        return 0.0;
+    }
+    return solver->determinant_coefficient_a;
+}
+
+double solver_mumps_get_det_exp_c(const struct InterfaceMUMPS *solver) {
+    if (solver == NULL) {
+        return 0.0;
+    }
+    return solver->determinant_exponent_c;
+}

russell_sparse/c_code/interface_mumps.h

@@ -1,49 +1,100 @@
 #pragma once
 
 #include "dmumps_c.h"
+#include "stdlib.h"
+
+#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
 
 const MUMPS_INT MUMPS_IGNORED = 0;
 
-const MUMPS_INT MUMPS_JOB_INITIALIZE = -1;
-const MUMPS_INT MUMPS_JOB_TERMINATE = -2;
-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_JOB_INITIALIZE = -1; // section 5.1.1, page 24
+const MUMPS_INT MUMPS_JOB_TERMINATE = -2;  // section 5.1.1, page 24
+const MUMPS_INT MUMPS_JOB_ANALYZE = 1;     // section 5.1.1, page 24
+const MUMPS_INT MUMPS_JOB_FACTORIZE = 2;   // section 5.1.1, page 25
+const MUMPS_INT MUMPS_JOB_SOLVE = 3;       // section 5.1.1, page 25
 
 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_SYMMETRY[3] = {
-    0, // Unsymmetric
-    1, // Positive-definite symmetric
-    2, // General symmetric
+/// @brief Holds the data for MUMPS
+struct InterfaceMUMPS {
+    DMUMPS_STRUC_C data;              // data structure
+    int32_t done_job_init;            // job init successfully
+    double determinant_coefficient_a; // if asked, stores the determinant coefficient a
+    double determinant_exponent_c;    // if asked, stores the determinant exponent c
 };
 
-const MUMPS_INT MUMPS_ORDERING[10] = {
-    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
-};
+/// @brief Allocates a new MUMPS interface
+struct InterfaceMUMPS *solver_mumps_new();
 
-const MUMPS_INT MUMPS_SCALING[9] = {
-    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
-};
+/// @brief Deallocates the MUMPS interface
+void solver_mumps_drop(struct InterfaceMUMPS *solver);
+
+/// @brief Performs the initialization
+/// @param solver Is a pointer to the solver interface
+/// @param n Is the number of rows and columns of the coefficient matrix
+/// @param nnz Is the number of non-zero values in the coefficient matrix
+/// @param symmetry Is the code for the kind of symmetry, if any
+/// @param ordering Is the russel_sparse ordering code
+/// @param scaling Is the russell_sparse scaling code
+/// @param pct_inc_workspace Is the allowed percentage increase of the workspace
+/// @param max_work_memory Is the allowed maximum memory
+/// @param openmp_num_threads Is the number of threads allowed for OpenMP
+/// @param compute_determinant Requests that determinant be computed
+/// @return A successful or error code
+int32_t solver_mumps_initialize(struct InterfaceMUMPS *solver,
+                                int32_t n,
+                                int32_t nnz,
+                                int32_t symmetry,
+                                int32_t ordering,
+                                int32_t scaling,
+                                int32_t pct_inc_workspace,
+                                int32_t max_work_memory,
+                                int32_t openmp_num_threads,
+                                int32_t compute_determinant);
+
+/// @brief Performs the factorization
+/// @param solver Is a pointer to the solver interface
+/// @param indices_i Are the CooMatrix row indices
+/// @param indices_j Are the CooMatrix column indices
+/// @param values_aij Are the CooMatrix values
+/// @param verbose Shows messages
+/// @return A successful or error code
+int32_t solver_mumps_factorize(struct InterfaceMUMPS *solver,
+                               int32_t const *indices_i,
+                               int32_t const *indices_j,
+                               double const *values_aij,
+                               int32_t verbose);
+
+/// @brief Computes the solution of the linear system
+/// @param solver Is a pointer to the solver interface
+/// @param rhs Is the right-hand side on the input and the vector of unknow values x on the output
+/// @param verbose Shows messages
+/// @return A successful or error code
+int32_t solver_mumps_solve(struct InterfaceMUMPS *solver, double *rhs, int32_t verbose);
+
+/// @brief Returns the effective ordering using during the computations
+/// @param solver Is a pointer to the solver
+/// @return The MUMPS ordering code
+int32_t solver_mumps_get_ordering(const struct InterfaceMUMPS *solver);
+
+/// @brief Returns the effective scaling using during the computations
+/// @param solver Is a pointer to the solver
+/// @return A success or failure code
+int32_t solver_mumps_get_scaling(const struct InterfaceMUMPS *solver);
+
+/// @brief Returns the coefficient needed to compute the determinant, if requested
+/// @param solver Is a pointer to the solver
+/// @return Coefficient a of a * 2^c
+double solver_mumps_get_det_coef_a(const struct InterfaceMUMPS *solver);
+
+/// @brief Returns the exponent needed to compute the determinant, if requested
+/// @param solver Is a pointer to the solver
+/// @return Exponent c of a * 2^c
+double solver_mumps_get_det_exp_c(const struct InterfaceMUMPS *solver);

russell_sparse/src/auxiliary.rs

@@ -0,0 +1,7 @@
+use crate::StrError;
+
+/// Converts usize to i32
+#[inline]
+pub(crate) fn to_i32(num: usize) -> Result<i32, StrError> {
+    i32::try_from(num).map_err(|_| "cannot downcast usize to i32")
+}

russell_sparse/src/bin/mem_check_build.rs

@@ -3,7 +3,6 @@ use russell_sparse::prelude::*;
 
 fn test_solver(name: LinSolKind) {
     match name {
-        LinSolKind::Mumps => println!("Testing MUMPS solver\n"),
         LinSolKind::Umfpack => println!("Testing UMFPACK solver\n"),
     }
 
@@ -73,7 +72,6 @@ fn test_solver(name: LinSolKind) {
 
 fn test_solver_singular(name: LinSolKind) {
     match name {
-        LinSolKind::Mumps => println!("Testing MUMPS solver\n"),
         LinSolKind::Umfpack => println!("Testing UMFPACK solver\n"),
     }
 
@@ -105,9 +103,7 @@ fn test_solver_singular(name: LinSolKind) {
 
 fn main() {
     println!("Running Mem Check\n");
-    test_solver(LinSolKind::Mumps);
     test_solver(LinSolKind::Umfpack);
-    test_solver_singular(LinSolKind::Mumps);
     test_solver_singular(LinSolKind::Umfpack);
     println!("Done\n");
 }

russell_sparse/src/bin/solve_matrix_market_build.rs

@@ -86,14 +86,14 @@ fn main() -> Result<(), StrError> {
     }
 
     // select linear solver
-    let name = if opt.mumps { LinSolKind::Mumps } else { LinSolKind::Umfpack };
+    let name = LinSolKind::Umfpack;
 
     // select the symmetric handling option
     let handling = match name {
-        LinSolKind::Mumps => {
-            // MUMPS uses the lower-diagonal if symmetric.
-            SymmetricHandling::LeaveAsLower
-        }
+        // LinSolKind::Mumps => {
+        //     // MUMPS uses the lower-diagonal if symmetric.
+        //     SymmetricHandling::LeaveAsLower
+        // }
         LinSolKind::Umfpack => {
             // UMFPACK uses the full matrix, if symmetric or not
             SymmetricHandling::MakeItFull
@@ -151,7 +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 (det_a, det_c) = (0.0, 0.0); //solver.get_determinant());
     let info = SolutionInfo {
         platform: "Russell".to_string(),
         blas_lib: "OpenBLAS".to_string(),

russell_sparse/src/config_solver.rs

@@ -90,13 +90,6 @@ impl ConfigSolver {
     /// Returns the name of the solver
     pub fn str_solver(&self) -> String {
         match self.lin_sol_kind {
-            LinSolKind::Mumps => {
-                if cfg!(local_mumps) {
-                    "MUMPS-local".to_string()
-                } else {
-                    "MUMPS".to_string()
-                }
-            }
             LinSolKind::Umfpack => "UMFPACK".to_string(),
         }
     }
@@ -106,7 +99,7 @@ impl ConfigSolver {
 
 #[cfg(test)]
 mod tests {
-    use super::{ConfigSolver, LinSolKind, Ordering, Scaling};
+    use super::{ConfigSolver, Ordering, Scaling};
 
     #[test]
     fn clone_copy_and_debug_work() {
@@ -130,18 +123,6 @@ mod tests {
         assert_eq!(config.verbose, 0);
     }
 
-    #[test]
-    fn set_solver_works() {
-        let mut config = ConfigSolver::new();
-        for name in [LinSolKind::Mumps, LinSolKind::Umfpack] {
-            config.lin_sol_kind(name);
-            match config.lin_sol_kind {
-                LinSolKind::Mumps => assert!(true),
-                LinSolKind::Umfpack => assert!(true),
-            }
-        }
-    }
-
     #[test]
     fn set_ordering_works() {
         let mut config = ConfigSolver::new();