Skip to content

Add preconditioner interface #405

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pelesh merged 14 commits into from
Feb 1, 2026
Merged

Add preconditioner interface #405

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
7ac61fd
Add preconditioner interface
KakeruUeda Dec 28, 2025
c10a6f2
Fix memory leaks
KakeruUeda Dec 28, 2025
7c4a3b6
Apply formatting
KakeruUeda Dec 28, 2025
63d713e
Temporarily remove reset function (to be implemented next)
KakeruUeda Dec 30, 2025
2a531c7
Implement setPreconditioner() in LinSolverIterative
KakeruUeda Dec 31, 2025
0383c0d
Add reset function
KakeruUeda Jan 30, 2026
cd21cf5
Apply formatting
KakeruUeda Jan 30, 2026
87f4dbd
Default matrix reset in base class
KakeruUeda Jan 30, 2026
7a8755d
Minor fix
KakeruUeda Jan 30, 2026
72e4a9f
Update CHANGELOG.md
KakeruUeda Jan 30, 2026
1c541e9
Override LinSolverDirectRocSparseILU0::reset()
KakeruUeda Jan 30, 2026
4499205
Add missing override
KakeruUeda Jan 30, 2026
ce8308f
Fix memory leak
KakeruUeda Jan 31, 2026
e9e46a2
Fix memory leak
KakeruUeda Jan 31, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
4 changes: 3 additions & 1 deletion examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#include <resolve/LinSolverDirectCuSolverRf.hpp>
#include <resolve/LinSolverDirectKLU.hpp>
#include <resolve/LinSolverIterativeFGMRES.hpp>
#include <resolve/PreconditionerLU.hpp>
#include <resolve/matrix/Coo.hpp>
#include <resolve/matrix/Csr.hpp>
#include <resolve/matrix/MatrixHandler.hpp>
Expand Down Expand Up @@ -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);
Expand Down
4 changes: 3 additions & 1 deletion examples/gpuRefactor.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
#include <resolve/GramSchmidt.hpp>
#include <resolve/LinSolverDirectKLU.hpp>
#include <resolve/LinSolverIterativeFGMRES.hpp>
#include <resolve/PreconditionerLU.hpp>
#include <resolve/Profiling.hpp>
#include <resolve/matrix/Csr.hpp>
#include <resolve/matrix/MatrixHandler.hpp>
Expand Down Expand Up @@ -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()))
Expand Down
4 changes: 3 additions & 1 deletion examples/kluFactor.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include <resolve/GramSchmidt.hpp>
#include <resolve/LinSolverDirectKLU.hpp>
#include <resolve/LinSolverIterativeFGMRES.hpp>
#include <resolve/PreconditionerLU.hpp>
#include <resolve/matrix/Csr.hpp>
#include <resolve/matrix/MatrixHandler.hpp>
#include <resolve/matrix/io.hpp>
Expand Down Expand Up @@ -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()))
Expand Down
4 changes: 3 additions & 1 deletion examples/kluRefactor.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include <resolve/GramSchmidt.hpp>
#include <resolve/LinSolverDirectKLU.hpp>
#include <resolve/LinSolverIterativeFGMRES.hpp>
#include <resolve/PreconditionerLU.hpp>
#include <resolve/matrix/Csr.hpp>
#include <resolve/matrix/MatrixHandler.hpp>
#include <resolve/matrix/io.hpp>
Expand Down Expand Up @@ -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()))
Expand Down
4 changes: 3 additions & 1 deletion examples/randGmres.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#include <resolve/LinSolverDirectCpuILU0.hpp>
#include <resolve/LinSolverDirectSerialILU0.hpp>
#include <resolve/LinSolverIterativeRandFGMRES.hpp>
#include <resolve/PreconditionerLU.hpp>
#include <resolve/matrix/Csr.hpp>
#include <resolve/matrix/MatrixHandler.hpp>
#include <resolve/matrix/io.hpp>
Expand Down Expand Up @@ -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);
Expand Down
4 changes: 4 additions & 0 deletions resolve/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@ set(ReSolve_SRC
LinSolverIterativeRandFGMRES.cpp
LinSolverDirectSerialILU0.cpp
SystemSolver.cpp
Preconditioner.cpp
PreconditionerLU.cpp
)

set(ReSolve_KLU_SRC LinSolverDirectKLU.cpp)
Expand Down Expand Up @@ -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)
Expand Down
17 changes: 17 additions & 0 deletions resolve/LinSolverDirect.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,23 @@ namespace ReSolve
return 0;
}

/**
* @brief Resets the matrix for the solver.
*
* @param[in] A - matrix to be reset
*
* @return int - error code, 0 if successful
*/
int LinSolverDirect::reset(matrix::Sparse* A)
{
if (A == nullptr)
{
return 1;
}
A_ = A;
return 0;
}

/**
* @brief Placeholder function for symbolic factorization.
*/
Expand Down
1 change: 1 addition & 0 deletions resolve/LinSolverDirect.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
2 changes: 1 addition & 1 deletion resolve/LinSolverDirectCpuILU0.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
2 changes: 1 addition & 1 deletion resolve/LinSolverDirectCuSolverRf.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
2 changes: 1 addition & 1 deletion resolve/LinSolverDirectCuSparseILU0.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
2 changes: 1 addition & 1 deletion resolve/LinSolverDirectRocSparseILU0.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion resolve/LinSolverDirectSerialILU0.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
10 changes: 10 additions & 0 deletions resolve/LinSolverIterative.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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";
Expand Down
7 changes: 5 additions & 2 deletions resolve/LinSolverIterative.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,22 +15,23 @@ namespace ReSolve
{
class GramSchmidt;
class LinSolverDirect;
class Preconditioner;

class LinSolverIterative : public LinSolver
{
public:
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 solve(vector_type* rhs, vector_type* init_guess) = 0;

virtual real_type getFinalResidualNorm() const;
virtual real_type getInitResidualNorm() const;
virtual index_type getNumIter() const;

virtual int setPreconditioner(Preconditioner* preconditioner);
virtual int setOrthogonalization(GramSchmidt* gs);

real_type getTol() const;
Expand All @@ -40,6 +41,8 @@ namespace ReSolve
void setMaxit(index_type new_maxit);

protected:
Preconditioner* preconditioner_{nullptr};

real_type initial_residual_norm_;
real_type final_residual_norm_;
index_type total_iters_;
Expand Down
17 changes: 2 additions & 15 deletions resolve/LinSolverIterativeFGMRES.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include <iostream>

#include <resolve/GramSchmidt.hpp>
#include <resolve/Preconditioner.hpp>
#include <resolve/matrix/MatrixHandler.hpp>
#include <resolve/matrix/Sparse.hpp>
#include <resolve/random/SketchingHandler.hpp>
Expand Down Expand Up @@ -317,20 +318,6 @@ namespace ReSolve
return 0;
}

int LinSolverIterativeFGMRES::setupPreconditioner(std::string type, LinSolverDirect* LU_solver)
{
if (type != "LU")
{
out::warning() << "Only LU-type solve can be used as a preconditioner at this time." << std::endl;
return 1;
}
else
{
LU_solver_ = LU_solver;
return 0;
}
}

int LinSolverIterativeFGMRES::resetMatrix(matrix::Sparse* new_matrix)
{
A_ = new_matrix;
Expand Down Expand Up @@ -598,7 +585,7 @@ namespace ReSolve

void LinSolverIterativeFGMRES::precV(vector_type* rhs, vector_type* x)
{
LU_solver_->solve(rhs, x);
preconditioner_->apply(rhs, x);
}

void LinSolverIterativeFGMRES::setMemorySpace()
Expand Down
9 changes: 4 additions & 5 deletions resolve/LinSolverIterativeFGMRES.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ namespace ReSolve
// Forward declarations
class SketchingHandler;
class GramSchmidt;
class Preconditioner;

namespace matrix
{
Expand Down Expand Up @@ -55,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 setupPreconditioner(std::string name, LinSolverDirect* LU_solver) override;
int setOrthogonalization(GramSchmidt* gs) override;

int setRestart(index_type restart);
Expand Down Expand Up @@ -103,10 +103,9 @@ namespace ReSolve
real_type* h_s_{nullptr};
real_type* h_rs_{nullptr};

GramSchmidt* GS_{nullptr};
LinSolverDirect* LU_solver_{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
};
Expand Down
17 changes: 2 additions & 15 deletions resolve/LinSolverIterativeRandFGMRES.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include <iostream>

#include <resolve/GramSchmidt.hpp>
#include <resolve/Preconditioner.hpp>
#include <resolve/matrix/MatrixHandler.hpp>
#include <resolve/matrix/Sparse.hpp>
#include <resolve/random/SketchingHandler.hpp>
Expand Down Expand Up @@ -406,20 +407,6 @@ namespace ReSolve
return 0;
}

int LinSolverIterativeRandFGMRES::setupPreconditioner(std::string type, LinSolverDirect* LU_solver)
{
if (type != "LU")
{
out::warning() << "Only cusolverRf tri solve can be used as a preconditioner at this time." << std::endl;
return 1;
}
else
{
LU_solver_ = LU_solver;
return 0;
}
}

index_type LinSolverIterativeRandFGMRES::getKrand()
{
return k_rand_;
Expand Down Expand Up @@ -788,7 +775,7 @@ namespace ReSolve

void LinSolverIterativeRandFGMRES::precV(vector_type* rhs, vector_type* x)
{
LU_solver_->solve(rhs, x);
preconditioner_->apply(rhs, x);
}

/**
Expand Down
9 changes: 4 additions & 5 deletions resolve/LinSolverIterativeRandFGMRES.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ namespace ReSolve
// Forward declarations
class SketchingHandler;
class GramSchmidt;
class Preconditioner;

namespace matrix
{
Expand Down Expand Up @@ -67,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 setupPreconditioner(std::string name, LinSolverDirect* LU_solver) override;
int setOrthogonalization(GramSchmidt* gs) override;

int setRestart(index_type restart);
Expand Down Expand Up @@ -123,10 +123,9 @@ namespace ReSolve
real_type* h_rs_{nullptr};
vector_type* vec_aux_{nullptr};

GramSchmidt* GS_{nullptr};
LinSolverDirect* LU_solver_{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
Expand Down
Loading
Loading