From 44d84abae0c4a3dfef79e9eca93bdef11a5042cb Mon Sep 17 00:00:00 2001 From: Pranav Jain Date: Thu, 10 Mar 2022 13:43:25 -0500 Subject: [PATCH 1/4] Added CHOLMOD Wrapper and added cholmod test --- src/CHOLMODSolver.cpp | 104 ++++++++++++++++++++++++++++++++++++++++++ src/CHOLMODSolver.hpp | 100 ++++++++++++++++++++++++++++++++++++++++ src/CMakeLists.txt | 2 + tests/test_solver.cpp | 26 +++++++++++ 4 files changed, 232 insertions(+) create mode 100644 src/CHOLMODSolver.cpp create mode 100644 src/CHOLMODSolver.hpp diff --git a/src/CHOLMODSolver.cpp b/src/CHOLMODSolver.cpp new file mode 100644 index 00000000..053f156c --- /dev/null +++ b/src/CHOLMODSolver.cpp @@ -0,0 +1,104 @@ +#include "CHOLMODSolver.hpp" +#include +#include + +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 ¶ms) 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 diff --git a/src/CHOLMODSolver.hpp b/src/CHOLMODSolver.hpp new file mode 100644 index 00000000..d7de5b5f --- /dev/null +++ b/src/CHOLMODSolver.hpp @@ -0,0 +1,100 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include "cholmod.h" +#include + +#include +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 StiffnessMatrixL; + /** + @brief Base class for cholmod solver. + */ + class CHOLMODSolver + { + + // public: + // Shortcut alias + // typedef Eigen::VectorXd VectorXd; + // template + // using Ref = Eigen::Ref; + + // public: + ////////////////// + // Constructors // + ////////////////// + + // Virtual destructor + + + // Static constructor + // + // @param[in] solver Solver type + // @param[in] precond Preconditioner for iterative solvers + // + // static std::unique_ptr create(const std::string &solver, const std::string &precond); + + // List available solvers + // static std::vector availableSolvers(); + // static std::string defaultSolver(); + + // List available preconditioners + // static std::vector 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 ¶ms) {} + + // Get info on the last solve step + void getInfo(json ¶ms) 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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ac7123a7..3f64fcf0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,6 +14,8 @@ set(SOURCES LinearSolverPardiso.hpp SaddlePointSolver.cpp SaddlePointSolver.hpp + CHOLMODSolver.cpp + CHOLMODSolver.hpp ) polysolve_prepend_current_path(SOURCES) diff --git a/tests/test_solver.cpp b/tests/test_solver.cpp index a554a1d4..e35e35e6 100644 --- a/tests/test_solver.cpp +++ b/tests/test_solver.cpp @@ -1,5 +1,6 @@ //////////////////////////////////////////////////////////////////////////////// #include +#include #include #include @@ -245,3 +246,28 @@ 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 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; + std::cout<<"here1\n"; + solver.analyzePattern(A); + std::cout<<"here2\n"; + solver.factorize(A); + std::cout<<"here3\n"; + solver.solve(b, x); + std::cout<<"here4\n"; + + const double err = (b - (A * x)).norm(); + REQUIRE(err < 1e-8); +} From 361d44125b844b8f50e9d2edd7a24b9618d48531 Mon Sep 17 00:00:00 2001 From: Pranav Jain Date: Thu, 10 Mar 2022 13:55:06 -0500 Subject: [PATCH 2/4] Cleaned the code --- tests/test_solver.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/test_solver.cpp b/tests/test_solver.cpp index e35e35e6..e7f73572 100644 --- a/tests/test_solver.cpp +++ b/tests/test_solver.cpp @@ -260,13 +260,9 @@ TEST_CASE("CHOLMOD", "[solver]") x.setZero(); CHOLMODSolver solver; - std::cout<<"here1\n"; solver.analyzePattern(A); - std::cout<<"here2\n"; solver.factorize(A); - std::cout<<"here3\n"; solver.solve(b, x); - std::cout<<"here4\n"; const double err = (b - (A * x)).norm(); REQUIRE(err < 1e-8); From aef3adf259c695bfce4c8dcc45d2ed4caa2f3fa2 Mon Sep 17 00:00:00 2001 From: Pranav Jain Date: Thu, 10 Mar 2022 14:15:52 -0500 Subject: [PATCH 3/4] Cleaned code and resized result --- src/CHOLMODSolver.cpp | 4 ++-- src/CHOLMODSolver.hpp | 44 +------------------------------------------ tests/test_solver.cpp | 2 +- 3 files changed, 4 insertions(+), 46 deletions(-) diff --git a/src/CHOLMODSolver.cpp b/src/CHOLMODSolver.cpp index 053f156c..c1fc0bd8 100644 --- a/src/CHOLMODSolver.cpp +++ b/src/CHOLMODSolver.cpp @@ -88,7 +88,7 @@ namespace polysolve { b = eigen2cholmod(rhs); x = cholmod_l_solve (CHOLMOD_A, L, &b, cm); - + result.conservativeResize(rhs.size()); memcpy(result.data(), x->x, result.size() * sizeof(result[0])); } @@ -96,7 +96,7 @@ namespace polysolve CHOLMODSolver::~CHOLMODSolver() { - cholmod_l_gpu_stats(cm); + // cholmod_l_gpu_stats(cm); cholmod_l_free_factor (&L, cm); cholmod_l_free_dense (&x, cm); } diff --git a/src/CHOLMODSolver.hpp b/src/CHOLMODSolver.hpp index d7de5b5f..16d6fa4c 100644 --- a/src/CHOLMODSolver.hpp +++ b/src/CHOLMODSolver.hpp @@ -10,19 +10,6 @@ #include 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 { @@ -32,36 +19,7 @@ namespace polysolve */ class CHOLMODSolver { - - // public: - // Shortcut alias - // typedef Eigen::VectorXd VectorXd; - // template - // using Ref = Eigen::Ref; - - // public: - ////////////////// - // Constructors // - ////////////////// - - // Virtual destructor - - - // Static constructor - // - // @param[in] solver Solver type - // @param[in] precond Preconditioner for iterative solvers - // - // static std::unique_ptr create(const std::string &solver, const std::string &precond); - - // List available solvers - // static std::vector availableSolvers(); - // static std::string defaultSolver(); - - // List available preconditioners - // static std::vector availablePrecond(); - // static std::string defaultPrecond(); - + protected: cholmod_common *cm; diff --git a/tests/test_solver.cpp b/tests/test_solver.cpp index e7f73572..d1213829 100644 --- a/tests/test_solver.cpp +++ b/tests/test_solver.cpp @@ -251,7 +251,7 @@ TEST_CASE("CHOLMOD", "[solver]") { const std::string path = POLYSOLVE_DATA_DIR; Eigen::SparseMatrix A; - bool ok = loadMarket(A, path + "/nd6k.mtx"); + bool ok = loadMarket(A, path + "/A_2.mat"); REQUIRE(ok); Eigen::VectorXd b(A.rows()); From ef98a7a3c9ca9c3852d3113176941c6b902e3002 Mon Sep 17 00:00:00 2001 From: Teseo Schneider Date: Tue, 29 Nov 2022 18:39:42 -0800 Subject: [PATCH 4/4] cleanup --- src/polysolve/CHOLMODSolver.cpp | 83 ++++++++++++++++++--------------- src/polysolve/CHOLMODSolver.hpp | 18 ++++--- src/polysolve/LinearSolver.cpp | 6 +++ tests/test_solver.cpp | 8 ++-- 4 files changed, 67 insertions(+), 48 deletions(-) diff --git a/src/polysolve/CHOLMODSolver.cpp b/src/polysolve/CHOLMODSolver.cpp index 4042aec9..b2396398 100644 --- a/src/polysolve/CHOLMODSolver.cpp +++ b/src/polysolve/CHOLMODSolver.cpp @@ -2,41 +2,49 @@ #include "CHOLMODSolver.hpp" #include + +#include #include +#include namespace polysolve { - - cholmod_sparse eigen2cholmod(StiffnessMatrixL &mat) + namespace { - 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 + cholmod_sparse eigen2cholmod(const StiffnessMatrix &mat) { - res.packed = 0; - res.nz = mat.innerNonZeroPtr(); + 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; } - - res.dtype = CHOLMOD_DOUBLE; - - res.itype = CHOLMOD_LONG; - res.stype = 1; - res.xtype = CHOLMOD_REAL; - - return res; - } + } // namespace cholmod_dense eigen2cholmod(Eigen::VectorXd &mat) { @@ -72,26 +80,27 @@ namespace polysolve // params["final_res_norm"] = final_res_norm; // } - void CHOLMODSolver::analyzePattern(StiffnessMatrixL &Ain) + void CHOLMODSolver::analyzePattern(const StiffnessMatrix &Ain, const int precond_num) { - if (!Ain.isCompressed()) - { - Ain.makeCompressed(); - } + assert(Ain.isCompressed()); + A = eigen2cholmod(Ain); L = cholmod_l_analyze(&A, cm); } - void CHOLMODSolver::factorize(StiffnessMatrixL &Ain) + void CHOLMODSolver::factorize(const StiffnessMatrix &Ain) { cholmod_l_factorize(&A, L, cm); } - void CHOLMODSolver::solve(Eigen::VectorXd &rhs, Eigen::VectorXd &result) + void CHOLMODSolver::solve(const Ref rhs, Ref result) { - b = eigen2cholmod(rhs); + assert(result.size() == rhs.size()); + + // b = eigen2cholmod(rhs); //TODO: fix me + x = cholmod_l_solve(CHOLMOD_A, L, &b, cm); - result.conservativeResize(rhs.size()); + memcpy(result.data(), x->x, result.size() * sizeof(result[0])); } diff --git a/src/polysolve/CHOLMODSolver.hpp b/src/polysolve/CHOLMODSolver.hpp index 5dc4ac70..b5134356 100644 --- a/src/polysolve/CHOLMODSolver.hpp +++ b/src/polysolve/CHOLMODSolver.hpp @@ -4,6 +4,8 @@ //////////////////////////////////////////////////////////////////////////////// +#include + #include #include @@ -15,11 +17,10 @@ using json = nlohmann::json; namespace polysolve { - typedef Eigen::SparseMatrix StiffnessMatrixL; /** @brief Base class for cholmod solver. - */ - class CHOLMODSolver + */ + class CHOLMODSolver : public LinearSolver { protected: @@ -32,17 +33,20 @@ namespace polysolve CHOLMODSolver(); ~CHOLMODSolver(); + private: + POLYSOLVE_DELETE_MOVE_COPY(CHOLMODSolver) + public: // Set solver parameters // virtual void setParameters(const json ¶ms) {} // Get info on the last solve step - void getInfo(json ¶ms) const; + virtual void getInfo(json ¶ms) const override; // Analyze sparsity pattern - void analyzePattern(StiffnessMatrixL &Ain); + virtual void analyzePattern(const StiffnessMatrix &A, const int precond_num) override; // Factorize system matrix - void factorize(StiffnessMatrixL &Ain); + void factorize(const StiffnessMatrix &A) override; // // @brief { Solve the linear system Ax = b } @@ -53,7 +57,7 @@ namespace polysolve // initial guess, and must thus be properly allocated // and initialized. } // - void solve(Eigen::VectorXd &rhs, Eigen::VectorXd &result); + virtual void solve(const Ref b, Ref x) override; }; } // namespace polysolve diff --git a/src/polysolve/LinearSolver.cpp b/src/polysolve/LinearSolver.cpp index 5dce34c8..0dd4d361 100644 --- a/src/polysolve/LinearSolver.cpp +++ b/src/polysolve/LinearSolver.cpp @@ -7,6 +7,7 @@ #include #ifdef POLYSOLVE_WITH_CHOLMOD #include +#include #endif #ifdef POLYSOLVE_WITH_UMFPACK #include @@ -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(); #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", diff --git a/tests/test_solver.cpp b/tests/test_solver.cpp index d9e38f26..2ba5178f 100644 --- a/tests/test_solver.cpp +++ b/tests/test_solver.cpp @@ -585,10 +585,10 @@ TEST_CASE("CHOLMOD", "[solver]") Eigen::VectorXd x(A.rows()); x.setZero(); - CHOLMODSolver solver; - solver.analyzePattern(A); - solver.factorize(A); - solver.solve(b, x); + 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);