diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e43672b2..46238e50 100644 --- a/src/CMakeLists.txt +++ b/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) diff --git a/src/polysolve/CHOLMODSolver.cpp b/src/polysolve/CHOLMODSolver.cpp new file mode 100644 index 00000000..b2396398 --- /dev/null +++ b/src/polysolve/CHOLMODSolver.cpp @@ -0,0 +1,118 @@ +#ifdef POLYSOLVE_WITH_CHOLMOD + +#include "CHOLMODSolver.hpp" +#include + +#include +#include +#include + +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 ¶ms) 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 rhs, Ref 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 diff --git a/src/polysolve/CHOLMODSolver.hpp b/src/polysolve/CHOLMODSolver.hpp new file mode 100644 index 00000000..b5134356 --- /dev/null +++ b/src/polysolve/CHOLMODSolver.hpp @@ -0,0 +1,65 @@ +#pragma once + +#ifdef POLYSOLVE_WITH_CHOLMOD + +//////////////////////////////////////////////////////////////////////////////// + +#include + +#include +#include + +#include +#include + +#include +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 ¶ms) {} + + // Get info on the last solve step + virtual void getInfo(json ¶ms) 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 b, Ref x) override; + }; + +} // namespace polysolve + +#endif 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 df5478b3..2ba5178f 100644 --- a/tests/test_solver.cpp +++ b/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 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