From 09d7a724f6b1c11cdcfd4192dfe3d89bbbd574b9 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Thu, 8 Oct 2026 13:45:51 +0900 Subject: [PATCH 1/5] Add identity preconditioner --- CHANGELOG.md | 1 + resolve/CMakeLists.txt | 2 + resolve/PreconditionerIdentity.cpp | 61 +++++++++++ resolve/PreconditionerIdentity.hpp | 35 ++++++ tests/unit/preconditioner/CMakeLists.txt | 10 ++ .../PreconditionerIdentityTests.hpp | 101 ++++++++++++++++++ .../runPreconditionerIdentityTests.cpp | 41 +++++++ 7 files changed, 251 insertions(+) create mode 100644 resolve/PreconditionerIdentity.cpp create mode 100644 resolve/PreconditionerIdentity.hpp create mode 100644 tests/unit/preconditioner/PreconditionerIdentityTests.hpp create mode 100644 tests/unit/preconditioner/runPreconditionerIdentityTests.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 72c0c4d7f..5f93f3b5a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Changes to Re::Solve since release 0.99.4 +- Added an identity preconditioner with CPU and GPU support. - Refactored system solver: Added helper functions to check and fix conflicts in configuration options and to reduce code duplication. - Code returns errors instead of warnings when encountering ambiguous erronneus inputs from the user. diff --git a/resolve/CMakeLists.txt b/resolve/CMakeLists.txt index 19eb51430..99761835d 100644 --- a/resolve/CMakeLists.txt +++ b/resolve/CMakeLists.txt @@ -20,6 +20,7 @@ set(ReSolve_SRC LinSolverIterativeRandFGMRES.cpp SystemSolver.cpp Preconditioner.cpp + PreconditionerIdentity.cpp PreconditionerLU.cpp PreconditionerUserMatrix.cpp ) @@ -55,6 +56,7 @@ set(ReSolve_HEADER_INSTALL GramSchmidt.hpp MemoryUtils.hpp Preconditioner.hpp + PreconditionerIdentity.hpp PreconditionerLU.hpp PreconditionerUserMatrix.hpp ) diff --git a/resolve/PreconditionerIdentity.cpp b/resolve/PreconditionerIdentity.cpp new file mode 100644 index 000000000..bc68bf90a --- /dev/null +++ b/resolve/PreconditionerIdentity.cpp @@ -0,0 +1,61 @@ +/** + * @file PreconditionerIdentity.cpp + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Implementation of identity preconditioner class. + * + */ + +#include "PreconditionerIdentity.hpp" + +#include + +#include +#include + +namespace ReSolve +{ + /** + * @brief Constructs an identity preconditioner for a memory space. + * + * @param[in] memspace Memory space used by input and output vectors + */ + PreconditionerIdentity::PreconditionerIdentity(memory::MemorySpace memspace) + : memspace_(memspace) + { + } + + /** + * @brief Sets up the identity preconditioner. + * + * The identity preconditioner has no matrix-dependent state, so setup is a + * no-op and accepts a null matrix pointer. + * + * @return Always returns zero. + */ + int PreconditionerIdentity::setup(matrix_type* /* A */) + { + return 0; + } + + /** + * @brief Applies the identity preconditioner. + * + * Copies rhs to x in the memory space specified at construction. + * + * @param[in] rhs Input vector + * @param[out] x Output vector + * + * @return Zero on success, one if the vectors are invalid or cannot be copied. + */ + int PreconditionerIdentity::apply(vector_type* rhs, vector_type* x) + { + assert(rhs != nullptr && "Input vector is null!"); + assert(x != nullptr && "Output vector is null!"); + assert(rhs->getSize() == x->getSize() && "Vector sizes do not match!"); + assert(rhs->getNumVectors() == x->getNumVectors() && "Vector counts do not match!"); + assert(rhs->isAllocated(memspace_) && "Input vector data is not allocated!"); + assert(x->isAllocated(memspace_) && "Output vector data is not allocated!"); + + return x->copyFromExternal(rhs, memspace_, memspace_) == 0 ? 0 : 1; + } +} // namespace ReSolve diff --git a/resolve/PreconditionerIdentity.hpp b/resolve/PreconditionerIdentity.hpp new file mode 100644 index 000000000..a223ab0fe --- /dev/null +++ b/resolve/PreconditionerIdentity.hpp @@ -0,0 +1,35 @@ +/** + * @file PreconditionerIdentity.hpp + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Declaration of identity preconditioner class. + * + */ + +#pragma once + +#include +#include + +namespace ReSolve +{ + /** + * @brief Identity preconditioner. + * + * Applying this preconditioner computes x = rhs. It can be used to run an + * iterative solver without changing its preconditioner-based code path. + */ + class PreconditionerIdentity : public Preconditioner + { + public: + using vector_type = vector::Vector; + using matrix_type = matrix::Sparse; + + explicit PreconditionerIdentity(memory::MemorySpace memspace); + + int setup(matrix_type* A) override; + int apply(vector_type* rhs, vector_type* x) override; + + private: + memory::MemorySpace memspace_; + }; +} // namespace ReSolve diff --git a/tests/unit/preconditioner/CMakeLists.txt b/tests/unit/preconditioner/CMakeLists.txt index 63df670ec..16062ff4f 100644 --- a/tests/unit/preconditioner/CMakeLists.txt +++ b/tests/unit/preconditioner/CMakeLists.txt @@ -7,14 +7,21 @@ add_executable( runPreconditionerUserMatrixTests.exe runPreconditionerUserMatrixTests.cpp ) +add_executable( + runPreconditionerIdentityTests.exe runPreconditionerIdentityTests.cpp +) add_executable(runPreconditionerLUTests.exe runPreconditionerLUTests.cpp) target_link_libraries(runPreconditionerUserMatrixTests.exe PRIVATE ReSolve) +target_link_libraries(runPreconditionerIdentityTests.exe PRIVATE ReSolve) target_link_libraries(runPreconditionerLUTests.exe PRIVATE ReSolve) install(TARGETS runPreconditionerUserMatrixTests.exe RUNTIME DESTINATION bin/resolve/tests/unit ) +install(TARGETS runPreconditionerIdentityTests.exe + RUNTIME DESTINATION bin/resolve/tests/unit +) install(TARGETS runPreconditionerLUTests.exe RUNTIME DESTINATION bin/resolve/tests/unit ) @@ -22,6 +29,9 @@ install(TARGETS runPreconditionerLUTests.exe add_test(NAME preconditioner_user_matrix_test COMMAND $ ) +add_test(NAME preconditioner_identity_test + COMMAND $ +) add_test(NAME preconditioner_lu_test COMMAND $ ) diff --git a/tests/unit/preconditioner/PreconditionerIdentityTests.hpp b/tests/unit/preconditioner/PreconditionerIdentityTests.hpp new file mode 100644 index 000000000..21ca19a7c --- /dev/null +++ b/tests/unit/preconditioner/PreconditionerIdentityTests.hpp @@ -0,0 +1,101 @@ +/** + * @file PreconditionerIdentityTests.hpp + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Tests for PreconditionerIdentity class. + * + */ + +#pragma once + +#include +#include + +#include +#include +#include + +namespace ReSolve +{ + namespace tests + { + /** + * @brief Unit tests for PreconditionerIdentity. + */ + class PreconditionerIdentityTests : public TestBase + { + public: + explicit PreconditionerIdentityTests(memory::MemorySpace memspace) + : memspace_(memspace) + { + } + + /** + * @brief Check that setup is independent of a system matrix. + */ + TestOutcome setup() + { + TestStatus status; + + PreconditionerIdentity preconditioner(memspace_); + status *= (preconditioner.setup(nullptr) == 0); + + return status.report(__func__); + } + + /** + * @brief Check the default and configurable preconditioning side. + */ + TestOutcome checkSide() + { + TestStatus status; + + PreconditionerIdentity preconditioner(memspace_); + status *= (preconditioner.getSide() == Preconditioner::Side::RIGHT); + status *= (preconditioner.setSide(Preconditioner::Side::LEFT) == 0); + status *= (preconditioner.getSide() == Preconditioner::Side::LEFT); + + return status.report(__func__); + } + + /** + * @brief Check that apply copies the input without changing its values. + */ + TestOutcome apply() + { + TestStatus status; + + const index_type n = 4; + const real_type values[n] = {1.0, -2.0, 3.5, 0.25}; + vector::Vector rhs(n); + vector::Vector x(n); + + rhs.allocate(memspace_); + x.allocate(memspace_); + rhs.copyFromExternal(values, memory::HOST, memspace_); + + PreconditionerIdentity preconditioner(memspace_); + status *= (preconditioner.apply(&rhs, &x) == 0); + + if (memspace_ == memory::DEVICE) + { + x.syncData(memory::HOST); + } + + for (index_type i = 0; i < n; ++i) + { + if (!isEqual(x.getData(memory::HOST)[i], values[i])) + { + std::cout << __func__ << ": x[" << i << "] = " + << x.getData(memory::HOST)[i] << ", expected " << values[i] << "\n"; + status *= false; + } + } + + return status.report(__func__); + } + + private: + memory::MemorySpace memspace_; + }; + } // namespace tests +} // namespace ReSolve diff --git a/tests/unit/preconditioner/runPreconditionerIdentityTests.cpp b/tests/unit/preconditioner/runPreconditionerIdentityTests.cpp new file mode 100644 index 000000000..195a0a72c --- /dev/null +++ b/tests/unit/preconditioner/runPreconditionerIdentityTests.cpp @@ -0,0 +1,41 @@ +/** + * @file runPreconditionerIdentityTests.cpp + * @brief Tests for PreconditionerIdentity class. + * + */ + +#include +#include + +#include "PreconditionerIdentityTests.hpp" + +/** + * @brief Run PreconditionerIdentity tests in a given memory space. + */ +void runTests(const std::string& backend, + ReSolve::memory::MemorySpace memspace, + ReSolve::tests::TestingResults& result) +{ + std::cout << "Running PreconditionerIdentity tests on " << backend << ":\n"; + + ReSolve::tests::PreconditionerIdentityTests test(memspace); + + result += test.setup(); + result += test.checkSide(); + result += test.apply(); + + std::cout << "\n"; +} + +int main(int, char**) +{ + ReSolve::tests::TestingResults result; + + runTests("CPU", ReSolve::memory::HOST, result); + +#if defined(RESOLVE_USE_CUDA) || defined(RESOLVE_USE_HIP) + runTests("GPU", ReSolve::memory::DEVICE, result); +#endif + + return result.summary(); +} From ee9c15531c15c9ea949a40062044fba3d98b7583 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Thu, 8 Oct 2026 20:42:53 +0900 Subject: [PATCH 2/5] Support unpreconditioned GMRES --- CHANGELOG.md | 2 +- examples/sysGmres.cpp | 41 ++++++++--- resolve/PreconditionerIdentity.cpp | 22 ++++-- resolve/PreconditionerIdentity.hpp | 1 + resolve/SystemSolver.cpp | 34 +++++++++- tests/functionality/CMakeLists.txt | 10 +++ tests/functionality/testSysGmres.cpp | 68 ++++++++++++++----- .../PreconditionerIdentityTests.hpp | 3 + 8 files changed, 145 insertions(+), 36 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5f93f3b5a..194a03548 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,9 +2,9 @@ ## Changes to Re::Solve since release 0.99.4 -- Added an identity preconditioner with CPU and GPU support. - Refactored system solver: Added helper functions to check and fix conflicts in configuration options and to reduce code duplication. - Code returns errors instead of warnings when encountering ambiguous erronneus inputs from the user. +- Added an identity preconditioner for unpreconditioned GMRES. ## HyKKT Release changes diff --git a/examples/sysGmres.cpp b/examples/sysGmres.cpp index 14d68bda4..d06e926bd 100644 --- a/examples/sysGmres.cpp +++ b/examples/sysGmres.cpp @@ -31,6 +31,7 @@ void printHelpInfo() std::cout << "sysGmres.exe -m -r \n\n"; std::cout << "Optional features:\n"; std::cout << "\t-b \tSelects hardware backend.\n"; + std::cout << "\t-c \tSelects preconditioner (default 'ilu0').\n"; std::cout << "\t-h \tPrints this message.\n"; std::cout << "\t-i \tIterative method: randgmres or fgmres (default 'randgmres').\n"; std::cout << "\t-g \tGram-Schmidt method: cgs1, cgs2, or mgs (default 'cgs2').\n"; @@ -63,7 +64,8 @@ static void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& flexible, - std::string& side); + std::string& side, + std::string& preconditioner); /// Processes backend-specific ILU0 zero-pivot options static int processILU0Inputs(SystemSolver& solver, @@ -179,10 +181,14 @@ int sysGmres(int argc, char* argv[]) opt = options.getParamFromKey("-p"); std::string side = opt ? (*opt).second : "right"; - processInputs(method, gs, sketch, flexible, side); + opt = options.getParamFromKey("-c"); + std::string preconditioner = opt ? (*opt).second : "ilu0"; + + processInputs(method, gs, sketch, flexible, side, preconditioner); std::cout << "Matrix file: " << matrix_pathname << "\n" - << "RHS file: " << rhs_pathname << "\n"; + << "RHS file: " << rhs_pathname << "\n" + << "Preconditioner: " << (preconditioner == "none" ? "none (identity)" : "ILU0") << "\n"; // Create workspace workspace_type workspace; @@ -201,15 +207,18 @@ int sysGmres(int argc, char* argv[]) "none", "none", method, - "ilu0", + preconditioner, "none"); solver.setGramSchmidtMethod(gs); - status = processILU0Inputs(solver, hw_backend, options); - if (status != 0) + if (preconditioner == "ilu0") { - return 1; + status = processILU0Inputs(solver, hw_backend, options); + if (status != 0) + { + return 1; + } } // Read and open matrix and right-hand-side vector @@ -268,8 +277,8 @@ int sysGmres(int argc, char* argv[]) solver.getIterativeSolver().setCliParam("flexible", flexible); solver.getIterativeSolver().setCliParam("restart", "200"); - // Set up the preconditioner - if (return_code == 0) + // Set up the requested non-identity preconditioner + if (return_code == 0 && preconditioner == "ilu0") { status = solver.preconditionerSetup(side); std::cout << "solver.preconditionerSetup returned status: " << status << "\n"; @@ -305,7 +314,12 @@ int sysGmres(int argc, char* argv[]) return return_code; } -void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& flexible, std::string& side) +void processInputs(std::string& method, + std::string& gs, + std::string& sketch, + std::string& flexible, + std::string& side, + std::string& preconditioner) { if (method == "randgmres") { @@ -344,6 +358,13 @@ void processInputs(std::string& method, std::string& gs, std::string& sketch, st std::cout << "Preconditioning side " << side << " not recognized.\n"; std::cout << "Setting preconditioning side to the default (right).\n\n"; } + + if (preconditioner != "ilu0" && preconditioner != "none") + { + std::cout << "Preconditioner " << preconditioner << " not recognized.\n"; + std::cout << "Setting preconditioner to the default (ILU0).\n\n"; + preconditioner = "ilu0"; + } } int processILU0Inputs(SystemSolver& solver, const std::string& hw_backend, const CliOptions& options) diff --git a/resolve/PreconditionerIdentity.cpp b/resolve/PreconditionerIdentity.cpp index bc68bf90a..3fb146e70 100644 --- a/resolve/PreconditionerIdentity.cpp +++ b/resolve/PreconditionerIdentity.cpp @@ -37,6 +37,19 @@ namespace ReSolve return 0; } + /** + * @brief Resets the identity preconditioner. + * + * The identity preconditioner has no matrix-dependent state, so reset is a + * no-op and accepts a null matrix pointer. + * + * @return Always returns zero. + */ + int PreconditionerIdentity::reset(matrix_type* /* A */) + { + return 0; + } + /** * @brief Applies the identity preconditioner. * @@ -49,12 +62,13 @@ namespace ReSolve */ int PreconditionerIdentity::apply(vector_type* rhs, vector_type* x) { - assert(rhs != nullptr && "Input vector is null!"); - assert(x != nullptr && "Output vector is null!"); + if (rhs == nullptr || x == nullptr) + { + return 1; + } + assert(rhs->getSize() == x->getSize() && "Vector sizes do not match!"); assert(rhs->getNumVectors() == x->getNumVectors() && "Vector counts do not match!"); - assert(rhs->isAllocated(memspace_) && "Input vector data is not allocated!"); - assert(x->isAllocated(memspace_) && "Output vector data is not allocated!"); return x->copyFromExternal(rhs, memspace_, memspace_) == 0 ? 0 : 1; } diff --git a/resolve/PreconditionerIdentity.hpp b/resolve/PreconditionerIdentity.hpp index a223ab0fe..9a6bebce3 100644 --- a/resolve/PreconditionerIdentity.hpp +++ b/resolve/PreconditionerIdentity.hpp @@ -27,6 +27,7 @@ namespace ReSolve explicit PreconditionerIdentity(memory::MemorySpace memspace); 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.cpp b/resolve/SystemSolver.cpp index b1ab0c826..a3f0ddf1e 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -236,6 +237,8 @@ namespace ReSolve // Make sure the method combination is consistent before creating objects validateConfiguration(); + const bool is_iterative_solve = (solve_method_ == "randgmres" || solve_method_ == "fgmres"); + // First delete old objects iterative_solver_.reset(); preconditioner_.reset(); @@ -269,7 +272,11 @@ namespace ReSolve // Create preconditioner if (precondition_method_ == "none") { - // do nothing + if (is_iterative_solve) + { + const memory::MemorySpace memspace = memspace_ == "cpu" ? memory::HOST : memory::DEVICE; + preconditioner_.reset(new PreconditionerIdentity(memspace)); + } } else if (precondition_method_ == "ilu0") { @@ -307,9 +314,17 @@ namespace ReSolve } // Create iterative solver - if (solve_method_ == "randgmres" || solve_method_ == "fgmres") + if (is_iterative_solve) { - createIterativeSolver(solve_method_); + if (createIterativeSolver(solve_method_) != 0) + { + return 1; + } + + if (precondition_method_ == "none") + { + return iterative_solver_->setPreconditioner(preconditioner_.get()); + } } return 0; @@ -723,12 +738,25 @@ namespace ReSolve ir_method_ = "none"; iterative_solver_.reset(); + if ((method == "fgmres" || method == "randgmres") && precondition_method_ == "none") + { + const memory::MemorySpace solve_memspace = + memspace_ == "cpu" ? memory::HOST : memory::DEVICE; + preconditioner_.reset(new PreconditionerIdentity(solve_memspace)); + } + if (createIterativeSolver(method) != 0) { out::error() << "Solve method " << solve_method_ << " not recognized ...\n"; return 1; } + + if (preconditioner_ != nullptr) + { + return iterative_solver_->setPreconditioner(preconditioner_.get()); + } + return 0; } diff --git a/tests/functionality/CMakeLists.txt b/tests/functionality/CMakeLists.txt index 38e664ca0..5ab7c8309 100644 --- a/tests/functionality/CMakeLists.txt +++ b/tests/functionality/CMakeLists.txt @@ -191,6 +191,16 @@ add_test(NAME sys_gmres_mgspm_test "-g" "mgs_pm" ) +# Unpreconditioned GMRES uses the identity preconditioner. +add_test(NAME sys_gmres_unpreconditioned_test + COMMAND $ "-x" "no" "-i" "fgmres" + "-g" "cgs2" "-c" "none" +) +add_test(NAME sys_rand_gmres_unpreconditioned_test + COMMAND $ "-x" "no" "-i" + "randgmres" "-g" "cgs2" "-s" "count" "-c" "none" +) + # Left-preconditioned Krylov solvers tests (GMRES) add_test(NAME sys_rand_count_gmres_left_cgs2_test COMMAND $ "-x" "no" "-i" diff --git a/tests/functionality/testSysGmres.cpp b/tests/functionality/testSysGmres.cpp index 4e2e04b5a..c494fe834 100644 --- a/tests/functionality/testSysGmres.cpp +++ b/tests/functionality/testSysGmres.cpp @@ -44,13 +44,18 @@ template static int test(int argc, char* argv[]); /// Checks if inputs are valid, otherwise sets defaults -static void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& side); +static void processInputs(std::string& method, + std::string& gs, + std::string& sketch, + std::string& side, + std::string& preconditioner); /// Creates string with test description static std::string headerInfo(const std::string& method, const std::string& gs, const std::string& sketch, - std::string flexible); + std::string flexible, + const std::string& preconditioner); /// Generates test system matrix static ReSolve::matrix::Csr* generateMatrix(const index_type N, MemorySpace memspace); @@ -108,7 +113,10 @@ int test(int argc, char* argv[]) opt = options.getParamFromKey("-p"); std::string side = opt ? (*opt).second : "right"; - processInputs(method, gs, sketch, side); + opt = options.getParamFromKey("-c"); + std::string preconditioner = opt ? (*opt).second : "ilu0"; + + processInputs(method, gs, sketch, side, preconditioner); // Create workspace and initialize its handles. workspace_type workspace; @@ -135,15 +143,15 @@ int test(int argc, char* argv[]) } // Create system solver - ReSolve::SystemSolver solver(&workspace, "none", "none", method, "ilu0", "none"); + ReSolve::SystemSolver solver(&workspace, "none", "none", method, preconditioner, "none"); solver.setGramSchmidtMethod(gs); // Configure ILU0 zero-pivot handling - if (hwbackend == "CPU") + if (preconditioner == "ilu0" && hwbackend == "CPU") { status = solver.getPreconditionerSolver().setCliParam("zero_diagonal", "1e-7"); } - else + else if (preconditioner == "ilu0") { status = solver.getPreconditionerSolver().setCliParam("numeric_boost", "yes"); status += solver.getPreconditionerSolver().setCliParam("boost_tolerance", "1e-7"); @@ -183,9 +191,12 @@ int test(int argc, char* argv[]) solver.getIterativeSolver().setCliParam("flexible", flexible); solver.getIterativeSolver().setCliParam("restart", "200"); - // Set preconditioner (default in this case ILU0) - status = solver.preconditionerSetup(side); - error_sum += status; + // Set up the requested non-identity preconditioner. + if (preconditioner == "ilu0") + { + status = solver.preconditionerSetup(side); + error_sum += status; + } // Solve system status = solver.solve(vec_rhs, &vec_x); @@ -203,7 +214,7 @@ int test(int argc, char* argv[]) bad_guess_x.allocateAll(memspace); bad_guess_x.setToConst(1.0e6, memspace); - ReSolve::SystemSolver bad_guess_solver(&workspace, "none", "none", method, "ilu0", "none"); + ReSolve::SystemSolver bad_guess_solver(&workspace, "none", "none", method, preconditioner, "none"); bad_guess_solver.setGramSchmidtMethod(gs); bad_guess_solver.getIterativeSolver().setCliParam("maxit", "2500"); bad_guess_solver.getIterativeSolver().setCliParam("tol", "1e-12"); @@ -221,8 +232,11 @@ int test(int argc, char* argv[]) bad_guess_solver.getIterativeSolver().setCliParam("flexible", flexible); bad_guess_solver.getIterativeSolver().setCliParam("restart", "200"); - status = bad_guess_solver.preconditionerSetup(side); - error_sum += status; + if (preconditioner == "ilu0") + { + status = bad_guess_solver.preconditionerSetup(side); + error_sum += status; + } const real_type bad_guess_rnorm = bad_guess_solver.getResidualNorm(vec_rhs, &bad_guess_x); @@ -261,7 +275,7 @@ int test(int argc, char* argv[]) vec_x_guess.copyFromExternal(&vec_x, memspace, memspace); vector_handler.scal(0.9, &vec_x_guess, memspace); - ReSolve::SystemSolver accepted_guess_solver(&workspace, "none", "none", method, "ilu0", "none"); + ReSolve::SystemSolver accepted_guess_solver(&workspace, "none", "none", method, preconditioner, "none"); accepted_guess_solver.setGramSchmidtMethod(gs); accepted_guess_solver.getIterativeSolver().setCliParam("maxit", "2500"); accepted_guess_solver.getIterativeSolver().setCliParam("tol", "1e-12"); @@ -279,8 +293,11 @@ int test(int argc, char* argv[]) accepted_guess_solver.getIterativeSolver().setCliParam("flexible", flexible); accepted_guess_solver.getIterativeSolver().setCliParam("restart", "200"); - status = accepted_guess_solver.preconditionerSetup(side); - error_sum += status; + if (preconditioner == "ilu0") + { + status = accepted_guess_solver.preconditionerSetup(side); + error_sum += status; + } const real_type initial_guess_rnorm = accepted_guess_solver.getResidualNorm(vec_rhs, &vec_x_guess); @@ -309,7 +326,7 @@ int test(int argc, char* argv[]) helper.setSystem(A, vec_rhs, &vec_x); std::cout << std::defaultfloat - << headerInfo(method, gs, sketch, flexible) + << headerInfo(method, gs, sketch, flexible, preconditioner) << "\t Hardware backend: " << hwbackend << "\n" << "\t Solver tolerance: " << tol_out << "\n"; helper.printIterativeSolverSummary(&(solver.getIterativeSolver())); @@ -328,7 +345,11 @@ int test(int argc, char* argv[]) // Definitions of helper functions // -void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& side) +void processInputs(std::string& method, + std::string& gs, + std::string& sketch, + std::string& side, + std::string& preconditioner) { if (method == "randgmres") { @@ -361,12 +382,20 @@ void processInputs(std::string& method, std::string& gs, std::string& sketch, st std::cout << "Setting preconditioning side to the default (right).\n\n"; side = "right"; } + + if (preconditioner != "ilu0" && preconditioner != "none") + { + std::cout << "Preconditioner " << preconditioner << " not recognized.\n"; + std::cout << "Setting preconditioner to the default (ILU0).\n\n"; + preconditioner = "ilu0"; + } } std::string headerInfo(const std::string& method, const std::string& gs, const std::string& sketch, - std::string flexible) + std::string flexible, + const std::string& preconditioner) { bool is_flexible = !(flexible == "no"); std::string header("Results for "); @@ -421,6 +450,9 @@ std::string headerInfo(const std::string& method, // do nothing } + header += "\t Preconditioner: "; + header += preconditioner == "none" ? "identity\n" : "ILU0\n"; + return header; } diff --git a/tests/unit/preconditioner/PreconditionerIdentityTests.hpp b/tests/unit/preconditioner/PreconditionerIdentityTests.hpp index 21ca19a7c..8cc2208c5 100644 --- a/tests/unit/preconditioner/PreconditionerIdentityTests.hpp +++ b/tests/unit/preconditioner/PreconditionerIdentityTests.hpp @@ -38,6 +38,7 @@ namespace ReSolve PreconditionerIdentity preconditioner(memspace_); status *= (preconditioner.setup(nullptr) == 0); + status *= (preconditioner.reset(nullptr) == 0); return status.report(__func__); } @@ -74,6 +75,8 @@ namespace ReSolve rhs.copyFromExternal(values, memory::HOST, memspace_); PreconditionerIdentity preconditioner(memspace_); + status *= (preconditioner.apply(nullptr, &x) == 1); + status *= (preconditioner.apply(&rhs, nullptr) == 1); status *= (preconditioner.apply(&rhs, &x) == 0); if (memspace_ == memory::DEVICE) From e32c09060f13e82c2a050c8e11c2e11a951afb09 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Thu, 8 Oct 2026 21:58:26 +0900 Subject: [PATCH 3/5] Update comments --- resolve/PreconditionerIdentity.cpp | 4 ++-- resolve/PreconditionerIdentity.hpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/resolve/PreconditionerIdentity.cpp b/resolve/PreconditionerIdentity.cpp index 3fb146e70..a141130d2 100644 --- a/resolve/PreconditionerIdentity.cpp +++ b/resolve/PreconditionerIdentity.cpp @@ -15,9 +15,9 @@ namespace ReSolve { /** - * @brief Constructs an identity preconditioner for a memory space. + * @brief Constructs an identity preconditioner. * - * @param[in] memspace Memory space used by input and output vectors + * @param[in] memspace Memory space used by input and output vectors. */ PreconditionerIdentity::PreconditionerIdentity(memory::MemorySpace memspace) : memspace_(memspace) diff --git a/resolve/PreconditionerIdentity.hpp b/resolve/PreconditionerIdentity.hpp index 9a6bebce3..1abf54f72 100644 --- a/resolve/PreconditionerIdentity.hpp +++ b/resolve/PreconditionerIdentity.hpp @@ -15,8 +15,8 @@ namespace ReSolve /** * @brief Identity preconditioner. * - * Applying this preconditioner computes x = rhs. It can be used to run an - * iterative solver without changing its preconditioner-based code path. + * This preconditioner uses the identity matrix, + * so applying it copies rhs to x. */ class PreconditionerIdentity : public Preconditioner { From c336c9d0ae9bd6054adb79a9d4f8b82229b6ca70 Mon Sep 17 00:00:00 2001 From: Kakeru Date: Thu, 8 Oct 2026 22:09:25 +0900 Subject: [PATCH 4/5] Refactor iterative solve method checks --- resolve/SystemSolver.cpp | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index a3f0ddf1e..2ef541cdc 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -133,6 +133,13 @@ namespace ReSolve return "hip"; } #endif + + /// Returns true when the method selects a supported iterative solver. + bool isIterativeSolve(const std::string& method) + { + return method == "fgmres" || method == "randgmres"; + } + } // namespace SystemSolver::SystemSolver(LinAlgWorkspaceCpu* workspace_cpu, @@ -211,7 +218,7 @@ namespace ReSolve } // If we use iterative solver, we can set it up here - if (solve_method_ == "randgmres" || solve_method_ == "fgmres") + if (isIterativeSolve(solve_method_)) { if (iterative_solver_ == nullptr) { @@ -237,8 +244,6 @@ namespace ReSolve // Make sure the method combination is consistent before creating objects validateConfiguration(); - const bool is_iterative_solve = (solve_method_ == "randgmres" || solve_method_ == "fgmres"); - // First delete old objects iterative_solver_.reset(); preconditioner_.reset(); @@ -272,7 +277,7 @@ namespace ReSolve // Create preconditioner if (precondition_method_ == "none") { - if (is_iterative_solve) + if (isIterativeSolve(solve_method_)) { const memory::MemorySpace memspace = memspace_ == "cpu" ? memory::HOST : memory::DEVICE; preconditioner_.reset(new PreconditionerIdentity(memspace)); @@ -314,7 +319,7 @@ namespace ReSolve } // Create iterative solver - if (is_iterative_solve) + if (isIterativeSolve(solve_method_)) { if (createIterativeSolver(solve_method_) != 0) { @@ -491,7 +496,7 @@ namespace ReSolve int status = 0; // Use Krylov solver if selected - if (solve_method_ == "randgmres" || solve_method_ == "fgmres") + if (isIterativeSolve(solve_method_)) { status += iterative_solver_->resetMatrix(A_); status += iterative_solver_->solve(rhs, x); @@ -666,7 +671,7 @@ namespace ReSolve */ int SystemSolver::setFactorizationMethod(std::string method) { - if (solve_method_ == "fgmres" || solve_method_ == "randgmres") + if (isIterativeSolve(solve_method_)) { out::error() << "Factorization method cannot be set while iterative solve method '" << solve_method_ << "' is active. Keeping '" << factorization_method_ << "'.\n"; @@ -738,7 +743,7 @@ namespace ReSolve ir_method_ = "none"; iterative_solver_.reset(); - if ((method == "fgmres" || method == "randgmres") && precondition_method_ == "none") + if (isIterativeSolve(method) && precondition_method_ == "none") { const memory::MemorySpace solve_memspace = memspace_ == "cpu" ? memory::HOST : memory::DEVICE; @@ -773,7 +778,7 @@ namespace ReSolve { // With an iterative solve method the Krylov solver belongs to the solve // path, not to iterative refinement, so leave it untouched. - if (solve_method_ == "randgmres" || solve_method_ == "fgmres") + if (isIterativeSolve(solve_method_)) { if (method != "none") { @@ -1092,9 +1097,7 @@ namespace ReSolve */ void SystemSolver::validateConfiguration() { - const bool is_iterative_solve = (solve_method_ == "fgmres" || solve_method_ == "randgmres"); - - if (is_iterative_solve) + if (isIterativeSolve(solve_method_)) { if (factorization_method_ != "none") { From 5fed171d7d0d2ce4781e6ce2d2041f380d59a97c Mon Sep 17 00:00:00 2001 From: Kakeru Date: Thu, 8 Oct 2026 22:28:35 +0900 Subject: [PATCH 5/5] Improve naming --- examples/randGmres.cpp | 4 ++-- resolve/SystemSolver.cpp | 22 +++++++++++----------- resolve/SystemSolver.hpp | 2 +- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/examples/randGmres.cpp b/examples/randGmres.cpp index b3092a920..1fd25d862 100644 --- a/examples/randGmres.cpp +++ b/examples/randGmres.cpp @@ -116,7 +116,7 @@ int runGmresExample(int argc, char* argv[]) GramSchmidt GS(&vector_handler, GramSchmidt::CGS2); - precon_type precondition_solver(&workspace); + precon_type preconditioner_solver(&workspace); LinSolverIterativeRandFGMRES FGMRES(&matrix_handler, &vector_handler, LinSolverIterativeRandFGMRES::cs, @@ -179,7 +179,7 @@ int runGmresExample(int argc, char* argv[]) FGMRES.setup(A); FGMRES.resetMatrix(A); - ReSolve::PreconditionerLU preconditioner(&precondition_solver); + ReSolve::PreconditionerLU preconditioner(&preconditioner_solver); preconditioner.setup(A); FGMRES.setPreconditioner(&preconditioner); diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 2ef541cdc..e0819b776 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -152,7 +152,7 @@ namespace ReSolve factorization_method_(factor), refactorization_method_(refactor), solve_method_(solve), - precondition_method_(precond), + preconditioner_method_(precond), ir_method_(ir) { completeSetup(workspace_cpu_); @@ -169,7 +169,7 @@ namespace ReSolve factorization_method_(factor), refactorization_method_(refactor), solve_method_(solve), - precondition_method_(precond), + preconditioner_method_(precond), ir_method_(ir) { completeSetup(workspace_cuda_); @@ -187,7 +187,7 @@ namespace ReSolve factorization_method_(factor), refactorization_method_(refactor), solve_method_(solve), - precondition_method_(precond), + preconditioner_method_(precond), ir_method_(ir) { completeSetup(workspace_hip_); @@ -275,7 +275,7 @@ namespace ReSolve } // Create preconditioner - if (precondition_method_ == "none") + if (preconditioner_method_ == "none") { if (isIterativeSolve(solve_method_)) { @@ -283,7 +283,7 @@ namespace ReSolve preconditioner_.reset(new PreconditionerIdentity(memspace)); } } - else if (precondition_method_ == "ilu0") + else if (preconditioner_method_ == "ilu0") { if (memspace_ == "cpu") { @@ -313,7 +313,7 @@ namespace ReSolve } else { - out::error() << "Preconditioner method " << precondition_method_ + out::error() << "Preconditioner method " << preconditioner_method_ << " not recognized ...\n"; return 1; } @@ -326,7 +326,7 @@ namespace ReSolve return 1; } - if (precondition_method_ == "none") + if (preconditioner_method_ == "none") { return iterative_solver_->setPreconditioner(preconditioner_.get()); } @@ -743,7 +743,7 @@ namespace ReSolve ir_method_ = "none"; iterative_solver_.reset(); - if (isIterativeSolve(method) && precondition_method_ == "none") + if (isIterativeSolve(method) && preconditioner_method_ == "none") { const memory::MemorySpace solve_memspace = memspace_ == "cpu" ? memory::HOST : memory::DEVICE; @@ -1076,7 +1076,7 @@ namespace ReSolve out::error() << "SystemSolver initialization failed with factorization '" << factorization_method_ << "', refactorization '" << refactorization_method_ << "', solve '" << solve_method_ - << "', preconditioner '" << precondition_method_ + << "', preconditioner '" << preconditioner_method_ << "', iterative refinement '" << ir_method_ << "'. Solver is not usable in this state.\n"; } @@ -1120,9 +1120,9 @@ namespace ReSolve } else { - if (precondition_method_ != "none") + if (preconditioner_method_ != "none") { - out::error() << "Incorrect input: preconditioner '" << precondition_method_ + out::error() << "Incorrect input: preconditioner '" << preconditioner_method_ << "' can only be used with an iterative solve method ('fgmres' or 'randgmres'). " << "\n"; } diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index b19eceacc..98ba6d9f1 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -133,7 +133,7 @@ namespace ReSolve std::string factorization_method_{"none"}; std::string refactorization_method_{"none"}; std::string solve_method_{"none"}; - std::string precondition_method_{"none"}; + std::string preconditioner_method_{"none"}; std::string ir_method_{"none"}; std::string gs_method_{"cgs2"}; std::string sketching_method_{"count"}; ///< @todo move this to LinSolverIterative class