CHOLMOD Wrapper Integration

src/CHOLMODSolver.cpp

@@ -0,0 +1,104 @@
+#include "CHOLMODSolver.hpp"
+#include <Eigen/CholmodSupport>
+#include <iostream>
+
+namespace polysolve
+{
+
+    cholmod_sparse eigen2cholmod(StiffnessMatrixL &mat)
+    {
+        cholmod_sparse res;
+        res.nzmax   = mat.nonZeros();
+        res.nrow    = mat.rows();
+        res.ncol    = mat.cols();
+        res.p       = mat.outerIndexPtr();
+        res.i       = mat.innerIndexPtr();
+        res.x       = mat.valuePtr();
+        res.z       = 0;
+        res.sorted  = 1;
+        if(mat.isCompressed())
+        {
+            res.packed  = 1;
+            res.nz = 0;
+        }
+        else
+        {
+            res.packed  = 0;
+            res.nz = mat.innerNonZeroPtr();
+        }
+
+        res.dtype = CHOLMOD_DOUBLE;
+
+        res.itype = CHOLMOD_LONG;
+        res.stype = 1;
+        res.xtype = CHOLMOD_REAL;
+
+        return res;
+    }
+
+    cholmod_dense eigen2cholmod(Eigen::VectorXd &mat)
+    {
+        cholmod_dense res;
+        res.nrow   = mat.size();
+        res.ncol   = 1;
+        res.nzmax  = res.nrow * res.ncol;
+        res.d      = mat.derived().size();
+        res.x      = (void*)(mat.derived().data());
+        res.z      = 0;
+        res.xtype = CHOLMOD_REAL;
+        res.dtype   = 0;
+
+        return res;
+    }
+    ////////////////////////////////////////////////////////////////////////////////
+
+    CHOLMODSolver::CHOLMODSolver()
+    {
+        cm = (cholmod_common*)malloc(sizeof(cholmod_common));
+        cholmod_l_start (cm);
+        L = NULL;
+
+        cm->useGPU = 1;
+        cm->maxGpuMemBytes = 2e9;
+        cm->print = 4;
+        cm->supernodal = CHOLMOD_SUPERNODAL;
+    }
+
+    // void CHOLMODSolver::getInfo(json &params) const
+    // {
+    //     params["num_iterations"] = num_iterations;
+    //     params["final_res_norm"] = final_res_norm;
+    // }
+
+    void CHOLMODSolver::analyzePattern(StiffnessMatrixL &Ain) 
+    {   
+        if(!Ain.isCompressed()) {
+            Ain.makeCompressed();
+        }
+        A = eigen2cholmod(Ain);
+        L = cholmod_l_analyze (&A, cm);
+    }
+
+    void CHOLMODSolver::factorize(StiffnessMatrixL &Ain)
+    {
+        cholmod_l_factorize (&A, L, cm);
+    }
+
+    void CHOLMODSolver::solve(Eigen::VectorXd &rhs, Eigen::VectorXd &result)
+    {
+        b = eigen2cholmod(rhs);
+        x = cholmod_l_solve (CHOLMOD_A, L, &b, cm);
+
+        memcpy(result.data(), x->x, result.size() * sizeof(result[0]));
+    }
+
+    ////////////////////////////////////////////////////////////////////////////////
+
+    CHOLMODSolver::~CHOLMODSolver()
+    {
+        cholmod_l_gpu_stats(cm);
+        cholmod_l_free_factor (&L, cm);
+        cholmod_l_free_dense (&x, cm);
+    }
+
+} // namespace polysolve

src/CHOLMODSolver.hpp

@@ -0,0 +1,100 @@
+#pragma once
+
+////////////////////////////////////////////////////////////////////////////////
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include "cholmod.h"
+#include <memory>
+
+#include <nlohmann/json.hpp>
+using json = nlohmann::json;
+
+// #define POLYSOLVE_DELETE_MOVE_COPY(Base) \
+//     Base(Base &&) = delete;                    \
+//     Base &operator=(Base &&) = delete;         \
+//     Base(const Base &) = delete;               \
+//     Base &operator=(const Base &) = delete;
+
+////////////////////////////////////////////////////////////////////////////////
+// TODO:
+// - [ ] Support both RowMajor + ColumnMajor sparse matrices
+// - [ ] Wrapper around MUMPS
+// - [ ] Wrapper around other iterative solvers (AMGCL, ViennaCL, etc.)
+// - [ ] Document the json parameters for each
+////////////////////////////////////////////////////////////////////////////////
+
+namespace polysolve
+{
+    typedef Eigen::SparseMatrix<double, Eigen::ColMajor, long int> StiffnessMatrixL;
+    /**
+     @brief      Base class for cholmod solver.
+ */
+    class CHOLMODSolver
+    {
+
+    // public:
+        // Shortcut alias
+        // typedef Eigen::VectorXd VectorXd;
+        // template <typename T>
+        // using Ref = Eigen::Ref<T>;
+
+    // public:
+        //////////////////
+        // Constructors //
+        //////////////////
+
+        // Virtual destructor
+
+
+        // Static constructor
+        //
+        // @param[in]  solver   Solver type
+        // @param[in]  precond  Preconditioner for iterative solvers
+        //
+        // static std::unique_ptr<LinearSolver> create(const std::string &solver, const std::string &precond);
+
+        // List available solvers
+        // static std::vector<std::string> availableSolvers();
+        // static std::string defaultSolver();
+
+        // List available preconditioners
+        // static std::vector<std::string> availablePrecond();
+        // static std::string defaultPrecond();
+
+    protected:
+
+        cholmod_common *cm;
+        cholmod_sparse A;
+        cholmod_dense *x, b;
+        cholmod_factor *L;
+
+    public:
+        CHOLMODSolver();
+        ~CHOLMODSolver();
+
+        // Set solver parameters
+        // virtual void setParameters(const json &params) {}
+
+        // Get info on the last solve step
+        void getInfo(json &params) const;
+
+        // Analyze sparsity pattern
+        void analyzePattern(StiffnessMatrixL &Ain);
+
+        // Factorize system matrix
+        void factorize(StiffnessMatrixL &Ain);
+
+        //
+        // @brief         { Solve the linear system Ax = b }
+        //
+        // @param[in]     b     { Right-hand side. }
+        // @param[in,out] x     { Unknown to compute. When using an iterative
+        //                      solver, the input unknown vector is used as an
+        //                      initial guess, and must thus be properly allocated
+        //                      and initialized. }
+        //
+        void solve(Eigen::VectorXd &rhs, Eigen::VectorXd &result);
+    };
+
+} // namespace polysolve

src/CMakeLists.txt

@@ -14,6 +14,8 @@ set(SOURCES
     LinearSolverPardiso.hpp
     SaddlePointSolver.cpp
     SaddlePointSolver.hpp
+    CHOLMODSolver.cpp
+    CHOLMODSolver.hpp
 )
 
 polysolve_prepend_current_path(SOURCES)

tests/test_solver.cpp

@@ -1,5 +1,6 @@
 ////////////////////////////////////////////////////////////////////////////////
 #include <polysolve/FEMSolver.hpp>
+#include <polysolve/CHOLMODSolver.hpp>
 
 #include <catch2/catch.hpp>
 #include <iostream>
@@ -245,3 +246,24 @@ TEST_CASE("saddle_point_test", "[solver]")
     const double err = (A * x - b).norm();
     REQUIRE(err < 1e-8);
 }
+
+TEST_CASE("CHOLMOD", "[solver]")
+{
+    const std::string path = POLYSOLVE_DATA_DIR;
+    Eigen::SparseMatrix<double, Eigen::ColMajor, long int> A;
+    bool ok = loadMarket(A, path + "/nd6k.mtx");
+    REQUIRE(ok);
+
+    Eigen::VectorXd b(A.rows());
+    b.setOnes();
+    Eigen::VectorXd x(A.rows());
+    x.setZero();
+
+    CHOLMODSolver solver;
+    solver.analyzePattern(A);
+    solver.factorize(A);
+    solver.solve(b, x);
+
+    const double err = (b - (A * x)).norm();
+    REQUIRE(err < 1e-8);
+}