CHOLMOD Wrapper Integration

src/CMakeLists.txt

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

src/polysolve/CHOLMODSolver.cpp

@@ -0,0 +1,118 @@
+#ifdef POLYSOLVE_WITH_CHOLMOD
+
+#include "CHOLMODSolver.hpp"
+#include <Eigen/CholmodSupport>
+
+#include <stdio.h>
+#include <iostream>
+#include <memory>
+
+namespace polysolve
+{
+    namespace
+    {
+        cholmod_sparse eigen2cholmod(const StiffnessMatrix &mat)
+        {
+            cholmod_sparse res;
+            res.nzmax = mat.nonZeros();
+            res.nrow = mat.rows();
+            res.ncol = mat.cols();
+
+            memcpy(res.p, mat.outerIndexPtr(), sizeof(mat.outerIndexPtr()));
+            memcpy(res.i, mat.innerIndexPtr(), sizeof(mat.innerIndexPtr()));
+            memcpy(res.x, mat.valuePtr(), sizeof(mat.valuePtr()));
+
+            res.z = 0;
+            res.sorted = 1;
+            // if (mat.isCompressed())
+            {
+                assert(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;
+        }
+    } // namespace
+
+    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(const StiffnessMatrix &Ain, const int precond_num)
+    {
+        assert(Ain.isCompressed());
+
+        A = eigen2cholmod(Ain);
+        L = cholmod_l_analyze(&A, cm);
+    }
+
+    void CHOLMODSolver::factorize(const StiffnessMatrix &Ain)
+    {
+        cholmod_l_factorize(&A, L, cm);
+    }
+
+    void CHOLMODSolver::solve(const Ref<const VectorXd> rhs, Ref<VectorXd> result)
+    {
+        assert(result.size() == rhs.size());
+
+        // b = eigen2cholmod(rhs); //TODO: fix me
+
+        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
+
+#endif

src/polysolve/CHOLMODSolver.hpp

@@ -0,0 +1,65 @@
+#pragma once
+
+#ifdef POLYSOLVE_WITH_CHOLMOD
+
+////////////////////////////////////////////////////////////////////////////////
+
+#include <polysolve/LinearSolver.hpp>
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+#include <cholmod.h>
+#include <memory>
+
+#include <nlohmann/json.hpp>
+using json = nlohmann::json;
+
+namespace polysolve
+{
+    /**
+     @brief      Base class for cholmod solver.
+     */
+    class CHOLMODSolver : public LinearSolver
+    {
+
+    protected:
+        cholmod_common *cm;
+        cholmod_sparse A;
+        cholmod_dense *x, b;
+        cholmod_factor *L;
+
+    public:
+        CHOLMODSolver();
+        ~CHOLMODSolver();
+
+    private:
+        POLYSOLVE_DELETE_MOVE_COPY(CHOLMODSolver)
+    public:
+        // Set solver parameters
+        // virtual void setParameters(const json &params) {}
+
+        // Get info on the last solve step
+        virtual void getInfo(json &params) const override;
+
+        // Analyze sparsity pattern
+        virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override;
+
+        // Factorize system matrix
+        void factorize(const StiffnessMatrix &A) override;
+
+        //
+        // @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. }
+        //
+        virtual void solve(const Ref<const VectorXd> b, Ref<VectorXd> x) override;
+    };
+
+} // namespace polysolve
+
+#endif

src/polysolve/LinearSolver.cpp

@@ -7,6 +7,7 @@
 #include <Eigen/Sparse>
 #ifdef POLYSOLVE_WITH_CHOLMOD
 #include <Eigen/CholmodSupport>
+#include <polysolve/CHOLMODSolver.hpp>
 #endif
 #ifdef POLYSOLVE_WITH_UMFPACK
 #include <Eigen/UmfPackSupport>
@@ -183,6 +184,10 @@ namespace polysolve
         else if (solver == "Eigen::CholmodSupernodalLLT")
         {
             RETURN_DIRECT_SOLVER_PTR(CholmodSupernodalLLT, "Eigen::CholmodSupernodalLLT");
+        }
+        else if (solver == "CholmodGPU")
+        {
+            std::make_unique<CHOLMODSolver>();
 #endif
 #ifdef POLYSOLVE_WITH_UMFPACK
 #ifndef POLYSOLVE_LARGE_INDEX
@@ -275,6 +280,7 @@ namespace polysolve
             "Eigen::SparseLU",
 #ifdef POLYSOLVE_WITH_CHOLMOD
             "Eigen::CholmodSupernodalLLT",
+            "CholmodGPU",
 #endif
 #ifdef POLYSOLVE_WITH_UMFPACK
             "Eigen::UmfPackLU",

tests/test_solver.cpp

@@ -459,7 +459,7 @@ TEST_CASE("amgcl_blocksolver_crystm03_CG", "[solver]")
         json solver_info;
         auto solver = LinearSolver::create("AMGCL", "");
         prof.tic("setup");
-                json params;
+        json params;
         params["AMGCL"]["tolerance"] = 1e-8;
         params["AMGCL"]["max_iter"] = 1000;
         params["AMGCL"]["block_size"] = 3;
@@ -481,7 +481,7 @@ TEST_CASE("amgcl_blocksolver_crystm03_CG", "[solver]")
         json solver_info;
         auto solver = LinearSolver::create("AMGCL", "");
         prof.tic("setup");
-                json params;
+        json params;
         params["AMGCL"]["tolerance"] = 1e-8;
         params["AMGCL"]["max_iter"] = 10000;
         solver->setParameters(params);
@@ -571,3 +571,26 @@ TEST_CASE("amgcl_blocksolver_crystm03_Bicgstab", "[solver]")
     REQUIRE((A * x_b - b).norm() / b.norm() < 1e-7);
 }
 #endif
+
+#ifdef POLYSOLVE_WITH_CHOLMOD
+TEST_CASE("CHOLMOD", "[solver]")
+{
+    const std::string path = POLYSOLVE_DATA_DIR;
+    Eigen::SparseMatrix<double, Eigen::ColMajor, long int> A;
+    bool ok = loadMarket(A, path + "/A_2.mat");
+    REQUIRE(ok);
+
+    Eigen::VectorXd b(A.rows());
+    b.setOnes();
+    Eigen::VectorXd x(A.rows());
+    x.setZero();
+
+    auto solver = LinearSolver::create("CholmodGPU", "");
+    solver->analyzePattern(A, A.rows());
+    solver->factorize(A);
+    solver->solve(b, x);
+
+    const double err = (b - (A * x)).norm();
+    REQUIRE(err < 1e-8);
+}
+#endif