From 7ac61fdd2f713f00aad34dfdb6e7d93982e907ac Mon Sep 17 00:00:00 2001 From: Kakeru Date: Sun, 28 Dec 2025 22:31:23 +0900 Subject: [PATCH 01/14] Add preconditioner interface --- .../r_KLU_rf_FGMRES_reuse_factorization.cpp | 4 +- examples/gpuRefactor.cpp | 4 +- examples/kluFactor.cpp | 4 +- examples/kluRefactor.cpp | 4 +- examples/randGmres.cpp | 4 +- resolve/CMakeLists.txt | 4 ++ resolve/LinSolverIterative.hpp | 13 +++-- resolve/LinSolverIterativeFGMRES.cpp | 24 +++++---- resolve/LinSolverIterativeFGMRES.hpp | 3 +- resolve/LinSolverIterativeRandFGMRES.cpp | 22 ++++---- resolve/LinSolverIterativeRandFGMRES.hpp | 3 +- resolve/Preconditioner.cpp | 20 +++++++ resolve/Preconditioner.hpp | 40 ++++++++++++++ resolve/PreconditionerLU.cpp | 53 +++++++++++++++++++ resolve/PreconditionerLU.hpp | 41 ++++++++++++++ resolve/SystemSolver.cpp | 30 ++++++++--- resolve/SystemSolver.hpp | 6 +-- tests/functionality/testKlu.cpp | 7 ++- tests/functionality/testRandGmres.cpp | 4 +- tests/functionality/testRefactor.cpp | 7 ++- 20 files changed, 251 insertions(+), 46 deletions(-) create mode 100644 resolve/Preconditioner.cpp create mode 100644 resolve/Preconditioner.hpp create mode 100644 resolve/PreconditionerLU.cpp create mode 100644 resolve/PreconditionerLU.hpp diff --git a/examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp b/examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp index 35baf7dc8..3baafd634 100644 --- a/examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp +++ b/examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -185,7 +186,8 @@ int main(int argc, char* argv[]) << status << std::endl; vec_rhs->copyDataFrom(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); status = Rf->solve(vec_rhs, vec_x); - FGMRES->setupPreconditioner("LU", Rf); + ReSolve::PreconditionerLU precond_lu(Rf); + FGMRES->setPreconditioner(&precond_lu); } // if (i%2!=0) vec_x->setToZero(ReSolve::memory::DEVICE); real_type norm_x = vector_handler->dot(vec_x, vec_x, ReSolve::memory::DEVICE); diff --git a/examples/gpuRefactor.cpp b/examples/gpuRefactor.cpp index 90f5ee06e..dd7d26c2b 100644 --- a/examples/gpuRefactor.cpp +++ b/examples/gpuRefactor.cpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -306,7 +307,8 @@ int gpuRefactor(int argc, char* argv[]) { // Setup iterative refinement FGMRES.resetMatrix(A); - FGMRES.setupPreconditioner("LU", &Rf); + ReSolve::PreconditionerLU precond_lu(&Rf); + FGMRES.setPreconditioner(&precond_lu); // If refactorization produced finite solution do iterative refinement if (std::isfinite(helper.getNormRelativeResidual())) diff --git a/examples/kluFactor.cpp b/examples/kluFactor.cpp index 78ad54d50..090ec127f 100644 --- a/examples/kluFactor.cpp +++ b/examples/kluFactor.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -191,7 +192,8 @@ int main(int argc, char* argv[]) { // Setup iterative refinement FGMRES.setup(A); - FGMRES.setupPreconditioner("LU", &KLU); + ReSolve::PreconditionerLU precond_lu(&KLU); + FGMRES.setPreconditioner(&precond_lu); // If refactorization produced finite solution do iterative refinement if (std::isfinite(helper.getNormRelativeResidual())) diff --git a/examples/kluRefactor.cpp b/examples/kluRefactor.cpp index ddf188d4a..992fc70fe 100644 --- a/examples/kluRefactor.cpp +++ b/examples/kluRefactor.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -197,7 +198,8 @@ int main(int argc, char* argv[]) { // Setup iterative refinement FGMRES.setup(A); - FGMRES.setupPreconditioner("LU", KLU); + ReSolve::PreconditionerLU precond_lu(KLU); + FGMRES.setPreconditioner(&precond_lu); // If refactorization produced finite solution do iterative refinement if (std::isfinite(helper.getNormRelativeResidual())) diff --git a/examples/randGmres.cpp b/examples/randGmres.cpp index 629cc2649..18845a075 100644 --- a/examples/randGmres.cpp +++ b/examples/randGmres.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -177,7 +178,8 @@ int runGmresExample(int argc, char* argv[]) FGMRES.setup(A); FGMRES.resetMatrix(A); - FGMRES.setupPreconditioner("LU", &Precond); + ReSolve::PreconditionerLU precond_lu(&Precond); + FGMRES.setPreconditioner(&precond_lu); FGMRES.setFlexible(1); FGMRES.solve(vec_rhs, vec_x); diff --git a/resolve/CMakeLists.txt b/resolve/CMakeLists.txt index ad45cf412..07dbf7c28 100644 --- a/resolve/CMakeLists.txt +++ b/resolve/CMakeLists.txt @@ -20,6 +20,8 @@ set(ReSolve_SRC LinSolverIterativeRandFGMRES.cpp LinSolverDirectSerialILU0.cpp SystemSolver.cpp + Preconditioner.cpp + PreconditionerLU.cpp ) set(ReSolve_KLU_SRC LinSolverDirectKLU.cpp) @@ -49,6 +51,8 @@ set(ReSolve_HEADER_INSTALL SystemSolver.hpp GramSchmidt.hpp MemoryUtils.hpp + Preconditioner.hpp + PreconditionerLU.hpp ) set(ReSolve_KLU_HEADER_INSTALL LinSolverDirectKLU.hpp) diff --git a/resolve/LinSolverIterative.hpp b/resolve/LinSolverIterative.hpp index 3cd1ef2ed..84eac784c 100644 --- a/resolve/LinSolverIterative.hpp +++ b/resolve/LinSolverIterative.hpp @@ -15,6 +15,7 @@ namespace ReSolve { class GramSchmidt; class LinSolverDirect; + class Preconditioner; class LinSolverIterative : public LinSolver { @@ -22,8 +23,8 @@ namespace ReSolve LinSolverIterative(); virtual ~LinSolverIterative(); virtual int setup(matrix::Sparse* A); - virtual int resetMatrix(matrix::Sparse* A) = 0; - virtual int setupPreconditioner(std::string type, LinSolverDirect* LU_solver) = 0; + virtual int resetMatrix(matrix::Sparse* A) = 0; + virtual int setPreconditioner(Preconditioner* preconditioner) = 0; virtual int solve(vector_type* rhs, vector_type* init_guess) = 0; @@ -40,9 +41,11 @@ namespace ReSolve void setMaxit(index_type new_maxit); protected: - real_type initial_residual_norm_; - real_type final_residual_norm_; - index_type total_iters_; + Preconditioner* preconditioner_{nullptr}; + + real_type initial_residual_norm_; + real_type final_residual_norm_; + index_type total_iters_; // Parameters common for all iterative solvers real_type tol_{1e-14}; ///< Solver tolerance diff --git a/resolve/LinSolverIterativeFGMRES.cpp b/resolve/LinSolverIterativeFGMRES.cpp index 643f23a46..fdbeaa1f5 100644 --- a/resolve/LinSolverIterativeFGMRES.cpp +++ b/resolve/LinSolverIterativeFGMRES.cpp @@ -17,6 +17,7 @@ #include #include #include +#include namespace ReSolve { @@ -317,18 +318,21 @@ namespace ReSolve return 0; } - int LinSolverIterativeFGMRES::setupPreconditioner(std::string type, LinSolverDirect* LU_solver) + /** + * @brief Sets pointer to Preconditioer. + * + * @param[in] precontitioner - pointer to Preconditioner class instance. + * @return 0 if successful, error code otherwise. + */ + int LinSolverIterativeFGMRES::setPreconditioner(Preconditioner* preconditioner) { - if (type != "LU") - { - out::warning() << "Only LU-type solve can be used as a preconditioner at this time." << std::endl; - return 1; - } - else + if (preconditioner == nullptr) { - LU_solver_ = LU_solver; - return 0; + out::warning() << "preconditioner pointer is null" << "\n"; + return 1; } + preconditioner_ = preconditioner; + return 0; } int LinSolverIterativeFGMRES::resetMatrix(matrix::Sparse* new_matrix) @@ -598,7 +602,7 @@ namespace ReSolve void LinSolverIterativeFGMRES::precV(vector_type* rhs, vector_type* x) { - LU_solver_->solve(rhs, x); + preconditioner_->apply(rhs, x); } void LinSolverIterativeFGMRES::setMemorySpace() diff --git a/resolve/LinSolverIterativeFGMRES.hpp b/resolve/LinSolverIterativeFGMRES.hpp index 388875793..89fca1479 100644 --- a/resolve/LinSolverIterativeFGMRES.hpp +++ b/resolve/LinSolverIterativeFGMRES.hpp @@ -16,6 +16,7 @@ namespace ReSolve // Forward declarations class SketchingHandler; class GramSchmidt; + class Preconditioner; namespace matrix { @@ -55,7 +56,7 @@ namespace ReSolve int solve(vector_type* rhs, vector_type* x) override; int setup(matrix::Sparse* A) override; int resetMatrix(matrix::Sparse* new_A) override; - int setupPreconditioner(std::string name, LinSolverDirect* LU_solver) override; + int setPreconditioner(Preconditioner* preconditioner) override; int setOrthogonalization(GramSchmidt* gs) override; int setRestart(index_type restart); diff --git a/resolve/LinSolverIterativeRandFGMRES.cpp b/resolve/LinSolverIterativeRandFGMRES.cpp index db26b0077..cdfd73863 100644 --- a/resolve/LinSolverIterativeRandFGMRES.cpp +++ b/resolve/LinSolverIterativeRandFGMRES.cpp @@ -18,6 +18,7 @@ #include #include #include +#include namespace ReSolve { @@ -406,18 +407,21 @@ namespace ReSolve return 0; } - int LinSolverIterativeRandFGMRES::setupPreconditioner(std::string type, LinSolverDirect* LU_solver) + /** + * @brief Sets pointer to Preconditioer. + * + * @param[in] precontitioner - pointer to Preconditioner class instance. + * @return 0 if successful, error code otherwise. + */ + int LinSolverIterativeRandFGMRES::setPreconditioner(Preconditioner* preconditioner) { - if (type != "LU") - { - out::warning() << "Only cusolverRf tri solve can be used as a preconditioner at this time." << std::endl; - return 1; - } - else + if (preconditioner == nullptr) { - LU_solver_ = LU_solver; - return 0; + out::warning() << "preconditioner pointer is null" << "\n"; + return 1; } + preconditioner_ = preconditioner; + return 0; } index_type LinSolverIterativeRandFGMRES::getKrand() diff --git a/resolve/LinSolverIterativeRandFGMRES.hpp b/resolve/LinSolverIterativeRandFGMRES.hpp index c16f77116..33bb97f47 100644 --- a/resolve/LinSolverIterativeRandFGMRES.hpp +++ b/resolve/LinSolverIterativeRandFGMRES.hpp @@ -16,6 +16,7 @@ namespace ReSolve // Forward declarations class SketchingHandler; class GramSchmidt; + class Preconditioner; namespace matrix { @@ -67,7 +68,7 @@ namespace ReSolve int solve(vector_type* rhs, vector_type* x) override; int setup(matrix::Sparse* A) override; int resetMatrix(matrix::Sparse* new_A) override; - int setupPreconditioner(std::string name, LinSolverDirect* LU_solver) override; + int setPreconditioner(Preconditioner* preconditoner) override; int setOrthogonalization(GramSchmidt* gs) override; int setRestart(index_type restart); diff --git a/resolve/Preconditioner.cpp b/resolve/Preconditioner.cpp new file mode 100644 index 000000000..8e2c38ab7 --- /dev/null +++ b/resolve/Preconditioner.cpp @@ -0,0 +1,20 @@ +/** + * @file Preconditioner.cpp + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Implementation of preconditioner base class. + * + */ + +#include "Preconditioner.hpp" + +namespace ReSolve +{ + Preconditioner::Preconditioner() + { + } + + Preconditioner::~Preconditioner() + { + } + +} // namespace ReSolve \ No newline at end of file diff --git a/resolve/Preconditioner.hpp b/resolve/Preconditioner.hpp new file mode 100644 index 000000000..c50955f12 --- /dev/null +++ b/resolve/Preconditioner.hpp @@ -0,0 +1,40 @@ +/** + * @file Preconditioner.hpp + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Declaration of preconditioner base class. + * + */ +#pragma once + +namespace ReSolve +{ + namespace matrix + { + class Sparse; + } // namespace matrix + + namespace vector + { + class Vector; + } // namespace vector + + /** + * @class Preconditioner + * + * @brief Interface for preconditioner. + */ + class Preconditioner + { + public: + using vector_type = vector::Vector; + using matrix_type = matrix::Sparse; + + Preconditioner(); + virtual ~Preconditioner(); + + virtual int setup(matrix_type* A) = 0; + virtual int reset(matrix_type* A) = 0; + virtual int apply(vector_type* rhs, vector_type* x) = 0; + }; +} // namespace ReSolve + diff --git a/resolve/PreconditionerLU.cpp b/resolve/PreconditionerLU.cpp new file mode 100644 index 000000000..0ea021981 --- /dev/null +++ b/resolve/PreconditionerLU.cpp @@ -0,0 +1,53 @@ +/** + * @file PreconditionerLU.cpp + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Declaration of preconditioner ILU0 class. + * + */ + +#include +#include "PreconditionerLU.hpp" + +namespace ReSolve +{ + PreconditionerLU::PreconditionerLU(LinSolverDirect* solver) + { + solver_ = solver; + } + + PreconditionerLU::~PreconditionerLU() + { + } + + int PreconditionerLU::setup(matrix::Sparse* A) + { + if (A == nullptr) + { + return 1; + } + solver_->setup(A); + + return 0; + } + + int PreconditionerLU::reset(matrix::Sparse* A) + { + if (solver_ == nullptr || A == nullptr) + { + return 1; + } + // LinSolverDirect doesn't have reset, so call setup instead + return solver_->setup(A); + } + + int PreconditionerLU::apply(vector_type* rhs, vector_type* x) + { + if (solver_ == nullptr) + { + return 1; + } + solver_->solve(rhs, x); + + return 0; + } +} // namespace ReSolve \ No newline at end of file diff --git a/resolve/PreconditionerLU.hpp b/resolve/PreconditionerLU.hpp new file mode 100644 index 000000000..d1c9fabdb --- /dev/null +++ b/resolve/PreconditionerLU.hpp @@ -0,0 +1,41 @@ +/** + * @file PreconditionerILU0.cpp + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Declaration of preconditioner ILU0 class. + * + */ + +#include + +namespace ReSolve +{ + // Forward declaration of workspace + class LinSolverDirect; + + namespace matrix + { + class Sparse; + } // namespace matrix + + namespace vector + { + class Vector; + } // namespace vector + + class PreconditionerLU : public Preconditioner + { + public: + using vector_type = vector::Vector; + using matrix_type = matrix::Sparse; + + PreconditionerLU(LinSolverDirect* solver); + ~PreconditionerLU(); + + int setup(matrix_type* A) override; + int reset(matrix_type* A) override; + int apply(vector_type* rhs, vector_type* x) override; + + private: + LinSolverDirect* solver_{nullptr}; + }; +} // namespace ReSolve diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 884c3fbd2..9e29caa19 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -171,6 +172,7 @@ namespace ReSolve { delete preconditioner_; } + delete matrixHandler_; delete vectorHandler_; @@ -229,6 +231,11 @@ namespace ReSolve delete refactorizationSolver_; refactorizationSolver_ = nullptr; } + if (preconditionSolver_) + { + delete preconditionSolver_; + preconditionSolver_ = nullptr; + } if (preconditioner_) { delete preconditioner_; @@ -312,19 +319,21 @@ namespace ReSolve { if (memspace_ == "cpu") { - // preconditioner_ = new LinSolverDirectSerialILU0(workspaceCpu_); - preconditioner_ = new LinSolverDirectCpuILU0(workspaceCpu_); + preconditionSolver_ = new LinSolverDirectSerialILU0(workspaceCpu_); + preconditioner_ = new PreconditionerLU(preconditionSolver_); #ifdef RESOLVE_USE_CUDA } else if (memspace_ == "cuda") { - preconditioner_ = new LinSolverDirectCuSparseILU0(workspaceCuda_); + preconditionSolver_ = new LinSolverDirectCuSparseILU0(workspaceCuda_); + preconditioner_ = new PreconditionerLU(preconditionSolver_); #endif #ifdef RESOLVE_USE_HIP } else if (memspace_ == "hip") { - preconditioner_ = new LinSolverDirectRocSparseILU0(workspaceHip_); + preconditionSolver_ = new LinSolverDirectRocSparseILU0(workspaceHip_); + preconditioner_ = new PreconditionerLU(preconditionSolver_); #endif } else @@ -482,9 +491,16 @@ namespace ReSolve if (irMethod_ == "fgmres") { status += iterativeSolver_->setup(A_); - status += iterativeSolver_->setupPreconditioner("LU", refactorizationSolver_); - } + if (preconditioner_) + { + delete preconditioner_; + preconditioner_ = nullptr; + } + + preconditioner_ = new PreconditionerLU(refactorizationSolver_); + status += iterativeSolver_->setPreconditioner(preconditioner_); + } return status; } @@ -550,7 +566,7 @@ namespace ReSolve { is_solve_on_device_ = true; } - iterativeSolver_->setupPreconditioner("LU", preconditioner_); + status += iterativeSolver_->setPreconditioner(preconditioner_); } return status; diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index 0aa5173c1..d81ea7fbd 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -10,6 +10,7 @@ namespace ReSolve class LinAlgWorkspaceCpu; class MatrixHandler; class VectorHandler; + class Preconditioner; namespace vector { @@ -27,8 +28,7 @@ namespace ReSolve using vector_type = vector::Vector; using matrix_type = matrix::Sparse; - /// @brief Temporary until abstract preconditioner class is created - using precond_type = LinSolverDirect; + using precond_type = Preconditioner; SystemSolver(LinAlgWorkspaceCpu* workspaceCpu, std::string factor = "klu", @@ -90,9 +90,9 @@ namespace ReSolve private: LinSolverDirect* factorizationSolver_{nullptr}; LinSolverDirect* refactorizationSolver_{nullptr}; + LinSolverDirect* preconditionSolver_{nullptr}; LinSolverIterative* iterativeSolver_{nullptr}; GramSchmidt* gs_{nullptr}; - precond_type* preconditioner_{nullptr}; LinAlgWorkspaceCUDA* workspaceCuda_{nullptr}; diff --git a/tests/functionality/testKlu.cpp b/tests/functionality/testKlu.cpp index 9bb6933b6..9ac09c3f7 100644 --- a/tests/functionality/testKlu.cpp +++ b/tests/functionality/testKlu.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -150,7 +151,8 @@ int runTest(int argc, char* argv[], std::string& solver_name) status = FGMRES.setup(A); error_sum += status; - status = FGMRES.setupPreconditioner("LU", &KLU); + ReSolve::PreconditionerLU precond_lu(&KLU); + status = FGMRES.setPreconditioner(&precond_lu); error_sum += status; status = FGMRES.solve(&vec_rhs, &vec_x); error_sum += status; @@ -197,7 +199,8 @@ int runTest(int argc, char* argv[], std::string& solver_name) if (is_ir) { FGMRES.resetMatrix(A); - status = FGMRES.setupPreconditioner("LU", &KLU); + ReSolve::PreconditionerLU precond_lu(&KLU); + status = FGMRES.setPreconditioner(&precond_lu); error_sum += status; status = FGMRES.solve(&vec_rhs, &vec_x); error_sum += status; diff --git a/tests/functionality/testRandGmres.cpp b/tests/functionality/testRandGmres.cpp index 35ed2e9c4..fe43b6c48 100644 --- a/tests/functionality/testRandGmres.cpp +++ b/tests/functionality/testRandGmres.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -132,7 +133,8 @@ int runTest(int argc, char* argv[]) FGMRES.setRestart(200); FGMRES.setSketchingMethod(LinSolverIterativeRandFGMRES::cs); - status = FGMRES.setupPreconditioner("LU", &ILU); + ReSolve::PreconditionerLU precond_lu(&ILU); + status = FGMRES.setPreconditioner(&precond_lu); error_sum += status; FGMRES.setFlexible(true); diff --git a/tests/functionality/testRefactor.cpp b/tests/functionality/testRefactor.cpp index aff7c2de5..c3e87c287 100644 --- a/tests/functionality/testRefactor.cpp +++ b/tests/functionality/testRefactor.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -190,7 +191,8 @@ int runTest(int argc, char* argv[], std::string& solver_name) status = FGMRES.setup(A); error_sum += status; - status = FGMRES.setupPreconditioner("LU", &Rf); + ReSolve::PreconditionerLU precond_lu(&Rf); + status = FGMRES.setPreconditioner(&precond_lu); error_sum += status; status = FGMRES.solve(&vec_rhs, &vec_x); @@ -239,7 +241,8 @@ int runTest(int argc, char* argv[], std::string& solver_name) if (is_ir) { FGMRES.resetMatrix(A); - status = FGMRES.setupPreconditioner("LU", &Rf); + ReSolve::PreconditionerLU precond_lu(&Rf); + status = FGMRES.setPreconditioner(&precond_lu); error_sum += status; status = FGMRES.solve(&vec_rhs, &vec_x); From c10a6f26fda9c9333b8f8382d9de2edfa49349a7 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Sun, 28 Dec 2025 22:38:19 +0900 Subject: [PATCH 02/14] Fix memory leaks --- resolve/LinSolverIterativeFGMRES.hpp | 1 - resolve/LinSolverIterativeRandFGMRES.cpp | 2 +- resolve/LinSolverIterativeRandFGMRES.hpp | 1 - 3 files changed, 1 insertion(+), 3 deletions(-) diff --git a/resolve/LinSolverIterativeFGMRES.hpp b/resolve/LinSolverIterativeFGMRES.hpp index 89fca1479..50a09ccfc 100644 --- a/resolve/LinSolverIterativeFGMRES.hpp +++ b/resolve/LinSolverIterativeFGMRES.hpp @@ -105,7 +105,6 @@ namespace ReSolve real_type* h_rs_{nullptr}; GramSchmidt* GS_{nullptr}; - LinSolverDirect* LU_solver_{nullptr}; index_type n_{0}; bool is_solver_set_{false}; diff --git a/resolve/LinSolverIterativeRandFGMRES.cpp b/resolve/LinSolverIterativeRandFGMRES.cpp index cdfd73863..0a76ff415 100644 --- a/resolve/LinSolverIterativeRandFGMRES.cpp +++ b/resolve/LinSolverIterativeRandFGMRES.cpp @@ -792,7 +792,7 @@ namespace ReSolve void LinSolverIterativeRandFGMRES::precV(vector_type* rhs, vector_type* x) { - LU_solver_->solve(rhs, x); + preconditioner_->apply(rhs, x); } /** diff --git a/resolve/LinSolverIterativeRandFGMRES.hpp b/resolve/LinSolverIterativeRandFGMRES.hpp index 33bb97f47..4160c7ce6 100644 --- a/resolve/LinSolverIterativeRandFGMRES.hpp +++ b/resolve/LinSolverIterativeRandFGMRES.hpp @@ -125,7 +125,6 @@ namespace ReSolve vector_type* vec_aux_{nullptr}; GramSchmidt* GS_{nullptr}; - LinSolverDirect* LU_solver_{nullptr}; index_type n_{0}; real_type one_over_k_{1.0}; From 7c4a3b6b5c16602bb6ec35dc324cc3cc80e0ed11 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Mon, 29 Dec 2025 08:29:17 +0900 Subject: [PATCH 03/14] Apply formatting --- resolve/LinSolverIterative.hpp | 6 +++--- resolve/LinSolverIterativeFGMRES.cpp | 10 +++++----- resolve/LinSolverIterativeFGMRES.hpp | 6 +++--- resolve/LinSolverIterativeRandFGMRES.cpp | 10 +++++----- resolve/LinSolverIterativeRandFGMRES.hpp | 6 +++--- resolve/Preconditioner.cpp | 2 +- resolve/Preconditioner.hpp | 19 +++++++++---------- resolve/PreconditionerLU.cpp | 5 +++-- resolve/SystemSolver.cpp | 1 - resolve/SystemSolver.hpp | 2 +- 10 files changed, 33 insertions(+), 34 deletions(-) diff --git a/resolve/LinSolverIterative.hpp b/resolve/LinSolverIterative.hpp index 84eac784c..7e38e9a9b 100644 --- a/resolve/LinSolverIterative.hpp +++ b/resolve/LinSolverIterative.hpp @@ -43,9 +43,9 @@ namespace ReSolve protected: Preconditioner* preconditioner_{nullptr}; - real_type initial_residual_norm_; - real_type final_residual_norm_; - index_type total_iters_; + real_type initial_residual_norm_; + real_type final_residual_norm_; + index_type total_iters_; // Parameters common for all iterative solvers real_type tol_{1e-14}; ///< Solver tolerance diff --git a/resolve/LinSolverIterativeFGMRES.cpp b/resolve/LinSolverIterativeFGMRES.cpp index fdbeaa1f5..7f7959e07 100644 --- a/resolve/LinSolverIterativeFGMRES.cpp +++ b/resolve/LinSolverIterativeFGMRES.cpp @@ -12,12 +12,12 @@ #include #include +#include #include #include #include #include #include -#include namespace ReSolve { @@ -320,7 +320,7 @@ namespace ReSolve /** * @brief Sets pointer to Preconditioer. - * + * * @param[in] precontitioner - pointer to Preconditioner class instance. * @return 0 if successful, error code otherwise. */ @@ -329,10 +329,10 @@ namespace ReSolve if (preconditioner == nullptr) { out::warning() << "preconditioner pointer is null" << "\n"; - return 1; + return 1; } - preconditioner_ = preconditioner; - return 0; + preconditioner_ = preconditioner; + return 0; } int LinSolverIterativeFGMRES::resetMatrix(matrix::Sparse* new_matrix) diff --git a/resolve/LinSolverIterativeFGMRES.hpp b/resolve/LinSolverIterativeFGMRES.hpp index 50a09ccfc..dfb115341 100644 --- a/resolve/LinSolverIterativeFGMRES.hpp +++ b/resolve/LinSolverIterativeFGMRES.hpp @@ -104,9 +104,9 @@ namespace ReSolve real_type* h_s_{nullptr}; real_type* h_rs_{nullptr}; - GramSchmidt* GS_{nullptr}; - index_type n_{0}; - bool is_solver_set_{false}; + GramSchmidt* GS_{nullptr}; + index_type n_{0}; + bool is_solver_set_{false}; MemoryHandler mem_; ///< Device memory manager object }; diff --git a/resolve/LinSolverIterativeRandFGMRES.cpp b/resolve/LinSolverIterativeRandFGMRES.cpp index 0a76ff415..ea44e4713 100644 --- a/resolve/LinSolverIterativeRandFGMRES.cpp +++ b/resolve/LinSolverIterativeRandFGMRES.cpp @@ -13,12 +13,12 @@ #include #include +#include #include #include #include #include #include -#include namespace ReSolve { @@ -409,7 +409,7 @@ namespace ReSolve /** * @brief Sets pointer to Preconditioer. - * + * * @param[in] precontitioner - pointer to Preconditioner class instance. * @return 0 if successful, error code otherwise. */ @@ -418,10 +418,10 @@ namespace ReSolve if (preconditioner == nullptr) { out::warning() << "preconditioner pointer is null" << "\n"; - return 1; + return 1; } - preconditioner_ = preconditioner; - return 0; + preconditioner_ = preconditioner; + return 0; } index_type LinSolverIterativeRandFGMRES::getKrand() diff --git a/resolve/LinSolverIterativeRandFGMRES.hpp b/resolve/LinSolverIterativeRandFGMRES.hpp index 4160c7ce6..d64b5f24f 100644 --- a/resolve/LinSolverIterativeRandFGMRES.hpp +++ b/resolve/LinSolverIterativeRandFGMRES.hpp @@ -124,9 +124,9 @@ namespace ReSolve real_type* h_rs_{nullptr}; vector_type* vec_aux_{nullptr}; - GramSchmidt* GS_{nullptr}; - index_type n_{0}; - real_type one_over_k_{1.0}; + GramSchmidt* GS_{nullptr}; + index_type n_{0}; + real_type one_over_k_{1.0}; index_type k_rand_{0}; ///< size of sketch space. We need to know it so we can allocate S! MemoryHandler mem_; ///< Device memory manager object diff --git a/resolve/Preconditioner.cpp b/resolve/Preconditioner.cpp index 8e2c38ab7..a97b9c795 100644 --- a/resolve/Preconditioner.cpp +++ b/resolve/Preconditioner.cpp @@ -17,4 +17,4 @@ namespace ReSolve { } -} // namespace ReSolve \ No newline at end of file +} // namespace ReSolve diff --git a/resolve/Preconditioner.hpp b/resolve/Preconditioner.hpp index c50955f12..8fe0dff82 100644 --- a/resolve/Preconditioner.hpp +++ b/resolve/Preconditioner.hpp @@ -20,21 +20,20 @@ namespace ReSolve /** * @class Preconditioner - * + * * @brief Interface for preconditioner. */ class Preconditioner { - public: - using vector_type = vector::Vector; - using matrix_type = matrix::Sparse; + public: + using vector_type = vector::Vector; + using matrix_type = matrix::Sparse; - Preconditioner(); - virtual ~Preconditioner(); + Preconditioner(); + virtual ~Preconditioner(); - virtual int setup(matrix_type* A) = 0; - virtual int reset(matrix_type* A) = 0; - virtual int apply(vector_type* rhs, vector_type* x) = 0; + virtual int setup(matrix_type* A) = 0; + virtual int reset(matrix_type* A) = 0; + virtual int apply(vector_type* rhs, vector_type* x) = 0; }; } // namespace ReSolve - diff --git a/resolve/PreconditionerLU.cpp b/resolve/PreconditionerLU.cpp index 0ea021981..66bb954bc 100644 --- a/resolve/PreconditionerLU.cpp +++ b/resolve/PreconditionerLU.cpp @@ -5,9 +5,10 @@ * */ -#include #include "PreconditionerLU.hpp" +#include + namespace ReSolve { PreconditionerLU::PreconditionerLU(LinSolverDirect* solver) @@ -50,4 +51,4 @@ namespace ReSolve return 0; } -} // namespace ReSolve \ No newline at end of file +} // namespace ReSolve diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 9e29caa19..18680e822 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -172,7 +172,6 @@ namespace ReSolve { delete preconditioner_; } - delete matrixHandler_; delete vectorHandler_; diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index d81ea7fbd..b6f57ebd0 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -93,7 +93,7 @@ namespace ReSolve LinSolverDirect* preconditionSolver_{nullptr}; LinSolverIterative* iterativeSolver_{nullptr}; GramSchmidt* gs_{nullptr}; - precond_type* preconditioner_{nullptr}; + precond_type* preconditioner_{nullptr}; LinAlgWorkspaceCUDA* workspaceCuda_{nullptr}; LinAlgWorkspaceHIP* workspaceHip_{nullptr}; From 63d713eee4b7a9243d3a934275c9543d91298c8e Mon Sep 17 00:00:00 2001 From: Kakeru Date: Wed, 31 Dec 2025 08:59:31 +0900 Subject: [PATCH 04/14] Temporarily remove reset function (to be implemented next) --- resolve/PreconditionerLU.cpp | 35 +++++++++++++++++++-------- resolve/PreconditionerLU.hpp | 1 - resolve/SystemSolver.hpp | 4 +-- tests/functionality/testRandGmres.cpp | 2 +- 4 files changed, 27 insertions(+), 15 deletions(-) diff --git a/resolve/PreconditionerLU.cpp b/resolve/PreconditionerLU.cpp index 66bb954bc..7e3b3d4b0 100644 --- a/resolve/PreconditionerLU.cpp +++ b/resolve/PreconditionerLU.cpp @@ -11,15 +11,30 @@ namespace ReSolve { + /** + * @brief Constructor for PreconditionerLU. + * + * @param[in] solver - Pointer to the LinSolverDirect object. + */ PreconditionerLU::PreconditionerLU(LinSolverDirect* solver) { solver_ = solver; } + /** + * @brief Destructor for PreconditionerLU + */ PreconditionerLU::~PreconditionerLU() { } + /** + * @brief Sets up the preconditioner with the given matrix + * + * @param[in] A - System matrix to set up the preconditioner with + * + * @return int 0 if successful, 1 if it fails + */ int PreconditionerLU::setup(matrix::Sparse* A) { if (A == nullptr) @@ -31,16 +46,16 @@ namespace ReSolve return 0; } - int PreconditionerLU::reset(matrix::Sparse* A) - { - if (solver_ == nullptr || A == nullptr) - { - return 1; - } - // LinSolverDirect doesn't have reset, so call setup instead - return solver_->setup(A); - } - + /** + * @brief Applies the preconditioner to solve the system Mx = rhs + * + * Computes x = M^(-1) * rhs where M is the preconditioner matrix. + * + * @param[in] rhs - Right-hand-side vector + * @param[in] x - Solution vector + * + * @return int 0 if successful, 1 if fails + */ int PreconditionerLU::apply(vector_type* rhs, vector_type* x) { if (solver_ == nullptr) diff --git a/resolve/PreconditionerLU.hpp b/resolve/PreconditionerLU.hpp index d1c9fabdb..5beb9f709 100644 --- a/resolve/PreconditionerLU.hpp +++ b/resolve/PreconditionerLU.hpp @@ -32,7 +32,6 @@ namespace ReSolve ~PreconditionerLU(); int setup(matrix_type* A) override; - int reset(matrix_type* A) override; int apply(vector_type* rhs, vector_type* x) override; private: diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index b6f57ebd0..6705bafc5 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -28,8 +28,6 @@ namespace ReSolve using vector_type = vector::Vector; using matrix_type = matrix::Sparse; - using precond_type = Preconditioner; - SystemSolver(LinAlgWorkspaceCpu* workspaceCpu, std::string factor = "klu", std::string refactor = "klu", @@ -93,7 +91,7 @@ namespace ReSolve LinSolverDirect* preconditionSolver_{nullptr}; LinSolverIterative* iterativeSolver_{nullptr}; GramSchmidt* gs_{nullptr}; - precond_type* preconditioner_{nullptr}; + Preconditioner* preconditioner_{nullptr}; LinAlgWorkspaceCUDA* workspaceCuda_{nullptr}; LinAlgWorkspaceHIP* workspaceHip_{nullptr}; diff --git a/tests/functionality/testRandGmres.cpp b/tests/functionality/testRandGmres.cpp index fe43b6c48..ca9adb464 100644 --- a/tests/functionality/testRandGmres.cpp +++ b/tests/functionality/testRandGmres.cpp @@ -133,7 +133,7 @@ int runTest(int argc, char* argv[]) FGMRES.setRestart(200); FGMRES.setSketchingMethod(LinSolverIterativeRandFGMRES::cs); - ReSolve::PreconditionerLU precond_lu(&ILU); + PreconditionerLU precond_lu(&ILU); status = FGMRES.setPreconditioner(&precond_lu); error_sum += status; From 2a531c7a05aea76691064ec0b07d12ec32700409 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Wed, 31 Dec 2025 11:16:19 +0900 Subject: [PATCH 05/14] Implement setPreconditioner() in LinSolverIterative --- resolve/LinSolverIterative.cpp | 10 ++++++++++ resolve/LinSolverIterative.hpp | 4 ++-- resolve/LinSolverIterativeFGMRES.cpp | 17 ----------------- resolve/LinSolverIterativeFGMRES.hpp | 1 - resolve/LinSolverIterativeRandFGMRES.cpp | 17 ----------------- resolve/LinSolverIterativeRandFGMRES.hpp | 1 - resolve/Preconditioner.hpp | 1 - 7 files changed, 12 insertions(+), 39 deletions(-) diff --git a/resolve/LinSolverIterative.cpp b/resolve/LinSolverIterative.cpp index c24af8f76..efc281b47 100644 --- a/resolve/LinSolverIterative.cpp +++ b/resolve/LinSolverIterative.cpp @@ -57,6 +57,16 @@ namespace ReSolve return maxit_; } + int LinSolverIterative::setPreconditioner(Preconditioner* preconditioner) + { + if (preconditioner == nullptr) + { + return 1; + } + preconditioner_ = preconditioner; + return 0; + } + int LinSolverIterative::setOrthogonalization(GramSchmidt* /* gs */) { out::error() << "Solver does not implement setting orthogonalization.\n"; diff --git a/resolve/LinSolverIterative.hpp b/resolve/LinSolverIterative.hpp index 7e38e9a9b..f881df97d 100644 --- a/resolve/LinSolverIterative.hpp +++ b/resolve/LinSolverIterative.hpp @@ -23,8 +23,7 @@ namespace ReSolve LinSolverIterative(); virtual ~LinSolverIterative(); virtual int setup(matrix::Sparse* A); - virtual int resetMatrix(matrix::Sparse* A) = 0; - virtual int setPreconditioner(Preconditioner* preconditioner) = 0; + virtual int resetMatrix(matrix::Sparse* A) = 0; virtual int solve(vector_type* rhs, vector_type* init_guess) = 0; @@ -32,6 +31,7 @@ namespace ReSolve virtual real_type getInitResidualNorm() const; virtual index_type getNumIter() const; + virtual int setPreconditioner(Preconditioner* preconditioner); virtual int setOrthogonalization(GramSchmidt* gs); real_type getTol() const; diff --git a/resolve/LinSolverIterativeFGMRES.cpp b/resolve/LinSolverIterativeFGMRES.cpp index 7f7959e07..871aa7c84 100644 --- a/resolve/LinSolverIterativeFGMRES.cpp +++ b/resolve/LinSolverIterativeFGMRES.cpp @@ -318,23 +318,6 @@ namespace ReSolve return 0; } - /** - * @brief Sets pointer to Preconditioer. - * - * @param[in] precontitioner - pointer to Preconditioner class instance. - * @return 0 if successful, error code otherwise. - */ - int LinSolverIterativeFGMRES::setPreconditioner(Preconditioner* preconditioner) - { - if (preconditioner == nullptr) - { - out::warning() << "preconditioner pointer is null" << "\n"; - return 1; - } - preconditioner_ = preconditioner; - return 0; - } - int LinSolverIterativeFGMRES::resetMatrix(matrix::Sparse* new_matrix) { A_ = new_matrix; diff --git a/resolve/LinSolverIterativeFGMRES.hpp b/resolve/LinSolverIterativeFGMRES.hpp index dfb115341..de9ebc8dd 100644 --- a/resolve/LinSolverIterativeFGMRES.hpp +++ b/resolve/LinSolverIterativeFGMRES.hpp @@ -56,7 +56,6 @@ namespace ReSolve int solve(vector_type* rhs, vector_type* x) override; int setup(matrix::Sparse* A) override; int resetMatrix(matrix::Sparse* new_A) override; - int setPreconditioner(Preconditioner* preconditioner) override; int setOrthogonalization(GramSchmidt* gs) override; int setRestart(index_type restart); diff --git a/resolve/LinSolverIterativeRandFGMRES.cpp b/resolve/LinSolverIterativeRandFGMRES.cpp index ea44e4713..e570cdaf3 100644 --- a/resolve/LinSolverIterativeRandFGMRES.cpp +++ b/resolve/LinSolverIterativeRandFGMRES.cpp @@ -407,23 +407,6 @@ namespace ReSolve return 0; } - /** - * @brief Sets pointer to Preconditioer. - * - * @param[in] precontitioner - pointer to Preconditioner class instance. - * @return 0 if successful, error code otherwise. - */ - int LinSolverIterativeRandFGMRES::setPreconditioner(Preconditioner* preconditioner) - { - if (preconditioner == nullptr) - { - out::warning() << "preconditioner pointer is null" << "\n"; - return 1; - } - preconditioner_ = preconditioner; - return 0; - } - index_type LinSolverIterativeRandFGMRES::getKrand() { return k_rand_; diff --git a/resolve/LinSolverIterativeRandFGMRES.hpp b/resolve/LinSolverIterativeRandFGMRES.hpp index d64b5f24f..e2da24b8e 100644 --- a/resolve/LinSolverIterativeRandFGMRES.hpp +++ b/resolve/LinSolverIterativeRandFGMRES.hpp @@ -68,7 +68,6 @@ namespace ReSolve int solve(vector_type* rhs, vector_type* x) override; int setup(matrix::Sparse* A) override; int resetMatrix(matrix::Sparse* new_A) override; - int setPreconditioner(Preconditioner* preconditoner) override; int setOrthogonalization(GramSchmidt* gs) override; int setRestart(index_type restart); diff --git a/resolve/Preconditioner.hpp b/resolve/Preconditioner.hpp index 8fe0dff82..63c867dde 100644 --- a/resolve/Preconditioner.hpp +++ b/resolve/Preconditioner.hpp @@ -33,7 +33,6 @@ namespace ReSolve virtual ~Preconditioner(); virtual int setup(matrix_type* A) = 0; - virtual int reset(matrix_type* A) = 0; virtual int apply(vector_type* rhs, vector_type* x) = 0; }; } // namespace ReSolve From 0383c0dfe89787bd12333e16239ec330bc0b5cf6 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Fri, 30 Jan 2026 09:50:43 -0500 Subject: [PATCH 06/14] Add reset function --- resolve/LinSolverDirect.cpp | 8 +++++++ resolve/LinSolverDirect.hpp | 1 + resolve/LinSolverDirectCpuILU0.hpp | 2 +- resolve/LinSolverDirectCuSparseILU0.hpp | 2 +- resolve/Preconditioner.cpp | 5 +++++ resolve/Preconditioner.hpp | 1 + resolve/PreconditionerLU.cpp | 20 ++++++++++++++++- resolve/PreconditionerLU.hpp | 1 + resolve/SystemSolver.cpp | 30 ++++++++++++++++++++++++- resolve/SystemSolver.hpp | 3 ++- 10 files changed, 68 insertions(+), 5 deletions(-) diff --git a/resolve/LinSolverDirect.cpp b/resolve/LinSolverDirect.cpp index d82e1046c..78a2c5727 100644 --- a/resolve/LinSolverDirect.cpp +++ b/resolve/LinSolverDirect.cpp @@ -77,6 +77,14 @@ namespace ReSolve return 1; } + /** + * @brief Placeholder function for resetting a matrix. + */ + int LinSolverDirect::reset(matrix::Sparse* A) + { + return 1; + } + /** * @brief Placeholder function for refactorization. */ diff --git a/resolve/LinSolverDirect.hpp b/resolve/LinSolverDirect.hpp index aafc03f8c..a7316d9e5 100644 --- a/resolve/LinSolverDirect.hpp +++ b/resolve/LinSolverDirect.hpp @@ -29,6 +29,7 @@ namespace ReSolve virtual int analyze(); // the same as symbolic factorization virtual int factorize(); virtual int refactorize(); + virtual int reset(matrix::Sparse* A); virtual int solve(vector_type* rhs, vector_type* x) = 0; virtual int solve(vector_type* x) = 0; diff --git a/resolve/LinSolverDirectCpuILU0.hpp b/resolve/LinSolverDirectCpuILU0.hpp index 5c8561bf1..d53367aad 100644 --- a/resolve/LinSolverDirectCpuILU0.hpp +++ b/resolve/LinSolverDirectCpuILU0.hpp @@ -56,7 +56,7 @@ namespace ReSolve index_type* Q = nullptr, vector_type* rhs = nullptr) override; // if values of A change, but the nnz pattern does not, redo the analysis only (reuse buffers though) - int reset(matrix::Sparse* A); + int reset(matrix::Sparse* A) override; int analyze() override; int factorize() override; diff --git a/resolve/LinSolverDirectCuSparseILU0.hpp b/resolve/LinSolverDirectCuSparseILU0.hpp index f4ae5ae6b..a56b23214 100644 --- a/resolve/LinSolverDirectCuSparseILU0.hpp +++ b/resolve/LinSolverDirectCuSparseILU0.hpp @@ -37,7 +37,7 @@ namespace ReSolve index_type* Q = nullptr, vector_type* rhs = nullptr) override; // if values of A change, but the nnz pattern does not, redo the analysis only (reuse buffers though) - int reset(matrix::Sparse* A); + int reset(matrix::Sparse* A) override; int solve(vector_type* rhs, vector_type* x) override; int solve(vector_type* rhs) override; diff --git a/resolve/Preconditioner.cpp b/resolve/Preconditioner.cpp index a97b9c795..2c6995606 100644 --- a/resolve/Preconditioner.cpp +++ b/resolve/Preconditioner.cpp @@ -17,4 +17,9 @@ namespace ReSolve { } + int Preconditioner::reset(matrix_type* A) + { + return 1; + } + } // namespace ReSolve diff --git a/resolve/Preconditioner.hpp b/resolve/Preconditioner.hpp index 63c867dde..9ea7a551b 100644 --- a/resolve/Preconditioner.hpp +++ b/resolve/Preconditioner.hpp @@ -34,5 +34,6 @@ namespace ReSolve virtual int setup(matrix_type* A) = 0; virtual int apply(vector_type* rhs, vector_type* x) = 0; + virtual int reset(matrix_type* A); }; } // namespace ReSolve diff --git a/resolve/PreconditionerLU.cpp b/resolve/PreconditionerLU.cpp index 7e3b3d4b0..fd6ec3c7c 100644 --- a/resolve/PreconditionerLU.cpp +++ b/resolve/PreconditionerLU.cpp @@ -35,7 +35,7 @@ namespace ReSolve * * @return int 0 if successful, 1 if it fails */ - int PreconditionerLU::setup(matrix::Sparse* A) + int PreconditionerLU::setup(matrix_type* A) { if (A == nullptr) { @@ -66,4 +66,22 @@ namespace ReSolve return 0; } + + /** + * @brief Resets the preconditioner with the given matrix + * + * @param[in] A - System matrix to reset the preconditioner with + * + * @return int 0 if successful, 1 if it fails + */ + int PreconditionerLU::reset(matrix_type* A) + { + if (A == nullptr) + { + return 1; + } + solver_->reset(A); + + return 0; + } } // namespace ReSolve diff --git a/resolve/PreconditionerLU.hpp b/resolve/PreconditionerLU.hpp index 5beb9f709..14ede135a 100644 --- a/resolve/PreconditionerLU.hpp +++ b/resolve/PreconditionerLU.hpp @@ -33,6 +33,7 @@ namespace ReSolve int setup(matrix_type* A) override; int apply(vector_type* rhs, vector_type* x) override; + int reset(matrix_type* A) override; private: LinSolverDirect* solver_{nullptr}; diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 18680e822..0153a44a6 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -177,7 +177,7 @@ namespace ReSolve delete vectorHandler_; } - int SystemSolver::setMatrix(matrix::Sparse* A) + int SystemSolver::setMatrix(matrix_type* A) { int status = 0; A_ = A; @@ -555,6 +555,13 @@ namespace ReSolve return status; } + /** + * @brief Sets up the preconditioner for the system solver + * + * Initializes and attaches the preconditioner to the iterative solver. + * + * @return int 0 if successful, 1 if it fails + */ int SystemSolver::preconditionerSetup() { int status = 0; @@ -571,6 +578,27 @@ namespace ReSolve return status; } + /** + * @brief Reset the preconditioner with a new matrix. + * + * Assumes the matrix sparsity pattern does not change. + * + * @param[in] A New sparse matrix (values updated). + * + * @return int 0 if successful, 1 if it fails + */ + int SystemSolver::resetPreconditioner(matrix_type* A) + { + int status = 0; + A_ = A; + if (precondition_method_ == "ilu0") + { + status += preconditioner_->reset(A); + } + + return status; + } + int SystemSolver::refine(vector_type* rhs, vector_type* x) { int status = 0; diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index 6705bafc5..ca0c313d8 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -50,12 +50,13 @@ namespace ReSolve ~SystemSolver(); int initialize(); - int setMatrix(matrix::Sparse* A); + int setMatrix(matrix_type* A); int analyze(); // symbolic part int factorize(); // numeric part int refactorize(); int refactorizationSetup(); int preconditionerSetup(); + int resetPreconditioner(matrix_type* A); int solve(vector_type* rhs, vector_type* x); // for direct and iterative int refine(vector_type* rhs, vector_type* x); // for iterative refinement From cd21cf5efffabc1605c994eb7a48acfd327e9cf6 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Fri, 30 Jan 2026 09:56:07 -0500 Subject: [PATCH 07/14] Apply formatting --- resolve/LinSolverDirect.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/LinSolverDirect.cpp b/resolve/LinSolverDirect.cpp index 78a2c5727..7b8d48f81 100644 --- a/resolve/LinSolverDirect.cpp +++ b/resolve/LinSolverDirect.cpp @@ -79,7 +79,7 @@ namespace ReSolve /** * @brief Placeholder function for resetting a matrix. - */ + */ int LinSolverDirect::reset(matrix::Sparse* A) { return 1; From 87f4dbd93d6328cf89e9c800214589987b69af8c Mon Sep 17 00:00:00 2001 From: Kakeru Date: Fri, 30 Jan 2026 10:17:12 -0500 Subject: [PATCH 08/14] Default matrix reset in base class --- resolve/LinSolverDirect.cpp | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/resolve/LinSolverDirect.cpp b/resolve/LinSolverDirect.cpp index 7b8d48f81..21d5183f7 100644 --- a/resolve/LinSolverDirect.cpp +++ b/resolve/LinSolverDirect.cpp @@ -62,25 +62,34 @@ namespace ReSolve } /** - * @brief Placeholder function for symbolic factorization. + * @brief Resets the matrix for the solver. + * + * @param[in] A - matrix to be reset + * + * @return int - error code, 0 if successful */ - int LinSolverDirect::analyze() + int LinSolverDirect::reset(matrix::Sparse* A) { - return 1; + if (A == nullptr) + { + return 1; + } + A_ = A; + return 0; } /** - * @brief Placeholder function for numeric factorization. + * @brief Placeholder function for symbolic factorization. */ - int LinSolverDirect::factorize() + int LinSolverDirect::analyze() { return 1; } /** - * @brief Placeholder function for resetting a matrix. + * @brief Placeholder function for numeric factorization. */ - int LinSolverDirect::reset(matrix::Sparse* A) + int LinSolverDirect::factorize() { return 1; } From 7a8755d48987399f4967d5b78e3b95c187514d9d Mon Sep 17 00:00:00 2001 From: Kakeru Date: Fri, 30 Jan 2026 10:41:40 -0500 Subject: [PATCH 09/14] Minor fix --- resolve/Preconditioner.cpp | 2 +- resolve/Preconditioner.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/resolve/Preconditioner.cpp b/resolve/Preconditioner.cpp index 2c6995606..4d85556d6 100644 --- a/resolve/Preconditioner.cpp +++ b/resolve/Preconditioner.cpp @@ -17,7 +17,7 @@ namespace ReSolve { } - int Preconditioner::reset(matrix_type* A) + int Preconditioner::reset(matrix_type* /* A */) { return 1; } diff --git a/resolve/Preconditioner.hpp b/resolve/Preconditioner.hpp index 9ea7a551b..4ed21cdad 100644 --- a/resolve/Preconditioner.hpp +++ b/resolve/Preconditioner.hpp @@ -34,6 +34,6 @@ namespace ReSolve virtual int setup(matrix_type* A) = 0; virtual int apply(vector_type* rhs, vector_type* x) = 0; - virtual int reset(matrix_type* A); + virtual int reset(matrix_type* /* A */); }; } // namespace ReSolve From 72e4a9f0a54b4230798c7fd198575b517f25acb3 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Fri, 30 Jan 2026 10:48:15 -0500 Subject: [PATCH 10/14] Update CHANGELOG.md --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f688ce34d..98cd0f408 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -67,3 +67,5 @@ It is seamless from the user perspective and fixed many bugs. 15. Added examples/sysGmres.cpp to demonstrate how to use SystemSolver with GMRES. 16. Updated MatrixHandler::addConst to return integer error codes instead of void. + +17. Added a preconditioner interface class so users can define thier own preconditioners. From 1c541e929efc8dd231623572c5b02d9bcd5e49cf Mon Sep 17 00:00:00 2001 From: Kakeru Date: Fri, 30 Jan 2026 14:24:10 -0500 Subject: [PATCH 11/14] Override LinSolverDirectRocSparseILU0::reset() --- resolve/LinSolverDirectRocSparseILU0.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/LinSolverDirectRocSparseILU0.hpp b/resolve/LinSolverDirectRocSparseILU0.hpp index 1a18a14d4..ccb376cd7 100644 --- a/resolve/LinSolverDirectRocSparseILU0.hpp +++ b/resolve/LinSolverDirectRocSparseILU0.hpp @@ -39,7 +39,7 @@ namespace ReSolve index_type* Q = nullptr, vector_type* rhs = nullptr) override; // if values of A change, but the nnz pattern does not, redo the analysis only (reuse buffers though) - int reset(matrix::Sparse* A); + int reset(matrix::Sparse* A) override; int solve(vector_type* rhs, vector_type* x) override; int solve(vector_type* rhs) override; // the solution is returned IN RHS (rhs is overwritten) From 4499205a8134bfbf6d04720b4fba2baaea91e495 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Fri, 30 Jan 2026 15:05:44 -0500 Subject: [PATCH 12/14] Add missing override --- resolve/LinSolverDirectCuSolverRf.hpp | 2 +- resolve/LinSolverDirectSerialILU0.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/resolve/LinSolverDirectCuSolverRf.hpp b/resolve/LinSolverDirectCuSolverRf.hpp index da51c1d88..bf9951ca4 100644 --- a/resolve/LinSolverDirectCuSolverRf.hpp +++ b/resolve/LinSolverDirectCuSolverRf.hpp @@ -37,7 +37,7 @@ namespace ReSolve matrix::Sparse* U, index_type* P, index_type* Q, - vector_type* rhs = nullptr); + vector_type* rhs = nullptr) override; int refactorize() override; int solve(vector_type* rhs, vector_type* x) override; diff --git a/resolve/LinSolverDirectSerialILU0.hpp b/resolve/LinSolverDirectSerialILU0.hpp index 64b96632a..f69aa2d93 100644 --- a/resolve/LinSolverDirectSerialILU0.hpp +++ b/resolve/LinSolverDirectSerialILU0.hpp @@ -34,7 +34,7 @@ namespace ReSolve index_type* P = nullptr, index_type* Q = nullptr, vector_type* rhs = nullptr) override; - int reset(matrix::Sparse* A); + int reset(matrix::Sparse* A) override; int solve(vector_type* rhs, vector_type* x) override; int solve(vector_type* rhs) override; // the solutuon is returned IN RHS (rhs is overwritten) From ce8308f7b2d5728cdec7fa2a9c0d926c47f0d2a9 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Fri, 30 Jan 2026 21:03:24 -0500 Subject: [PATCH 13/14] Fix memory leak --- resolve/PreconditionerLU.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/resolve/PreconditionerLU.cpp b/resolve/PreconditionerLU.cpp index fd6ec3c7c..d6fbfbe8a 100644 --- a/resolve/PreconditionerLU.cpp +++ b/resolve/PreconditionerLU.cpp @@ -26,6 +26,7 @@ namespace ReSolve */ PreconditionerLU::~PreconditionerLU() { + delete solver_; } /** From e9e46a2c6999bcb7d39a2f12bb5621b5c438882f Mon Sep 17 00:00:00 2001 From: Kakeru Date: Sat, 31 Jan 2026 08:33:55 -0500 Subject: [PATCH 14/14] Fix memory leak --- resolve/PreconditionerLU.cpp | 1 - resolve/SystemSolver.cpp | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/resolve/PreconditionerLU.cpp b/resolve/PreconditionerLU.cpp index d6fbfbe8a..fd6ec3c7c 100644 --- a/resolve/PreconditionerLU.cpp +++ b/resolve/PreconditionerLU.cpp @@ -26,7 +26,6 @@ namespace ReSolve */ PreconditionerLU::~PreconditionerLU() { - delete solver_; } /** diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 0153a44a6..ecf927c2a 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -171,6 +171,7 @@ namespace ReSolve if (precondition_method_ != "none") { delete preconditioner_; + delete preconditionSolver_; } delete matrixHandler_;