From 3527b317a29ac6f05e03b766b4337b63863b8fec Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 16:50:04 -0400 Subject: [PATCH 01/16] Use unique pointers for objects owned by SystemSolver. --- resolve/SystemSolver.cpp | 237 +++++++++++++-------------------------- resolve/SystemSolver.hpp | 24 ++-- 2 files changed, 90 insertions(+), 171 deletions(-) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 9aa7f8dd8..1a571d0f5 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -75,8 +75,8 @@ namespace ReSolve } // Instantiate handlers - matrixHandler_ = new MatrixHandler(workspaceCpu_); - vectorHandler_ = new VectorHandler(workspaceCpu_); + matrixHandler_.reset(new MatrixHandler(workspaceCpu_)); + vectorHandler_.reset(new VectorHandler(workspaceCpu_)); memspace_ = "cpu"; @@ -115,8 +115,8 @@ namespace ReSolve } // Instantiate handlers - matrixHandler_ = new MatrixHandler(workspaceCuda_); - vectorHandler_ = new VectorHandler(workspaceCuda_); + matrixHandler_.reset(new MatrixHandler(workspaceCuda_)); + vectorHandler_.reset(new VectorHandler(workspaceCuda_)); memspace_ = "cuda"; @@ -156,8 +156,8 @@ namespace ReSolve } // Instantiate handlers - matrixHandler_ = new MatrixHandler(workspaceHip_); - vectorHandler_ = new VectorHandler(workspaceHip_); + matrixHandler_.reset(new MatrixHandler(workspaceHip_)); + vectorHandler_.reset(new VectorHandler(workspaceHip_)); memspace_ = "hip"; @@ -165,57 +165,19 @@ namespace ReSolve } #endif - SystemSolver::~SystemSolver() - { - if (resVector_ != nullptr) - { - delete resVector_; - } - - if (factorizationMethod_ != "none") - { - delete factorizationSolver_; - } - - if (refactorizationMethod_ != "none") - { - delete refactorizationSolver_; - } - - if (solveMethod_ == "randgmres" || solveMethod_ == "fgmres") - { - delete iterativeSolver_; - } - - if (gs_ != nullptr) - { - delete gs_; - } - - if (irMethod_ != "none") - { - delete iterativeSolver_; - } - - if (precondition_method_ != "none") - { - delete preconditioner_; - delete preconditionerSolver_; - } - - delete matrixHandler_; - delete vectorHandler_; - } + /** + * @brief Destructor + * + * All owned components are held in `std::unique_ptr` and are released + * automatically in reverse declaration order. + */ + SystemSolver::~SystemSolver() = default; int SystemSolver::setMatrix(matrix_type* A) { int status = 0; A_ = A; - if (resVector_) - { - delete resVector_; - } - resVector_ = new vector_type(A->getNumRows()); + resVector_.reset(new vector_type(A->getNumRows())); if (memspace_ == "cpu") { resVector_->allocate(memory::HOST); @@ -229,12 +191,12 @@ namespace ReSolve // If we use iterative solver, we can set it up here if (solveMethod_ == "randgmres") { - auto* rgmres = dynamic_cast(iterativeSolver_); + auto* rgmres = dynamic_cast(iterativeSolver_.get()); status += rgmres->setup(A_); } else if (solveMethod_ == "fgmres") { - auto* fgmres = dynamic_cast(iterativeSolver_); + auto* fgmres = dynamic_cast(iterativeSolver_.get()); status += fgmres->setup(A_); } else @@ -254,36 +216,12 @@ namespace ReSolve int SystemSolver::initialize() { // First delete old objects - if (factorizationSolver_) - { - delete factorizationSolver_; - factorizationSolver_ = nullptr; - } - if (refactorizationSolver_) - { - delete refactorizationSolver_; - refactorizationSolver_ = nullptr; - } - if (preconditionerSolver_) - { - delete preconditionerSolver_; - preconditionerSolver_ = nullptr; - } - if (preconditioner_) - { - delete preconditioner_; - preconditioner_ = nullptr; - } - if (iterativeSolver_) - { - delete iterativeSolver_; - iterativeSolver_ = nullptr; - } - if (gs_) - { - delete gs_; - gs_ = nullptr; - } + iterativeSolver_.reset(); + preconditioner_.reset(); + preconditionerSolver_.reset(); + refactorizationSolver_.reset(); + factorizationSolver_.reset(); + gs_.reset(); // Create factorization solver if (factorizationMethod_ == "none") @@ -293,7 +231,7 @@ namespace ReSolve } else if (factorizationMethod_ == "klu") { - factorizationSolver_ = new ReSolve::LinSolverDirectKLU(); + factorizationSolver_.reset(new ReSolve::LinSolverDirectKLU()); #endif } else @@ -314,23 +252,23 @@ namespace ReSolve } else if (refactorizationMethod_ == "glu") { - refactorizationSolver_ = new ReSolve::LinSolverDirectCuSolverGLU(workspaceCuda_); + refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverGLU(workspaceCuda_)); } else if (refactorizationMethod_ == "cusolverrf") { - refactorizationSolver_ = new ReSolve::LinSolverDirectCuSolverRf(); + refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverRf()); #ifdef RESOLVE_USE_CUDSS } else if (refactorizationMethod_ == "cudssrf") { - refactorizationSolver_ = new ReSolve::LinSolverDirectCuDssRf(); + refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuDssRf()); #endif #endif #ifdef RESOLVE_USE_HIP } else if (refactorizationMethod_ == "rocsolverrf") { - refactorizationSolver_ = new ReSolve::LinSolverDirectRocSolverRf(workspaceHip_); + refactorizationSolver_.reset(new ReSolve::LinSolverDirectRocSolverRf(workspaceHip_)); #endif } else @@ -344,9 +282,9 @@ namespace ReSolve if (irMethod_ == "fgmres") { setGramSchmidtMethod(gsMethod_); - iterativeSolver_ = new LinSolverIterativeFGMRES(matrixHandler_, - vectorHandler_, - gs_); + iterativeSolver_.reset(new LinSolverIterativeFGMRES(matrixHandler_.get(), + vectorHandler_.get(), + gs_.get())); } // Create preconditioner @@ -358,21 +296,21 @@ namespace ReSolve { if (memspace_ == "cpu") { - preconditionerSolver_ = new LinSolverDirectCpuILU0(workspaceCpu_); - preconditioner_ = new PreconditionerLU(preconditionerSolver_); + preconditionerSolver_.reset(new LinSolverDirectCpuILU0(workspaceCpu_)); + preconditioner_.reset(new PreconditionerLU(preconditionerSolver_.get())); #ifdef RESOLVE_USE_CUDA } else if (memspace_ == "cuda") { - preconditionerSolver_ = new LinSolverDirectCuSparseILU0(workspaceCuda_); - preconditioner_ = new PreconditionerLU(preconditionerSolver_); + preconditionerSolver_.reset(new LinSolverDirectCuSparseILU0(workspaceCuda_)); + preconditioner_.reset(new PreconditionerLU(preconditionerSolver_.get())); #endif #ifdef RESOLVE_USE_HIP } else if (memspace_ == "hip") { - preconditionerSolver_ = new LinSolverDirectRocSparseILU0(workspaceHip_); - preconditioner_ = new PreconditionerLU(preconditionerSolver_); + preconditionerSolver_.reset(new LinSolverDirectRocSparseILU0(workspaceHip_)); + preconditioner_.reset(new PreconditionerLU(preconditionerSolver_.get())); #endif } else @@ -408,17 +346,17 @@ namespace ReSolve sketch = LinSolverIterativeRandFGMRES::cs; } setGramSchmidtMethod(gsMethod_); - iterativeSolver_ = new LinSolverIterativeRandFGMRES(matrixHandler_, - vectorHandler_, - sketch, - gs_); + iterativeSolver_.reset(new LinSolverIterativeRandFGMRES(matrixHandler_.get(), + vectorHandler_.get(), + sketch, + gs_.get())); } else if (solveMethod_ == "fgmres") { setGramSchmidtMethod(gsMethod_); - iterativeSolver_ = new LinSolverIterativeFGMRES(matrixHandler_, - vectorHandler_, - gs_); + iterativeSolver_.reset(new LinSolverIterativeFGMRES(matrixHandler_.get(), + vectorHandler_.get(), + gs_.get())); } else { @@ -516,7 +454,7 @@ namespace ReSolve { status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_); - LinSolverDirectCuSolverRf* Rf = dynamic_cast(refactorizationSolver_); + LinSolverDirectCuSolverRf* Rf = dynamic_cast(refactorizationSolver_.get()); Rf->setNumericalProperties(1e-14, 1e-1); is_solve_on_device_ = false; @@ -524,7 +462,7 @@ namespace ReSolve #ifdef RESOLVE_USE_CUDSS else if (refactorizationMethod_ == "cudssrf") { - LinSolverDirectCuDssRf* Rf = dynamic_cast(refactorizationSolver_); + LinSolverDirectCuDssRf* Rf = dynamic_cast(refactorizationSolver_.get()); Rf->setNumericalProperties(1e-14, 1e-1); status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_); @@ -538,7 +476,7 @@ namespace ReSolve if (refactorizationMethod_ == "rocsolverrf") { is_solve_on_device_ = false; - status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_, resVector_); + status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_, resVector_.get()); } #endif @@ -546,14 +484,8 @@ namespace ReSolve { status += iterativeSolver_->setup(A_); - if (preconditioner_) - { - delete preconditioner_; - preconditioner_ = nullptr; - } - - preconditioner_ = new PreconditionerLU(refactorizationSolver_); - status += iterativeSolver_->setPreconditioner(preconditioner_); + preconditioner_.reset(new PreconditionerLU(refactorizationSolver_.get())); + status += iterativeSolver_->setPreconditioner(preconditioner_.get()); } return status; } @@ -665,7 +597,7 @@ namespace ReSolve { is_solve_on_device_ = true; } - status += iterativeSolver_->setPreconditioner(preconditioner_); + status += iterativeSolver_->setPreconditioner(preconditioner_.get()); return status; } @@ -748,11 +680,7 @@ namespace ReSolve void SystemSolver::setRefactorizationMethod(std::string method) { refactorizationMethod_ = method; - if (refactorizationSolver_) - { - delete refactorizationSolver_; - refactorizationSolver_ = nullptr; - } + refactorizationSolver_.reset(); // Create refactorization solver if (refactorizationMethod_ == "klu") @@ -762,17 +690,17 @@ namespace ReSolve } else if (refactorizationMethod_ == "glu") { - refactorizationSolver_ = new ReSolve::LinSolverDirectCuSolverGLU(workspaceCuda_); + refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverGLU(workspaceCuda_)); } else if (refactorizationMethod_ == "cusolverrf") { - refactorizationSolver_ = new ReSolve::LinSolverDirectCuSolverRf(); + refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverRf()); #endif #ifdef RESOLVE_USE_HIP } else if (refactorizationMethod_ == "rocsolverrf") { - refactorizationSolver_ = new ReSolve::LinSolverDirectRocSolverRf(workspaceHip_); + refactorizationSolver_.reset(new ReSolve::LinSolverDirectRocSolverRf(workspaceHip_)); #endif } else @@ -794,11 +722,7 @@ namespace ReSolve // Remove existing iterative solver and set IR to "none". irMethod_ = "none"; - if (iterativeSolver_) - { - delete iterativeSolver_; - iterativeSolver_ = nullptr; - } + iterativeSolver_.reset(); if (method == "randgmres") { @@ -819,17 +743,17 @@ namespace ReSolve } setGramSchmidtMethod(gsMethod_); - iterativeSolver_ = new LinSolverIterativeRandFGMRES(matrixHandler_, - vectorHandler_, - sketch, - gs_); + iterativeSolver_.reset(new LinSolverIterativeRandFGMRES(matrixHandler_.get(), + vectorHandler_.get(), + sketch, + gs_.get())); } else if (solveMethod_ == "fgmres") { setGramSchmidtMethod(gsMethod_); - iterativeSolver_ = new LinSolverIterativeFGMRES(matrixHandler_, - vectorHandler_, - gs_); + iterativeSolver_.reset(new LinSolverIterativeFGMRES(matrixHandler_.get(), + vectorHandler_.get(), + gs_.get())); } else { @@ -850,17 +774,8 @@ namespace ReSolve */ void SystemSolver::setRefinementMethod(std::string method, std::string gsMethod) { - if (iterativeSolver_ != nullptr) - { - delete iterativeSolver_; - iterativeSolver_ = nullptr; - } - - if (gs_ != nullptr) - { - delete gs_; - gs_ = nullptr; - } + iterativeSolver_.reset(); + gs_.reset(); irMethod_ = "none"; @@ -887,10 +802,10 @@ namespace ReSolve if (method == "fgmres") { setGramSchmidtMethod(gsMethod); - iterativeSolver_ = new LinSolverIterativeFGMRES(matrixHandler_, - vectorHandler_, - gs_); - irMethod_ = method; + iterativeSolver_.reset(new LinSolverIterativeFGMRES(matrixHandler_.get(), + vectorHandler_.get(), + gs_.get())); + irMethod_ = method; } else { @@ -937,7 +852,7 @@ namespace ReSolve if (memspace_ == "cpu") { resVector_->copyFromExternal(rhs, memory::HOST, memory::HOST); - norm_b = std::sqrt(vectorHandler_->dot(resVector_, resVector_, memory::HOST)); + norm_b = std::sqrt(vectorHandler_->dot(resVector_.get(), resVector_.get(), memory::HOST)); #if defined(RESOLVE_USE_HIP) || defined(RESOLVE_USE_CUDA) } else if (memspace_ == "cuda" || memspace_ == "hip") @@ -945,13 +860,13 @@ namespace ReSolve if (is_solve_on_device_) { resVector_->copyFromExternal(rhs, memory::DEVICE, memory::DEVICE); - norm_b = std::sqrt(vectorHandler_->dot(resVector_, resVector_, memory::DEVICE)); + norm_b = std::sqrt(vectorHandler_->dot(resVector_.get(), resVector_.get(), memory::DEVICE)); } else { resVector_->copyFromExternal(rhs, memory::HOST, memory::DEVICE); resVector_->syncData(memory::HOST); - norm_b = std::sqrt(vectorHandler_->dot(resVector_, resVector_, memory::HOST)); + norm_b = std::sqrt(vectorHandler_->dot(resVector_.get(), resVector_.get(), memory::HOST)); // ms = memory::HOST; } ms = memory::DEVICE; @@ -963,8 +878,8 @@ namespace ReSolve return -1.0; } matrixHandler_->setValuesChanged(true, ms); - matrixHandler_->matvec(A_, x, resVector_, &ONE, &MINUS_ONE, ms); - resnorm = std::sqrt(vectorHandler_->dot(resVector_, resVector_, ms)); + matrixHandler_->matvec(A_, x, resVector_.get(), &ONE, &MINUS_ONE, ms); + resnorm = std::sqrt(vectorHandler_->dot(resVector_.get(), resVector_.get(), ms)); return resnorm / norm_b; } @@ -1000,8 +915,8 @@ namespace ReSolve return -1.0; } matrixHandler_->setValuesChanged(true, ms); - matrixHandler_->matvec(A_, x, resVector_, &ONE, &MINUS_ONE, ms); - resnorm = vectorHandler_->amax(resVector_, ms); + matrixHandler_->matvec(A_, x, resVector_.get(), &ONE, &MINUS_ONE, ms); + resnorm = vectorHandler_->amax(resVector_.get(), ms); norm_x = vectorHandler_->amax(x, ms); matrixHandler_->matrixInfNorm(A_, &norm_A, ms); return resnorm / (norm_x * norm_A); @@ -1073,7 +988,7 @@ namespace ReSolve if (iterativeSolver_) { // TODO: Use cast here as a temporary solution; will be replaced by parameter setting framework - auto* sol = dynamic_cast(iterativeSolver_); + auto* sol = dynamic_cast(iterativeSolver_.get()); sol->setSketchingMethod(tmp); } @@ -1121,7 +1036,7 @@ namespace ReSolve } else { - gs_ = new GramSchmidt(vectorHandler_, gs_variant); + gs_.reset(new GramSchmidt(vectorHandler_.get(), gs_variant)); } return 0; diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index 86dfc3a4d..f68e99099 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include @@ -95,20 +96,23 @@ namespace ReSolve int setGramSchmidtMethod(std::string gs_method); private: - LinSolverDirect* factorizationSolver_{nullptr}; - LinSolverDirect* refactorizationSolver_{nullptr}; - LinSolverDirect* preconditionerSolver_{nullptr}; - LinSolverIterative* iterativeSolver_{nullptr}; - GramSchmidt* gs_{nullptr}; - Preconditioner* preconditioner_{nullptr}; + std::unique_ptr matrixHandler_; + std::unique_ptr vectorHandler_; + + // Owned solver components. Declaration order matters for destruction: + // members are destroyed in reverse order, so dependents are listed + // after the objects they reference. + std::unique_ptr factorizationSolver_; + std::unique_ptr refactorizationSolver_; + std::unique_ptr preconditionerSolver_; + std::unique_ptr gs_; + std::unique_ptr preconditioner_; + std::unique_ptr iterativeSolver_; LinAlgWorkspaceCUDA* workspaceCuda_{nullptr}; LinAlgWorkspaceHIP* workspaceHip_{nullptr}; LinAlgWorkspaceCpu* workspaceCpu_{nullptr}; - MatrixHandler* matrixHandler_{nullptr}; - VectorHandler* vectorHandler_{nullptr}; - bool is_solve_on_device_{false}; matrix_type* L_{nullptr}; @@ -117,7 +121,7 @@ namespace ReSolve index_type* P_{nullptr}; index_type* Q_{nullptr}; - vector_type* resVector_{nullptr}; + std::unique_ptr resVector_; matrix::Sparse* A_{nullptr}; From 1ffd33491a4b796e1e05d56f863357e96b9e5e5d Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 16:57:44 -0400 Subject: [PATCH 02/16] Consolidate refactorization method configuration. --- resolve/SystemSolver.cpp | 125 ++++++++++++++++++++------------------- resolve/SystemSolver.hpp | 4 +- 2 files changed, 66 insertions(+), 63 deletions(-) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 1a571d0f5..6d4702595 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -241,40 +241,8 @@ namespace ReSolve } // Create refactorization solver - if (refactorizationMethod_ == "none") - { - // do nothing - } - else if (refactorizationMethod_ == "klu") - { - // do nothing for now, KLU is the only factorization solver available -#ifdef RESOLVE_USE_CUDA - } - else if (refactorizationMethod_ == "glu") - { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverGLU(workspaceCuda_)); - } - else if (refactorizationMethod_ == "cusolverrf") - { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverRf()); -#ifdef RESOLVE_USE_CUDSS - } - else if (refactorizationMethod_ == "cudssrf") - { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuDssRf()); -#endif -#endif -#ifdef RESOLVE_USE_HIP - } - else if (refactorizationMethod_ == "rocsolverrf") + if (createRefactorizationSolver() != 0) { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectRocSolverRf(workspaceHip_)); -#endif - } - else - { - out::error() << "Refactorization method " << refactorizationMethod_ - << " not recognized ...\n"; return 1; } @@ -673,41 +641,20 @@ namespace ReSolve * @param[in] method - ID for the refactorization method * * @post Destroys whatever refactorization solver existed before - * and sets `refactorization_solver_` pointer to the new + * and sets `refactorizationSolver_` pointer to the new * refactorization object. Sets refactorization method ID - * to the value in input parameter `method`. + * to the value in input parameter `method`. Resets + * `is_solve_on_device_` since the new solver has not been set up yet. + * + * @return int 0 if successful, 1 if method is not recognized */ - void SystemSolver::setRefactorizationMethod(std::string method) + int SystemSolver::setRefactorizationMethod(std::string method) { refactorizationMethod_ = method; refactorizationSolver_.reset(); + is_solve_on_device_ = false; - // Create refactorization solver - if (refactorizationMethod_ == "klu") - { - // do nothing for now -#ifdef RESOLVE_USE_CUDA - } - else if (refactorizationMethod_ == "glu") - { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverGLU(workspaceCuda_)); - } - else if (refactorizationMethod_ == "cusolverrf") - { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverRf()); -#endif -#ifdef RESOLVE_USE_HIP - } - else if (refactorizationMethod_ == "rocsolverrf") - { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectRocSolverRf(workspaceHip_)); -#endif - } - else - { - out::error() << "Refactorization method " << refactorizationMethod_ - << " not recognized ...\n"; - } + return createRefactorizationSolver(); } /** @@ -999,6 +946,60 @@ namespace ReSolve // Private methods // + /** + * @brief Instantiates refactorization solver selected by `refactorizationMethod_`. + * + * Shared by `initialize()` and `setRefactorizationMethod()` so the list of + * supported backends is maintained in one place. + * + * @pre `refactorizationSolver_` is null. + * @post `refactorizationSolver_` points to a new solver, or stays null for + * methods "none" and "klu" (KLU refactorization reuses the factorization + * solver). + * + * @return int 0 if successful, 1 if method is not recognized + */ + int SystemSolver::createRefactorizationSolver() + { + if (refactorizationMethod_ == "none") + { + // do nothing + } + else if (refactorizationMethod_ == "klu") + { + // do nothing for now, KLU is the only factorization solver available +#ifdef RESOLVE_USE_CUDA + } + else if (refactorizationMethod_ == "glu") + { + refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverGLU(workspaceCuda_)); + } + else if (refactorizationMethod_ == "cusolverrf") + { + refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverRf()); +#ifdef RESOLVE_USE_CUDSS + } + else if (refactorizationMethod_ == "cudssrf") + { + refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuDssRf()); +#endif +#endif +#ifdef RESOLVE_USE_HIP + } + else if (refactorizationMethod_ == "rocsolverrf") + { + refactorizationSolver_.reset(new ReSolve::LinSolverDirectRocSolverRf(workspaceHip_)); +#endif + } + else + { + out::error() << "Refactorization method " << refactorizationMethod_ + << " not recognized ...\n"; + return 1; + } + return 0; + } + int SystemSolver::setGramSchmidtMethod(std::string variant) { // Map string input to the Gram-Schmidt variant enum diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index f68e99099..ae20b2d90 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -89,13 +89,15 @@ namespace ReSolve // Set solver parameters void setFactorizationMethod(std::string method); - void setRefactorizationMethod(std::string method); + int setRefactorizationMethod(std::string method); int setSolveMethod(std::string method); void setRefinementMethod(std::string method, std::string gs = "cgs2"); int setSketchingMethod(std::string method); int setGramSchmidtMethod(std::string gs_method); private: + int createRefactorizationSolver(); + std::unique_ptr matrixHandler_; std::unique_ptr vectorHandler_; From 0fe11819e64ec345dd4b384d25a2a9256ff471d2 Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 17:22:06 -0400 Subject: [PATCH 03/16] Reuse code through helper functions. --- resolve/SystemSolver.cpp | 648 +++++++++++++++++++-------------------- resolve/SystemSolver.hpp | 42 +-- 2 files changed, 338 insertions(+), 352 deletions(-) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 6d4702595..040a899a2 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -44,124 +44,148 @@ namespace ReSolve // Create a shortcut name for Logger static class using out = io::Logger; - SystemSolver::SystemSolver(LinAlgWorkspaceCpu* workspaceCpu, - std::string factor, - std::string refactor, - std::string solve, - std::string precond, - std::string ir) - : workspaceCpu_(workspaceCpu), - factorizationMethod_(factor), - refactorizationMethod_(refactor), - solveMethod_(solve), - precondition_method_(precond), - irMethod_(ir) + // Helpers below map user-facing string IDs to solver enums. They are kept + // in an anonymous namespace (still inside ReSolve) so they have internal + // linkage: they are implementation details of SystemSolver, are not part + // of the library's exported symbols, and cannot collide with same-named + // functions in other translation units. If these mappings become useful + // elsewhere, move them to the classes that own the enums (GramSchmidt and + // LinSolverIterativeRandFGMRES) as public static methods. + namespace { - if ((refactor != "none") && (precond != "none")) + /// Maps string ID to the sketching method enum; warns and defaults to count sketch. + LinSolverIterativeRandFGMRES::SketchingMethod sketchingMethodFromString(const std::string& method) { - out::warning() << "Incorrect input: " - << "Refactorization and preconditioning cannot both be enabled.\n" - << "Setting both to 'none' ...\n"; - refactorizationMethod_ = "none"; - precondition_method_ = "none"; + if (method == "count") + { + return LinSolverIterativeRandFGMRES::cs; + } + if (method == "fwht") + { + return LinSolverIterativeRandFGMRES::fwht; + } + out::warning() << "Sketching method " << method << " not recognized!\n" + << "Using default (count sketch).\n"; + return LinSolverIterativeRandFGMRES::cs; } - if ((solve == "randgmres" || solve == "fgmres") && (ir != "none")) + /// Maps string ID to the Gram-Schmidt variant enum; warns and defaults to CGS2. + GramSchmidt::GSVariant gsVariantFromString(const std::string& variant) { - out::warning() << "Incorrect input: " - << "Iterative refinement cannot be enabled together with an iterative solve method.\n" - << "Setting refinement method to 'none' ...\n"; - irMethod_ = "none"; + if (variant == "cgs2") + { + return GramSchmidt::CGS2; + } + if (variant == "mgs") + { + return GramSchmidt::MGS; + } + if (variant == "mgs_two_sync") + { + return GramSchmidt::MGS_TWO_SYNC; + } + if (variant == "mgs_pm") + { + return GramSchmidt::MGS_PM; + } + if (variant == "cgs1") + { + return GramSchmidt::CGS1; + } + out::warning() << "Gram-Schmidt variant " << variant << " not recognized.\n"; + out::warning() << "Using default CGS2 Gram-Schmidt variant.\n"; + return GramSchmidt::CGS2; } - // Instantiate handlers - matrixHandler_.reset(new MatrixHandler(workspaceCpu_)); - vectorHandler_.reset(new VectorHandler(workspaceCpu_)); + /// Memory space ID associated with each workspace type. + const char* memorySpaceName(LinAlgWorkspaceCpu*) + { + return "cpu"; + } +#ifdef RESOLVE_USE_CUDA + const char* memorySpaceName(LinAlgWorkspaceCUDA*) + { + return "cuda"; + } +#endif +#ifdef RESOLVE_USE_HIP + const char* memorySpaceName(LinAlgWorkspaceHIP*) + { + return "hip"; + } +#endif + } // namespace - memspace_ = "cpu"; + /** + * @brief Shared tail of all constructors. + * + * Creates matrix and vector handlers for the given workspace, derives the + * memory space ID from the workspace type, validates the user-selected + * method combination and instantiates solver components. + * + * @tparam Workspace - one of LinAlgWorkspaceCpu, LinAlgWorkspaceCUDA, LinAlgWorkspaceHIP + * @param[in] workspace - pointer to the workspace (not owned) + */ + template + void SystemSolver::completeSetup(Workspace* workspace) + { + matrix_handler_.reset(new MatrixHandler(workspace)); + vector_handler_.reset(new VectorHandler(workspace)); + memspace_ = memorySpaceName(workspace); + validateConfiguration(); initialize(); } + SystemSolver::SystemSolver(LinAlgWorkspaceCpu* workspace_cpu, + std::string factor, + std::string refactor, + std::string solve, + std::string precond, + std::string ir) + : workspace_cpu_(workspace_cpu), + factorization_method_(factor), + refactorization_method_(refactor), + solve_method_(solve), + precondition_method_(precond), + ir_method_(ir) + { + completeSetup(workspace_cpu_); + } + #ifdef RESOLVE_USE_CUDA - SystemSolver::SystemSolver(LinAlgWorkspaceCUDA* workspaceCuda, + SystemSolver::SystemSolver(LinAlgWorkspaceCUDA* workspace_cuda, std::string factor, std::string refactor, std::string solve, std::string precond, std::string ir) - : workspaceCuda_(workspaceCuda), - factorizationMethod_(factor), - refactorizationMethod_(refactor), - solveMethod_(solve), + : workspace_cuda_(workspace_cuda), + factorization_method_(factor), + refactorization_method_(refactor), + solve_method_(solve), precondition_method_(precond), - irMethod_(ir) + ir_method_(ir) { - if ((refactor != "none") && (precond != "none")) - { - out::warning() << "Incorrect input: " - << "Refactorization and preconditioning cannot both be enabled.\n" - << "Setting both to 'none' ...\n"; - refactorizationMethod_ = "none"; - precondition_method_ = "none"; - } - - if ((solve == "randgmres" || solve == "fgmres") && (ir != "none")) - { - out::warning() << "Incorrect input: " - << "Iterative refinement cannot be enabled together with an iterative solve method.\n" - << "Setting refinement method to 'none' ...\n"; - irMethod_ = "none"; - } - - // Instantiate handlers - matrixHandler_.reset(new MatrixHandler(workspaceCuda_)); - vectorHandler_.reset(new VectorHandler(workspaceCuda_)); - - memspace_ = "cuda"; - - initialize(); + completeSetup(workspace_cuda_); } #endif #ifdef RESOLVE_USE_HIP - SystemSolver::SystemSolver(LinAlgWorkspaceHIP* workspaceHip, + SystemSolver::SystemSolver(LinAlgWorkspaceHIP* workspace_hip, std::string factor, std::string refactor, std::string solve, std::string precond, std::string ir) - : workspaceHip_(workspaceHip), - factorizationMethod_(factor), - refactorizationMethod_(refactor), - solveMethod_(solve), + : workspace_hip_(workspace_hip), + factorization_method_(factor), + refactorization_method_(refactor), + solve_method_(solve), precondition_method_(precond), - irMethod_(ir) + ir_method_(ir) { - if ((refactor != "none") && (precond != "none")) - { - out::warning() << "Incorrect input: " - << "Refactorization and preconditioning cannot both be enabled.\n" - << "Setting both to 'none' ...\n"; - refactorizationMethod_ = "none"; - precondition_method_ = "none"; - } - - if ((solve == "randgmres" || solve == "fgmres") && (ir != "none")) - { - out::warning() << "Incorrect input: " - << "Iterative refinement cannot be enabled together with an iterative solve method.\n" - << "Setting refinement method to 'none' ...\n"; - irMethod_ = "none"; - } - - // Instantiate handlers - matrixHandler_.reset(new MatrixHandler(workspaceHip_)); - vectorHandler_.reset(new VectorHandler(workspaceHip_)); - - memspace_ = "hip"; - - initialize(); + completeSetup(workspace_hip_); } #endif @@ -177,26 +201,26 @@ namespace ReSolve { int status = 0; A_ = A; - resVector_.reset(new vector_type(A->getNumRows())); + res_vector_.reset(new vector_type(A->getNumRows())); if (memspace_ == "cpu") { - resVector_->allocate(memory::HOST); + res_vector_->allocate(memory::HOST); } else { - resVector_->allocate(memory::DEVICE); - matrixHandler_->setValuesChanged(true, memory::DEVICE); + res_vector_->allocate(memory::DEVICE); + matrix_handler_->setValuesChanged(true, memory::DEVICE); } // If we use iterative solver, we can set it up here - if (solveMethod_ == "randgmres") + if (solve_method_ == "randgmres") { - auto* rgmres = dynamic_cast(iterativeSolver_.get()); + auto* rgmres = dynamic_cast(iterative_solver_.get()); status += rgmres->setup(A_); } - else if (solveMethod_ == "fgmres") + else if (solve_method_ == "fgmres") { - auto* fgmres = dynamic_cast(iterativeSolver_.get()); + auto* fgmres = dynamic_cast(iterative_solver_.get()); status += fgmres->setup(A_); } else @@ -216,27 +240,27 @@ namespace ReSolve int SystemSolver::initialize() { // First delete old objects - iterativeSolver_.reset(); + iterative_solver_.reset(); preconditioner_.reset(); - preconditionerSolver_.reset(); - refactorizationSolver_.reset(); - factorizationSolver_.reset(); + preconditioner_solver_.reset(); + refactorization_solver_.reset(); + factorization_solver_.reset(); gs_.reset(); // Create factorization solver - if (factorizationMethod_ == "none") + if (factorization_method_ == "none") { // do nothing #ifdef RESOLVE_USE_KLU } - else if (factorizationMethod_ == "klu") + else if (factorization_method_ == "klu") { - factorizationSolver_.reset(new ReSolve::LinSolverDirectKLU()); + factorization_solver_.reset(new ReSolve::LinSolverDirectKLU()); #endif } else { - out::error() << "Unrecognized factorization " << factorizationMethod_ << "\n"; + out::error() << "Unrecognized factorization " << factorization_method_ << "\n"; return 1; } @@ -247,12 +271,13 @@ namespace ReSolve } // Create iterative refinement - if (irMethod_ == "fgmres") + if (ir_method_ != "none") { - setGramSchmidtMethod(gsMethod_); - iterativeSolver_.reset(new LinSolverIterativeFGMRES(matrixHandler_.get(), - vectorHandler_.get(), - gs_.get())); + if (createIterativeSolver(ir_method_) != 0) + { + out::error() << "Iterative refinement method " << ir_method_ << " not recognized.\n"; + return 1; + } } // Create preconditioner @@ -264,21 +289,21 @@ namespace ReSolve { if (memspace_ == "cpu") { - preconditionerSolver_.reset(new LinSolverDirectCpuILU0(workspaceCpu_)); - preconditioner_.reset(new PreconditionerLU(preconditionerSolver_.get())); + preconditioner_solver_.reset(new LinSolverDirectCpuILU0(workspace_cpu_)); + preconditioner_.reset(new PreconditionerLU(preconditioner_solver_.get())); #ifdef RESOLVE_USE_CUDA } else if (memspace_ == "cuda") { - preconditionerSolver_.reset(new LinSolverDirectCuSparseILU0(workspaceCuda_)); - preconditioner_.reset(new PreconditionerLU(preconditionerSolver_.get())); + preconditioner_solver_.reset(new LinSolverDirectCuSparseILU0(workspace_cuda_)); + preconditioner_.reset(new PreconditionerLU(preconditioner_solver_.get())); #endif #ifdef RESOLVE_USE_HIP } else if (memspace_ == "hip") { - preconditionerSolver_.reset(new LinSolverDirectRocSparseILU0(workspaceHip_)); - preconditioner_.reset(new PreconditionerLU(preconditionerSolver_.get())); + preconditioner_solver_.reset(new LinSolverDirectRocSparseILU0(workspace_hip_)); + preconditioner_.reset(new PreconditionerLU(preconditioner_solver_.get())); #endif } else @@ -296,39 +321,9 @@ namespace ReSolve } // Create iterative solver - if (solveMethod_ == "randgmres") + if (solve_method_ == "randgmres" || solve_method_ == "fgmres") { - LinSolverIterativeRandFGMRES::SketchingMethod sketch; - if (sketching_method_ == "count") - { - sketch = LinSolverIterativeRandFGMRES::cs; - } - else if (sketching_method_ == "fwht") - { - sketch = LinSolverIterativeRandFGMRES::fwht; - } - else - { - out::warning() << "Sketching method " << sketching_method_ << " not recognized!\n" - << "Using default.\n"; - sketch = LinSolverIterativeRandFGMRES::cs; - } - setGramSchmidtMethod(gsMethod_); - iterativeSolver_.reset(new LinSolverIterativeRandFGMRES(matrixHandler_.get(), - vectorHandler_.get(), - sketch, - gs_.get())); - } - else if (solveMethod_ == "fgmres") - { - setGramSchmidtMethod(gsMethod_); - iterativeSolver_.reset(new LinSolverIterativeFGMRES(matrixHandler_.get(), - vectorHandler_.get(), - gs_.get())); - } - else - { - // do nothing + createIterativeSolver(solve_method_); } return 0; @@ -342,39 +337,39 @@ namespace ReSolve return 1; } - if (factorizationMethod_ == "klu") + if (factorization_method_ == "klu") { - factorizationSolver_->setup(A_); - return factorizationSolver_->analyze(); + factorization_solver_->setup(A_); + return factorization_solver_->analyze(); } return 1; } int SystemSolver::factorize() { - if (factorizationMethod_ == "klu") + if (factorization_method_ == "klu") { is_solve_on_device_ = false; - return factorizationSolver_->factorize(); + return factorization_solver_->factorize(); } return 1; } int SystemSolver::refactorize() { - if (refactorizationMethod_ == "klu") + if (refactorization_method_ == "klu") { - return factorizationSolver_->refactorize(); + return factorization_solver_->refactorize(); } - if (refactorizationMethod_ == "glu" || refactorizationMethod_ == "cusolverrf" || refactorizationMethod_ == "rocsolverrf" + if (refactorization_method_ == "glu" || refactorization_method_ == "cusolverrf" || refactorization_method_ == "rocsolverrf" #ifdef RESOLVE_USE_CUDSS - || refactorizationMethod_ == "cudssrf" + || refactorization_method_ == "cudssrf" #endif ) { is_solve_on_device_ = true; - return refactorizationSolver_->refactorize(); + return refactorization_solver_->refactorize(); } return 1; @@ -401,10 +396,10 @@ namespace ReSolve { int status = 0; - L_ = factorizationSolver_->getLFactor(); - U_ = factorizationSolver_->getUFactor(); - P_ = factorizationSolver_->getPOrdering(); - Q_ = factorizationSolver_->getQOrdering(); + L_ = factorization_solver_->getLFactor(); + U_ = factorization_solver_->getUFactor(); + P_ = factorization_solver_->getPOrdering(); + Q_ = factorization_solver_->getQOrdering(); if (L_ == nullptr) { @@ -413,27 +408,27 @@ namespace ReSolve } #ifdef RESOLVE_USE_CUDA - if (refactorizationMethod_ == "glu") + if (refactorization_method_ == "glu") { is_solve_on_device_ = true; - status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_); + status += refactorization_solver_->setup(A_, L_, U_, P_, Q_); } - else if (refactorizationMethod_ == "cusolverrf") + else if (refactorization_method_ == "cusolverrf") { - status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_); + status += refactorization_solver_->setup(A_, L_, U_, P_, Q_); - LinSolverDirectCuSolverRf* Rf = dynamic_cast(refactorizationSolver_.get()); + LinSolverDirectCuSolverRf* Rf = dynamic_cast(refactorization_solver_.get()); Rf->setNumericalProperties(1e-14, 1e-1); is_solve_on_device_ = false; } #ifdef RESOLVE_USE_CUDSS - else if (refactorizationMethod_ == "cudssrf") + else if (refactorization_method_ == "cudssrf") { - LinSolverDirectCuDssRf* Rf = dynamic_cast(refactorizationSolver_.get()); + LinSolverDirectCuDssRf* Rf = dynamic_cast(refactorization_solver_.get()); Rf->setNumericalProperties(1e-14, 1e-1); - status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_); + status += refactorization_solver_->setup(A_, L_, U_, P_, Q_); is_solve_on_device_ = false; } @@ -441,19 +436,19 @@ namespace ReSolve #endif #ifdef RESOLVE_USE_HIP - if (refactorizationMethod_ == "rocsolverrf") + if (refactorization_method_ == "rocsolverrf") { is_solve_on_device_ = false; - status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_, resVector_.get()); + status += refactorization_solver_->setup(A_, L_, U_, P_, Q_, res_vector_.get()); } #endif - if (irMethod_ == "fgmres") + if (ir_method_ == "fgmres") { - status += iterativeSolver_->setup(A_); + status += iterative_solver_->setup(A_); - preconditioner_.reset(new PreconditionerLU(refactorizationSolver_.get())); - status += iterativeSolver_->setPreconditioner(preconditioner_.get()); + preconditioner_.reset(new PreconditionerLU(refactorization_solver_.get())); + status += iterative_solver_->setPreconditioner(preconditioner_.get()); } return status; } @@ -476,35 +471,35 @@ namespace ReSolve int status = 0; // Use Krylov solver if selected - if (solveMethod_ == "randgmres" || solveMethod_ == "fgmres") + if (solve_method_ == "randgmres" || solve_method_ == "fgmres") { - status += iterativeSolver_->resetMatrix(A_); - status += iterativeSolver_->solve(rhs, x); + status += iterative_solver_->resetMatrix(A_); + status += iterative_solver_->solve(rhs, x); return status; } - if (solveMethod_ == "klu") + if (solve_method_ == "klu") { - status += factorizationSolver_->solve(rhs, x); + status += factorization_solver_->solve(rhs, x); } - if (solveMethod_ == "glu" || solveMethod_ == "cusolverrf" || solveMethod_ == "rocsolverrf" + if (solve_method_ == "glu" || solve_method_ == "cusolverrf" || solve_method_ == "rocsolverrf" #ifdef RESOLVE_USE_CUDSS - || solveMethod_ == "cudssrf" + || solve_method_ == "cudssrf" #endif ) { if (is_solve_on_device_) { - status += refactorizationSolver_->solve(rhs, x); + status += refactorization_solver_->solve(rhs, x); } else { - status += factorizationSolver_->solve(rhs, x); + status += factorization_solver_->solve(rhs, x); } } - if (irMethod_ == "fgmres") + if (ir_method_ == "fgmres") { if (is_solve_on_device_) { @@ -531,7 +526,7 @@ namespace ReSolve status += 1; } - if (iterativeSolver_ == nullptr) + if (iterative_solver_ == nullptr) { out::error() << "Iterative solver not initialized!\n"; status += 1; @@ -565,7 +560,7 @@ namespace ReSolve { is_solve_on_device_ = true; } - status += iterativeSolver_->setPreconditioner(preconditioner_.get()); + status += iterative_solver_->setPreconditioner(preconditioner_.get()); return status; } @@ -599,30 +594,30 @@ namespace ReSolve { int status = 0; - status += iterativeSolver_->resetMatrix(A_); - status += iterativeSolver_->solve(rhs, x); + status += iterative_solver_->resetMatrix(A_); + status += iterative_solver_->solve(rhs, x); return status; } LinSolverDirect& SystemSolver::getFactorizationSolver() { - return *factorizationSolver_; + return *factorization_solver_; } LinSolverDirect& SystemSolver::getRefactorizationSolver() { - return *refactorizationSolver_; + return *refactorization_solver_; } LinSolverDirect& SystemSolver::getPreconditionerSolver() { - return *preconditionerSolver_; + return *preconditioner_solver_; } LinSolverIterative& SystemSolver::getIterativeSolver() { - return *iterativeSolver_; + return *iterative_solver_; } Preconditioner& SystemSolver::getPreconditioner() @@ -632,7 +627,7 @@ namespace ReSolve void SystemSolver::setFactorizationMethod(std::string method) { - factorizationMethod_ = method; + factorization_method_ = method; } /** @@ -641,7 +636,7 @@ namespace ReSolve * @param[in] method - ID for the refactorization method * * @post Destroys whatever refactorization solver existed before - * and sets `refactorizationSolver_` pointer to the new + * and sets `refactorization_solver_` pointer to the new * refactorization object. Sets refactorization method ID * to the value in input parameter `method`. Resets * `is_solve_on_device_` since the new solver has not been set up yet. @@ -650,8 +645,8 @@ namespace ReSolve */ int SystemSolver::setRefactorizationMethod(std::string method) { - refactorizationMethod_ = method; - refactorizationSolver_.reset(); + refactorization_method_ = method; + refactorization_solver_.reset(); is_solve_on_device_ = false; return createRefactorizationSolver(); @@ -665,46 +660,15 @@ namespace ReSolve */ int SystemSolver::setSolveMethod(std::string method) { - solveMethod_ = method; + solve_method_ = method; // Remove existing iterative solver and set IR to "none". - irMethod_ = "none"; - iterativeSolver_.reset(); + ir_method_ = "none"; + iterative_solver_.reset(); - if (method == "randgmres") + if (createIterativeSolver(method) != 0) { - LinSolverIterativeRandFGMRES::SketchingMethod sketch; - if (sketching_method_ == "count") - { - sketch = LinSolverIterativeRandFGMRES::cs; - } - else if (sketching_method_ == "fwht") - { - sketch = LinSolverIterativeRandFGMRES::fwht; - } - else - { - out::warning() << "Sketching method " << sketching_method_ << " not recognized!\n" - << "Using default.\n"; - sketch = LinSolverIterativeRandFGMRES::cs; - } - - setGramSchmidtMethod(gsMethod_); - iterativeSolver_.reset(new LinSolverIterativeRandFGMRES(matrixHandler_.get(), - vectorHandler_.get(), - sketch, - gs_.get())); - } - else if (solveMethod_ == "fgmres") - { - setGramSchmidtMethod(gsMethod_); - iterativeSolver_.reset(new LinSolverIterativeFGMRES(matrixHandler_.get(), - vectorHandler_.get(), - gs_.get())); - } - else - { - out::error() << "Solve method " << solveMethod_ + out::error() << "Solve method " << solve_method_ << " not recognized ...\n"; return 1; } @@ -715,21 +679,21 @@ namespace ReSolve * @brief Sets iterative refinement method and related orthogonalization. * * @param[in] method - string ID for the iterative refinement method - * @param[in] gsMethod - string ID for the orthogonalization method to be used + * @param[in] gs_method - string ID for the orthogonalization method to be used * * @todo Iterative refinement temporarily disabled on CPU. Need to fix that. */ - void SystemSolver::setRefinementMethod(std::string method, std::string gsMethod) + void SystemSolver::setRefinementMethod(std::string method, std::string gs_method) { - iterativeSolver_.reset(); + iterative_solver_.reset(); gs_.reset(); - irMethod_ = "none"; + ir_method_ = "none"; if (method == "none") return; - if (solveMethod_ == "randgmres" || solveMethod_ == "fgmres") + if (solve_method_ == "randgmres" || solve_method_ == "fgmres") { out::warning() << "Iterative refinement cannot be enabled together with an " << "iterative solve method ('randgmres' or 'fgmres'). " @@ -744,15 +708,11 @@ namespace ReSolve return; } - gsMethod_ = gsMethod; + gs_method_ = gs_method; - if (method == "fgmres") + if (method == "fgmres" && createIterativeSolver(method) == 0) { - setGramSchmidtMethod(gsMethod); - iterativeSolver_.reset(new LinSolverIterativeFGMRES(matrixHandler_.get(), - vectorHandler_.get(), - gs_.get())); - irMethod_ = method; + ir_method_ = method; } else { @@ -766,18 +726,18 @@ namespace ReSolve real_type norm_b = 0.0; if (memspace_ == "cpu") { - norm_b = std::sqrt(vectorHandler_->dot(rhs, rhs, memory::HOST)); + norm_b = std::sqrt(vector_handler_->dot(rhs, rhs, memory::HOST)); #if defined(RESOLVE_USE_HIP) || defined(RESOLVE_USE_CUDA) } else if (memspace_ == "cuda" || memspace_ == "hip") { if (is_solve_on_device_) { - norm_b = std::sqrt(vectorHandler_->dot(rhs, rhs, memory::DEVICE)); + norm_b = std::sqrt(vector_handler_->dot(rhs, rhs, memory::DEVICE)); } else { - norm_b = std::sqrt(vectorHandler_->dot(rhs, rhs, memory::HOST)); + norm_b = std::sqrt(vector_handler_->dot(rhs, rhs, memory::HOST)); } #endif } @@ -792,28 +752,28 @@ namespace ReSolve real_type SystemSolver::getResidualNorm(vector_type* rhs, vector_type* x) { using namespace ReSolve::constants; - assert(rhs->getSize() == resVector_->getSize()); + assert(rhs->getSize() == res_vector_->getSize()); real_type norm_b = 0.0; real_type resnorm = 0.0; memory::MemorySpace ms = memory::HOST; if (memspace_ == "cpu") { - resVector_->copyFromExternal(rhs, memory::HOST, memory::HOST); - norm_b = std::sqrt(vectorHandler_->dot(resVector_.get(), resVector_.get(), memory::HOST)); + res_vector_->copyFromExternal(rhs, memory::HOST, memory::HOST); + norm_b = std::sqrt(vector_handler_->dot(res_vector_.get(), res_vector_.get(), memory::HOST)); #if defined(RESOLVE_USE_HIP) || defined(RESOLVE_USE_CUDA) } else if (memspace_ == "cuda" || memspace_ == "hip") { if (is_solve_on_device_) { - resVector_->copyFromExternal(rhs, memory::DEVICE, memory::DEVICE); - norm_b = std::sqrt(vectorHandler_->dot(resVector_.get(), resVector_.get(), memory::DEVICE)); + res_vector_->copyFromExternal(rhs, memory::DEVICE, memory::DEVICE); + norm_b = std::sqrt(vector_handler_->dot(res_vector_.get(), res_vector_.get(), memory::DEVICE)); } else { - resVector_->copyFromExternal(rhs, memory::HOST, memory::DEVICE); - resVector_->syncData(memory::HOST); - norm_b = std::sqrt(vectorHandler_->dot(resVector_.get(), resVector_.get(), memory::HOST)); + res_vector_->copyFromExternal(rhs, memory::HOST, memory::DEVICE); + res_vector_->syncData(memory::HOST); + norm_b = std::sqrt(vector_handler_->dot(res_vector_.get(), res_vector_.get(), memory::HOST)); // ms = memory::HOST; } ms = memory::DEVICE; @@ -824,34 +784,34 @@ namespace ReSolve out::error() << "Unrecognized device " << memspace_ << "\n"; return -1.0; } - matrixHandler_->setValuesChanged(true, ms); - matrixHandler_->matvec(A_, x, resVector_.get(), &ONE, &MINUS_ONE, ms); - resnorm = std::sqrt(vectorHandler_->dot(resVector_.get(), resVector_.get(), ms)); + matrix_handler_->setValuesChanged(true, ms); + matrix_handler_->matvec(A_, x, res_vector_.get(), &ONE, &MINUS_ONE, ms); + resnorm = std::sqrt(vector_handler_->dot(res_vector_.get(), res_vector_.get(), ms)); return resnorm / norm_b; } real_type SystemSolver::getNormOfScaledResiduals(vector_type* rhs, vector_type* x) { using namespace ReSolve::constants; - assert(rhs->getSize() == resVector_->getSize()); + assert(rhs->getSize() == res_vector_->getSize()); real_type norm_x = 0.0; real_type norm_A = 0.0; real_type resnorm = 0.0; memory::MemorySpace ms = memory::HOST; if (memspace_ == "cpu") { - resVector_->copyFromExternal(rhs, memory::HOST, memory::HOST); + res_vector_->copyFromExternal(rhs, memory::HOST, memory::HOST); #if defined(RESOLVE_USE_HIP) || defined(RESOLVE_USE_CUDA) } else if (memspace_ == "cuda" || memspace_ == "hip") { if (is_solve_on_device_) { - resVector_->copyFromExternal(rhs, memory::DEVICE, memory::DEVICE); + res_vector_->copyFromExternal(rhs, memory::DEVICE, memory::DEVICE); } else { - resVector_->copyFromExternal(rhs, memory::HOST, memory::DEVICE); + res_vector_->copyFromExternal(rhs, memory::HOST, memory::DEVICE); } ms = memory::DEVICE; #endif @@ -861,37 +821,37 @@ namespace ReSolve out::error() << "Unrecognized device " << memspace_ << "\n"; return -1.0; } - matrixHandler_->setValuesChanged(true, ms); - matrixHandler_->matvec(A_, x, resVector_.get(), &ONE, &MINUS_ONE, ms); - resnorm = vectorHandler_->amax(resVector_.get(), ms); - norm_x = vectorHandler_->amax(x, ms); - matrixHandler_->matrixInfNorm(A_, &norm_A, ms); + matrix_handler_->setValuesChanged(true, ms); + matrix_handler_->matvec(A_, x, res_vector_.get(), &ONE, &MINUS_ONE, ms); + resnorm = vector_handler_->amax(res_vector_.get(), ms); + norm_x = vector_handler_->amax(x, ms); + matrix_handler_->matrixInfNorm(A_, &norm_A, ms); return resnorm / (norm_x * norm_A); } const std::string SystemSolver::getFactorizationMethod() const { - return factorizationMethod_; + return factorization_method_; } const std::string SystemSolver::getRefactorizationMethod() const { - return refactorizationMethod_; + return refactorization_method_; } const std::string SystemSolver::getSolveMethod() const { - return solveMethod_; + return solve_method_; } const std::string SystemSolver::getRefinementMethod() const { - return irMethod_; + return ir_method_; } const std::string SystemSolver::getOrthogonalizationMethod() const { - return gsMethod_; + return gs_method_; } /** @@ -906,36 +866,22 @@ namespace ReSolve */ int SystemSolver::setSketchingMethod(std::string sketching_method) { - if (solveMethod_ != "randgmres") + if (solve_method_ != "randgmres") { out::warning() << "Trying to set sketching method to an incompatible solver.\n"; out::warning() << "The setting will be ignored.\n"; return 1; } - LinSolverIterativeRandFGMRES::SketchingMethod tmp; - if (sketching_method == "count") - { - tmp = LinSolverIterativeRandFGMRES::cs; - } - else if (sketching_method == "fwht") - { - tmp = LinSolverIterativeRandFGMRES::fwht; - } - else - { - out::warning() << "Sketching method " << sketching_method << " not recognized!\n" - << "Using default (count sketch).\n"; - tmp = LinSolverIterativeRandFGMRES::cs; - } + LinSolverIterativeRandFGMRES::SketchingMethod tmp = sketchingMethodFromString(sketching_method); sketching_method_ = sketching_method; // At this point iterative solver, if created, can only be LinSolverIterativeRandFGMRES - if (iterativeSolver_) + if (iterative_solver_) { // TODO: Use cast here as a temporary solution; will be replaced by parameter setting framework - auto* sol = dynamic_cast(iterativeSolver_.get()); + auto* sol = dynamic_cast(iterative_solver_.get()); sol->setSketchingMethod(tmp); } @@ -947,13 +893,76 @@ namespace ReSolve // /** - * @brief Instantiates refactorization solver selected by `refactorizationMethod_`. + * @brief Checks for unsupported combinations of user-selected methods. + * + * Refactorization and preconditioning are mutually exclusive, as are + * iterative refinement and an iterative solve method. Offending settings + * are reset to "none" with a warning. + */ + void SystemSolver::validateConfiguration() + { + if ((refactorization_method_ != "none") && (precondition_method_ != "none")) + { + out::warning() << "Incorrect input: " + << "Refactorization and preconditioning cannot both be enabled.\n" + << "Setting both to 'none' ...\n"; + refactorization_method_ = "none"; + precondition_method_ = "none"; + } + + if ((solve_method_ == "randgmres" || solve_method_ == "fgmres") && (ir_method_ != "none")) + { + out::warning() << "Incorrect input: " + << "Iterative refinement cannot be enabled together with an iterative solve method.\n" + << "Setting refinement method to 'none' ...\n"; + ir_method_ = "none"; + } + } + + /** + * @brief Instantiates Krylov solver of the requested type. + * + * Used both for iterative solve methods and for iterative refinement. + * Configures Gram-Schmidt orthogonalization according to `gs_method_` + * and, for randomized GMRES, sketching according to `sketching_method_`. + * + * @param[in] method - "fgmres" or "randgmres" + * @post `iterative_solver_` points to a new solver. + * + * @return int 0 if successful, 1 if method is not recognized + */ + int SystemSolver::createIterativeSolver(const std::string& method) + { + if (method == "randgmres") + { + setGramSchmidtMethod(gs_method_); + iterative_solver_.reset(new LinSolverIterativeRandFGMRES(matrix_handler_.get(), + vector_handler_.get(), + sketchingMethodFromString(sketching_method_), + gs_.get())); + } + else if (method == "fgmres") + { + setGramSchmidtMethod(gs_method_); + iterative_solver_.reset(new LinSolverIterativeFGMRES(matrix_handler_.get(), + vector_handler_.get(), + gs_.get())); + } + else + { + return 1; + } + return 0; + } + + /** + * @brief Instantiates refactorization solver selected by `refactorization_method_`. * * Shared by `initialize()` and `setRefactorizationMethod()` so the list of * supported backends is maintained in one place. * - * @pre `refactorizationSolver_` is null. - * @post `refactorizationSolver_` points to a new solver, or stays null for + * @pre `refactorization_solver_` is null. + * @post `refactorization_solver_` points to a new solver, or stays null for * methods "none" and "klu" (KLU refactorization reuses the factorization * solver). * @@ -961,39 +970,39 @@ namespace ReSolve */ int SystemSolver::createRefactorizationSolver() { - if (refactorizationMethod_ == "none") + if (refactorization_method_ == "none") { // do nothing } - else if (refactorizationMethod_ == "klu") + else if (refactorization_method_ == "klu") { // do nothing for now, KLU is the only factorization solver available #ifdef RESOLVE_USE_CUDA } - else if (refactorizationMethod_ == "glu") + else if (refactorization_method_ == "glu") { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverGLU(workspaceCuda_)); + refactorization_solver_.reset(new ReSolve::LinSolverDirectCuSolverGLU(workspace_cuda_)); } - else if (refactorizationMethod_ == "cusolverrf") + else if (refactorization_method_ == "cusolverrf") { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuSolverRf()); + refactorization_solver_.reset(new ReSolve::LinSolverDirectCuSolverRf()); #ifdef RESOLVE_USE_CUDSS } - else if (refactorizationMethod_ == "cudssrf") + else if (refactorization_method_ == "cudssrf") { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectCuDssRf()); + refactorization_solver_.reset(new ReSolve::LinSolverDirectCuDssRf()); #endif #endif #ifdef RESOLVE_USE_HIP } - else if (refactorizationMethod_ == "rocsolverrf") + else if (refactorization_method_ == "rocsolverrf") { - refactorizationSolver_.reset(new ReSolve::LinSolverDirectRocSolverRf(workspaceHip_)); + refactorization_solver_.reset(new ReSolve::LinSolverDirectRocSolverRf(workspace_hip_)); #endif } else { - out::error() << "Refactorization method " << refactorizationMethod_ + out::error() << "Refactorization method " << refactorization_method_ << " not recognized ...\n"; return 1; } @@ -1002,34 +1011,7 @@ namespace ReSolve int SystemSolver::setGramSchmidtMethod(std::string variant) { - // Map string input to the Gram-Schmidt variant enum - GramSchmidt::GSVariant gs_variant; - if (variant == "cgs2") - { - gs_variant = GramSchmidt::CGS2; - } - else if (variant == "mgs") - { - gs_variant = GramSchmidt::MGS; - } - else if (variant == "mgs_two_sync") - { - gs_variant = GramSchmidt::MGS_TWO_SYNC; - } - else if (variant == "mgs_pm") - { - gs_variant = GramSchmidt::MGS_PM; - } - else if (variant == "cgs1") - { - gs_variant = GramSchmidt::CGS1; - } - else - { - out::warning() << "Gram-Schmidt variant " << variant << " not recognized.\n"; - out::warning() << "Using default CGS2 Gram-Schmidt variant.\n"; - gs_variant = GramSchmidt::CGS2; - } + GramSchmidt::GSVariant gs_variant = gsVariantFromString(variant); if (gs_) { @@ -1037,7 +1019,7 @@ namespace ReSolve } else { - gs_.reset(new GramSchmidt(vectorHandler_.get(), gs_variant)); + gs_.reset(new GramSchmidt(vector_handler_.get(), gs_variant)); } return 0; diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index ae20b2d90..49ababaa9 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -35,19 +35,19 @@ namespace ReSolve using vector_type = vector::Vector; using matrix_type = matrix::Sparse; - SystemSolver(LinAlgWorkspaceCpu* workspaceCpu, + SystemSolver(LinAlgWorkspaceCpu* workspace_cpu, std::string factor = "klu", std::string refactor = "klu", std::string solve = "klu", std::string precond = "none", std::string ir = "none"); - SystemSolver(LinAlgWorkspaceCUDA* workspaceCuda, + SystemSolver(LinAlgWorkspaceCUDA* workspace_cuda, std::string factor = "klu", std::string refactor = "cusolverrf", std::string solve = "cusolverrf", std::string precond = "none", std::string ir = "none"); - SystemSolver(LinAlgWorkspaceHIP* workspaceHip, + SystemSolver(LinAlgWorkspaceHIP* workspace_hip, std::string factor = "klu", std::string refactor = "rocsolverrf", std::string solve = "rocsolverrf", @@ -96,24 +96,28 @@ namespace ReSolve int setGramSchmidtMethod(std::string gs_method); private: - int createRefactorizationSolver(); + template + void completeSetup(Workspace* workspace); + void validateConfiguration(); + int createRefactorizationSolver(); + int createIterativeSolver(const std::string& method); - std::unique_ptr matrixHandler_; - std::unique_ptr vectorHandler_; + std::unique_ptr matrix_handler_; + std::unique_ptr vector_handler_; // Owned solver components. Declaration order matters for destruction: // members are destroyed in reverse order, so dependents are listed // after the objects they reference. - std::unique_ptr factorizationSolver_; - std::unique_ptr refactorizationSolver_; - std::unique_ptr preconditionerSolver_; + std::unique_ptr factorization_solver_; + std::unique_ptr refactorization_solver_; + std::unique_ptr preconditioner_solver_; std::unique_ptr gs_; std::unique_ptr preconditioner_; - std::unique_ptr iterativeSolver_; + std::unique_ptr iterative_solver_; - LinAlgWorkspaceCUDA* workspaceCuda_{nullptr}; - LinAlgWorkspaceHIP* workspaceHip_{nullptr}; - LinAlgWorkspaceCpu* workspaceCpu_{nullptr}; + LinAlgWorkspaceCUDA* workspace_cuda_{nullptr}; + LinAlgWorkspaceHIP* workspace_hip_{nullptr}; + LinAlgWorkspaceCpu* workspace_cpu_{nullptr}; bool is_solve_on_device_{false}; @@ -123,17 +127,17 @@ namespace ReSolve index_type* P_{nullptr}; index_type* Q_{nullptr}; - std::unique_ptr resVector_; + std::unique_ptr res_vector_; matrix::Sparse* A_{nullptr}; // Configuration parameters - std::string factorizationMethod_{"none"}; - std::string refactorizationMethod_{"none"}; - std::string solveMethod_{"none"}; + std::string factorization_method_{"none"}; + std::string refactorization_method_{"none"}; + std::string solve_method_{"none"}; std::string precondition_method_{"none"}; - std::string irMethod_{"none"}; - std::string gsMethod_{"cgs2"}; + std::string ir_method_{"none"}; + std::string gs_method_{"cgs2"}; std::string sketching_method_{"count"}; ///< @todo move this to LinSolverIterative class std::string memspace_; From c9252719cbe78f26d22614f52ab528d1ba82053d Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 17:40:40 -0400 Subject: [PATCH 04/16] Configuration distinguishes between iterative and direct solvers. --- resolve/SystemSolver.cpp | 94 ++++++++++++++++++++++++++++------------ 1 file changed, 66 insertions(+), 28 deletions(-) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 040a899a2..068de0018 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -120,8 +120,8 @@ namespace ReSolve * @brief Shared tail of all constructors. * * Creates matrix and vector handlers for the given workspace, derives the - * memory space ID from the workspace type, validates the user-selected - * method combination and instantiates solver components. + * memory space ID from the workspace type and instantiates solver + * components via `initialize()`. * * @tparam Workspace - one of LinAlgWorkspaceCpu, LinAlgWorkspaceCUDA, LinAlgWorkspaceHIP * @param[in] workspace - pointer to the workspace (not owned) @@ -133,7 +133,6 @@ namespace ReSolve vector_handler_.reset(new VectorHandler(workspace)); memspace_ = memorySpaceName(workspace); - validateConfiguration(); initialize(); } @@ -234,11 +233,16 @@ namespace ReSolve /** * @brief Sets up the system solver * - * This method instantiates components of the system solver based on - * user inputs. + * This method validates the selected method combination (see + * `validateConfiguration()`) and then instantiates components of the + * system solver based on user inputs. It can be called again after the + * configuration has been changed through setters. */ int SystemSolver::initialize() { + // Make sure the method combination is consistent before creating objects + validateConfiguration(); + // First delete old objects iterative_solver_.reset(); preconditioner_.reset(); @@ -685,19 +689,27 @@ namespace ReSolve */ void SystemSolver::setRefinementMethod(std::string method, std::string gs_method) { + // 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 (method != "none") + { + out::warning() << "Iterative refinement cannot be enabled together with an " + << "iterative solve method ('randgmres' or 'fgmres'). " + << "Keeping refinement method 'none'.\n"; + } + ir_method_ = "none"; + return; + } + + // Direct solve path: the Krylov solver, if any, is the refinement solver. iterative_solver_.reset(); gs_.reset(); - ir_method_ = "none"; if (method == "none") - return; - - if (solve_method_ == "randgmres" || solve_method_ == "fgmres") { - out::warning() << "Iterative refinement cannot be enabled together with an " - << "iterative solve method ('randgmres' or 'fgmres'). " - << "Keeping refinement method 'none'.\n"; return; } @@ -893,29 +905,55 @@ namespace ReSolve // /** - * @brief Checks for unsupported combinations of user-selected methods. + * @brief Enforces supported combinations of user-selected methods. + * + * Two configurations are supported: + * - Iterative solver ("fgmres" or "randgmres" as the solve method): + * preconditioner may be set; factorization, refactorization and + * iterative refinement must be "none". + * - Direct solver (any other solve method): factorization and + * refactorization are set, iterative refinement is optional, and the + * preconditioner must be "none". * - * Refactorization and preconditioning are mutually exclusive, as are - * iterative refinement and an iterative solve method. Offending settings - * are reset to "none" with a warning. + * Offending settings are reset to "none" with a warning. */ void SystemSolver::validateConfiguration() { - if ((refactorization_method_ != "none") && (precondition_method_ != "none")) + const bool is_iterative_solve = (solve_method_ == "fgmres" || solve_method_ == "randgmres"); + + if (is_iterative_solve) { - out::warning() << "Incorrect input: " - << "Refactorization and preconditioning cannot both be enabled.\n" - << "Setting both to 'none' ...\n"; - refactorization_method_ = "none"; - precondition_method_ = "none"; + if (factorization_method_ != "none") + { + out::warning() << "Incorrect input: factorization method '" << factorization_method_ + << "' cannot be used with iterative solve method '" << solve_method_ + << "'. Setting factorization to 'none' ...\n"; + factorization_method_ = "none"; + } + if (refactorization_method_ != "none") + { + out::warning() << "Incorrect input: refactorization method '" << refactorization_method_ + << "' cannot be used with iterative solve method '" << solve_method_ + << "'. Setting refactorization to 'none' ...\n"; + refactorization_method_ = "none"; + } + if (ir_method_ != "none") + { + out::warning() << "Incorrect input: iterative refinement cannot be enabled " + << "together with iterative solve method '" << solve_method_ + << "'. Setting refinement method to 'none' ...\n"; + ir_method_ = "none"; + } } - - if ((solve_method_ == "randgmres" || solve_method_ == "fgmres") && (ir_method_ != "none")) + else { - out::warning() << "Incorrect input: " - << "Iterative refinement cannot be enabled together with an iterative solve method.\n" - << "Setting refinement method to 'none' ...\n"; - ir_method_ = "none"; + if (precondition_method_ != "none") + { + out::warning() << "Incorrect input: preconditioner '" << precondition_method_ + << "' can only be used with an iterative solve method ('fgmres' or 'randgmres'). " + << "Setting preconditioner to 'none' ...\n"; + precondition_method_ = "none"; + } } } From c27f63df8a48661a61b0076814c1dcae67f31e9d Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 17:46:32 -0400 Subject: [PATCH 05/16] Fix possibly dangling preconditioner pointer. --- resolve/SystemSolver.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 068de0018..025ba8c6e 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -644,12 +644,22 @@ namespace ReSolve * refactorization object. Sets refactorization method ID * to the value in input parameter `method`. Resets * `is_solve_on_device_` since the new solver has not been set up yet. + * If iterative refinement is enabled, its preconditioner wraps the + * refactorization solver and is destroyed as well; `refactorizationSetup()` + * recreates it. * * @return int 0 if successful, 1 if method is not recognized */ int SystemSolver::setRefactorizationMethod(std::string method) { refactorization_method_ = method; + + // Iterative refinement preconditioner holds a non-owning pointer to the + // refactorization solver, so it must not outlive it. + if (ir_method_ != "none") + { + preconditioner_.reset(); + } refactorization_solver_.reset(); is_solve_on_device_ = false; From d7045e7ac452481ad1ef28120aad96ca40f549ad Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 18:01:14 -0400 Subject: [PATCH 06/16] Fix issue in GS setting. --- resolve/SystemSolver.cpp | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 025ba8c6e..845d5e43a 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -97,6 +97,25 @@ namespace ReSolve return GramSchmidt::CGS2; } + /// Canonical string ID for a Gram-Schmidt variant (inverse of gsVariantFromString). + const char* gsVariantName(GramSchmidt::GSVariant variant) + { + switch (variant) + { + case GramSchmidt::CGS2: + return "cgs2"; + case GramSchmidt::MGS: + return "mgs"; + case GramSchmidt::MGS_TWO_SYNC: + return "mgs_two_sync"; + case GramSchmidt::MGS_PM: + return "mgs_pm"; + case GramSchmidt::CGS1: + return "cgs1"; + } + return "cgs2"; + } + /// Memory space ID associated with each workspace type. const char* memorySpaceName(LinAlgWorkspaceCpu*) { @@ -1057,10 +1076,26 @@ namespace ReSolve return 0; } + /** + * @brief Sets Gram-Schmidt orthogonalization variant. + * + * Records the variant in `gs_method_` so that it survives re-creation of + * the Krylov solver, and applies it to the existing `GramSchmidt` object + * or creates one if none exists yet. An unrecognized string ID falls back + * to CGS2 with a warning, and `gs_method_` is set to "cgs2" accordingly. + * + * @param[in] variant - string ID of the Gram-Schmidt variant + * + * @return int 0 on success + */ int SystemSolver::setGramSchmidtMethod(std::string variant) { GramSchmidt::GSVariant gs_variant = gsVariantFromString(variant); + // Store the canonical name so that the stored ID always matches the + // variant actually in use. + gs_method_ = gsVariantName(gs_variant); + if (gs_) { gs_->setVariant(gs_variant); From 38136be5a208c3b3c2e952cebbd040a17eef8a2a Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 18:25:31 -0400 Subject: [PATCH 07/16] Make KLU configuration more robust. --- resolve/SystemSolver.cpp | 163 ++++++++++++++++++++++++++++----------- resolve/SystemSolver.hpp | 8 +- 2 files changed, 119 insertions(+), 52 deletions(-) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 845d5e43a..8cc89c7ed 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -152,7 +152,15 @@ namespace ReSolve vector_handler_.reset(new VectorHandler(workspace)); memspace_ = memorySpaceName(workspace); - initialize(); + if (initialize() != 0) + { + out::error() << "SystemSolver initialization failed with factorization '" << factorization_method_ + << "', refactorization '" << refactorization_method_ + << "', solve '" << solve_method_ + << "', preconditioner '" << precondition_method_ + << "', iterative refinement '" << ir_method_ + << "'. Solver is not usable in this state.\n"; + } } SystemSolver::SystemSolver(LinAlgWorkspaceCpu* workspace_cpu, @@ -271,19 +279,8 @@ namespace ReSolve gs_.reset(); // Create factorization solver - if (factorization_method_ == "none") - { - // do nothing -#ifdef RESOLVE_USE_KLU - } - else if (factorization_method_ == "klu") + if (createFactorizationSolver() != 0) { - factorization_solver_.reset(new ReSolve::LinSolverDirectKLU()); -#endif - } - else - { - out::error() << "Unrecognized factorization " << factorization_method_ << "\n"; return 1; } @@ -648,9 +645,48 @@ namespace ReSolve return *preconditioner_; } - void SystemSolver::setFactorizationMethod(std::string method) + /** + * @brief Sets factorization method to use + * + * @param[in] method - ID for the factorization method + * + * @post Destroys the existing factorization solver together with the + * factors and permutation vectors it owned, and creates a new solver + * of the requested type. If refactorization method is "klu", the + * refactorization path reuses the factorization solver, so + * `refactorization_solver_` is reset as well. The caller must run + * `analyze()`, `factorize()` and, if applicable, `refactorizationSetup()` + * again before the next solve. + * + * @return int 0 if successful, 1 if method is not recognized or an + * iterative solve method is active + */ + int SystemSolver::setFactorizationMethod(std::string method) { + if (solve_method_ == "fgmres" || solve_method_ == "randgmres") + { + out::error() << "Factorization method cannot be set while iterative solve method '" + << solve_method_ << "' is active. Keeping '" << factorization_method_ << "'.\n"; + return 1; + } + factorization_method_ = method; + factorization_solver_.reset(); + + // Factors and permutations were owned by the old solver. + L_ = nullptr; + U_ = nullptr; + P_ = nullptr; + Q_ = nullptr; + is_solve_on_device_ = false; + + // KLU refactorization reuses the factorization solver; drop any stale state. + if (refactorization_method_ == "klu") + { + refactorization_solver_.reset(); + } + + return createFactorizationSolver(); } /** @@ -890,7 +926,7 @@ namespace ReSolve return ir_method_; } - const std::string SystemSolver::getOrthogonalizationMethod() const + const std::string SystemSolver::getGramSchmidtMethod() const { return gs_method_; } @@ -929,6 +965,38 @@ namespace ReSolve return 0; } + /** + * @brief Sets Gram-Schmidt orthogonalization variant. + * + * Records the variant in `gs_method_` so that it survives re-creation of + * the Krylov solver, and applies it to the existing `GramSchmidt` object + * or creates one if none exists yet. An unrecognized string ID falls back + * to CGS2 with a warning, and `gs_method_` is set to "cgs2" accordingly. + * + * @param[in] variant - string ID of the Gram-Schmidt variant + * + * @return int 0 on success + */ + int SystemSolver::setGramSchmidtMethod(std::string variant) + { + GramSchmidt::GSVariant gs_variant = gsVariantFromString(variant); + + // Store the canonical name so that the stored ID always matches the + // variant actually in use. + gs_method_ = gsVariantName(gs_variant); + + if (gs_) + { + gs_->setVariant(gs_variant); + } + else + { + gs_.reset(new GramSchmidt(vector_handler_.get(), gs_variant)); + } + + return 0; + } + // // Private methods // @@ -1022,6 +1090,39 @@ namespace ReSolve return 0; } + /** + * @brief Instantiates factorization solver selected by `factorization_method_`. + * + * Shared by `initialize()` and `setFactorizationMethod()` so the list of + * supported backends is maintained in one place. + * + * @pre `factorization_solver_` is null. + * @post `factorization_solver_` points to a new solver, or stays null for + * method "none". + * + * @return int 0 if successful, 1 if method is not recognized + */ + int SystemSolver::createFactorizationSolver() + { + if (factorization_method_ == "none") + { + // do nothing +#ifdef RESOLVE_USE_KLU + } + else if (factorization_method_ == "klu") + { + factorization_solver_.reset(new ReSolve::LinSolverDirectKLU()); +#endif + } + else + { + out::error() << "Factorization method " << factorization_method_ + << " not recognized ...\n"; + return 1; + } + return 0; + } + /** * @brief Instantiates refactorization solver selected by `refactorization_method_`. * @@ -1076,36 +1177,4 @@ namespace ReSolve return 0; } - /** - * @brief Sets Gram-Schmidt orthogonalization variant. - * - * Records the variant in `gs_method_` so that it survives re-creation of - * the Krylov solver, and applies it to the existing `GramSchmidt` object - * or creates one if none exists yet. An unrecognized string ID falls back - * to CGS2 with a warning, and `gs_method_` is set to "cgs2" accordingly. - * - * @param[in] variant - string ID of the Gram-Schmidt variant - * - * @return int 0 on success - */ - int SystemSolver::setGramSchmidtMethod(std::string variant) - { - GramSchmidt::GSVariant gs_variant = gsVariantFromString(variant); - - // Store the canonical name so that the stored ID always matches the - // variant actually in use. - gs_method_ = gsVariantName(gs_variant); - - if (gs_) - { - gs_->setVariant(gs_variant); - } - else - { - gs_.reset(new GramSchmidt(vector_handler_.get(), gs_variant)); - } - - return 0; - } - } // namespace ReSolve diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index 49ababaa9..f660d6ad3 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -67,9 +67,6 @@ namespace ReSolve int solve(vector_type* rhs, vector_type* x); // for direct and iterative int refine(vector_type* rhs, vector_type* x); // for iterative refinement - // we update the matrix once it changed - int updateMatrix(std::string format, int* ia, int* ja, double* a); - LinSolverDirect& getFactorizationSolver(); LinSolverDirect& getRefactorizationSolver(); LinSolverDirect& getPreconditionerSolver(); @@ -85,10 +82,10 @@ namespace ReSolve const std::string getRefactorizationMethod() const; const std::string getSolveMethod() const; const std::string getRefinementMethod() const; - const std::string getOrthogonalizationMethod() const; + const std::string getGramSchmidtMethod() const; // Set solver parameters - void setFactorizationMethod(std::string method); + int setFactorizationMethod(std::string method); int setRefactorizationMethod(std::string method); int setSolveMethod(std::string method); void setRefinementMethod(std::string method, std::string gs = "cgs2"); @@ -99,6 +96,7 @@ namespace ReSolve template void completeSetup(Workspace* workspace); void validateConfiguration(); + int createFactorizationSolver(); int createRefactorizationSolver(); int createIterativeSolver(const std::string& method); From 543b193a101fc1717a9562299f94ee4be2424425 Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 18:44:20 -0400 Subject: [PATCH 08/16] Fix SystemSolver IR when refactorization is with KLU. --- examples/sysRefactor.cpp | 6 ------ resolve/SystemSolver.cpp | 23 ++++++++++---------- tests/functionality/CMakeLists.txt | 28 +++++++++++++++++++------ tests/functionality/testSysRefactor.cpp | 3 +-- 4 files changed, 34 insertions(+), 26 deletions(-) diff --git a/examples/sysRefactor.cpp b/examples/sysRefactor.cpp index f306275cf..3ef2c4d44 100644 --- a/examples/sysRefactor.cpp +++ b/examples/sysRefactor.cpp @@ -243,12 +243,6 @@ int sysRefactor(int argc, char* argv[]) "none", // preconditioner (always 'none' here) "none"); // iterative refinement - // Disable iterative refinement temporarily for CPU backend - if (hw_backend == "CPU") - { - is_iterative_refinement = false; - } - if (is_iterative_refinement) { solver.setRefinementMethod("fgmres", "cgs2"); diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 8cc89c7ed..bacf6ffc6 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -467,7 +467,12 @@ namespace ReSolve { status += iterative_solver_->setup(A_); - preconditioner_.reset(new PreconditionerLU(refactorization_solver_.get())); + // The refinement preconditioner applies the LU factors. With "klu" + // refactorization there is no separate refactorization solver and the + // factorization solver provides the factors. + LinSolverDirect* lu_solver = refactorization_solver_ ? refactorization_solver_.get() + : factorization_solver_.get(); + preconditioner_.reset(new PreconditionerLU(lu_solver)); status += iterative_solver_->setPreconditioner(preconditioner_.get()); } return status; @@ -519,9 +524,12 @@ namespace ReSolve } } - if (ir_method_ == "fgmres") + // Iterative refinement is applied once the LU solver wrapped by its + // preconditioner is ready: after refactorizationSetup() on the host KLU + // path, or once the device refactorization solver has taken over. + if (ir_method_ == "fgmres" && preconditioner_) { - if (is_solve_on_device_) + if (is_solve_on_device_ || refactorization_method_ == "klu") { status += refine(rhs, x); } @@ -749,8 +757,6 @@ namespace ReSolve * * @param[in] method - string ID for the iterative refinement method * @param[in] gs_method - string ID for the orthogonalization method to be used - * - * @todo Iterative refinement temporarily disabled on CPU. Need to fix that. */ void SystemSolver::setRefinementMethod(std::string method, std::string gs_method) { @@ -778,13 +784,6 @@ namespace ReSolve return; } - if (memspace_ == "cpu") - { - out::warning() << "Iterative refinement not supported on CPU. " - << "Turning off ...\n"; - return; - } - gs_method_ = gs_method; if (method == "fgmres" && createIterativeSolver(method) == 0) diff --git a/tests/functionality/CMakeLists.txt b/tests/functionality/CMakeLists.txt index 3770dd3c3..a90c52bc9 100644 --- a/tests/functionality/CMakeLists.txt +++ b/tests/functionality/CMakeLists.txt @@ -21,14 +21,16 @@ if(RESOLVE_USE_KLU) # Build KLU+KLU test add_executable(klu_klu_test.exe testKlu.cpp) target_link_libraries(klu_klu_test.exe PRIVATE ReSolve) + + # System solver test with refactorization (KLU on CPU, cusolver/rocsolver Rf + # on GPU) and iterative refinement + add_executable(sys_refactor_test.exe testSysRefactor.cpp) + target_link_libraries(sys_refactor_test.exe PRIVATE ReSolve) + if(RESOLVE_USE_GPU) # Build KLU+Rf tests add_executable(klu_rf_test.exe testRefactor.cpp) target_link_libraries(klu_rf_test.exe PRIVATE ReSolve) - - # System solver test with cusolver rf and iterative refinement - add_executable(sys_refactor_test.exe testSysRefactor.cpp) - target_link_libraries(sys_refactor_test.exe PRIVATE ReSolve) endif(RESOLVE_USE_GPU) if(RESOLVE_USE_CUDA) @@ -56,9 +58,9 @@ list(APPEND installable_tests rand_gmres_test.exe) list(APPEND installable_tests sys_rand_gmres_test.exe) if(RESOLVE_USE_KLU) - list(APPEND installable_tests klu_klu_test.exe) + list(APPEND installable_tests klu_klu_test.exe sys_refactor_test.exe) if(RESOLVE_USE_GPU) - list(APPEND installable_tests klu_rf_test.exe sys_refactor_test.exe) + list(APPEND installable_tests klu_rf_test.exe) endif(RESOLVE_USE_GPU) if(RESOLVE_USE_CUDA) @@ -243,6 +245,20 @@ add_test(NAME sys_gmres_left_mgspm_test add_test(NAME rand_gmres_test COMMAND $) if(RESOLVE_USE_KLU) + # System solver with KLU refactorization and iterative refinement on CPU + add_test( + NAME sys_refactor_cpu_test + COMMAND + $ "${test_data_dir}" "-d" + "${test_data_dir}" "-m" "${sym_matrix_path}" "-r" "${sym_rhs_path}" + ) + add_test( + NAME sys_refactor_cpu_asym_test + COMMAND + $ "${test_data_dir}" "-d" + "${test_data_dir}" "-m" "${asym_matrix_path}" "-r" "${asym_rhs_path}" + ) + # Using KLU as refactorization solver add_test(NAME klu_klu_test COMMAND $ "-d" "${test_data_dir}" "-m" diff --git a/tests/functionality/testSysRefactor.cpp b/tests/functionality/testSysRefactor.cpp index 86d629057..4f927ce72 100644 --- a/tests/functionality/testSysRefactor.cpp +++ b/tests/functionality/testSysRefactor.cpp @@ -42,8 +42,7 @@ int main(int argc, char* argv[]) { int error_sum = 0; - // Refactorization on CPU not currently supported in SystemSolver class - // error_sum += runTest(argc, argv, "cpu"); + error_sum += runTest(argc, argv, "cpu"); #ifdef RESOLVE_USE_CUDA error_sum += runTest(argc, argv, "cuda"); From e8a227d4ba1f1329311101cce93f980ec33a3442 Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 22:51:38 +0000 Subject: [PATCH 09/16] Apply pre-commmit fixes --- tests/functionality/CMakeLists.txt | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tests/functionality/CMakeLists.txt b/tests/functionality/CMakeLists.txt index a90c52bc9..38e664ca0 100644 --- a/tests/functionality/CMakeLists.txt +++ b/tests/functionality/CMakeLists.txt @@ -248,9 +248,8 @@ if(RESOLVE_USE_KLU) # System solver with KLU refactorization and iterative refinement on CPU add_test( NAME sys_refactor_cpu_test - COMMAND - $ "${test_data_dir}" "-d" - "${test_data_dir}" "-m" "${sym_matrix_path}" "-r" "${sym_rhs_path}" + COMMAND $ "${test_data_dir}" "-d" + "${test_data_dir}" "-m" "${sym_matrix_path}" "-r" "${sym_rhs_path}" ) add_test( NAME sys_refactor_cpu_asym_test From 372c743473b74570ebef35d411f7b2f9339ea157 Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 18:46:06 -0400 Subject: [PATCH 10/16] Update version to 0.99.5 --- CMakeLists.txt | 2 +- tests/functionality/testVersion.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 89fbcc29a..3471e4252 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,7 +9,7 @@ cmake_minimum_required(VERSION 3.22) # Adds version settings and set variable CMAKE_PROJECT_VERSION -project(ReSolve VERSION "0.99.4") +project(ReSolve VERSION "0.99.5") set(CMAKE_CXX_STANDARD 11) diff --git a/tests/functionality/testVersion.cpp b/tests/functionality/testVersion.cpp index 2e49a066d..571db7229 100644 --- a/tests/functionality/testVersion.cpp +++ b/tests/functionality/testVersion.cpp @@ -19,7 +19,7 @@ int main() { using namespace ReSolve::colors; - std::string answer("0.99.4"); + std::string answer("0.99.5"); std::string versionstr; ReSolve::VersionGetVersionStr(versionstr); std::cout << "ReSolveVersionGetVersionStr Test: " << versionstr << std::endl From cc348ca1e89508c9baf2ce5e189355dc18c5ef05 Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 19:16:10 -0400 Subject: [PATCH 11/16] Get string IDs by const reference. --- docs/images/ReSolve-Icon-Red-GitHub.png | Bin 0 -> 26185 bytes resolve/SystemSolver.cpp | 10 +++++----- resolve/SystemSolver.hpp | 10 +++++----- 3 files changed, 10 insertions(+), 10 deletions(-) create mode 100644 docs/images/ReSolve-Icon-Red-GitHub.png diff --git a/docs/images/ReSolve-Icon-Red-GitHub.png b/docs/images/ReSolve-Icon-Red-GitHub.png new file mode 100644 index 0000000000000000000000000000000000000000..6218d36466540e945048dbef3895d54aa1ad2d84 GIT binary patch literal 26185 zcmY(q1zcP)(=Uv>Lvfel?(XhTptw6rad&rjcZcF$+})wLySsC@&->i>e&6ozY;v-b z$;nAFnfcEoTuDI^5e^p)1Ox<8T1rd>1O$`~c+tT?0ZZf|6fJ=XsH2Lc2uRHo{t56V z$n>YQnY=s*4KNP_0tt!+0{)K#_yqyQ1%doe9t1=N6z_lYDxg&Vr2zv03AF@)_%DqX zF#VUM`}gy|DR>^}|Kpek_J65C+48{uH~+5&NawTkzZBL^O4AVpgn{he1?pBT=mKP= zY5C))(@%Lh9^+pC21AoyMy3pI0K0#+K=|BvfLVa4lOeGiz{b{*$Bm!#KNLK`{J&yG zQsV!RI9cA^tj4UklKni+CcUvbzH+ow~vj1W7fBA@+IvP7z+BsSNvL*h< z*U;#fvlBll>A!~l_wzr`>11j4|5~zj{BO5_9c29X4ds{(L*jQ{)E1mJM=cA!B(gg~Ulgnzh! zUU4IY@N`J~qIVfN5g^vWVov=yxUJQqE8O7Zab7C=iigaVPL=(g^AS&m`P2kAfi-i)W_E1 zECnu!4|E!r-9m}Kh?dIeZ)dsPyLr*LZ196x8jl!~0rUOTfNf#p0<^snwWzw^ck9jGO?FTH+?% z)vs(*^nz~fel=T?Vtdtm+{>Ij`PIjahgKhs?CdfCAw^<|mqoc~U>lSmY^98wKj1 zokF$@=N9H7e%?4W&SKr}HY`3~s>zI2I6}mG!GV1vCJs^uwdnsw!W^IXnSF^%T*X#~ zP=OQu^>KYUq^9GxQ6e|a3!W}7V$ttU`vmDvf^y8+6EZ^Cf<`=31ZKp;vGKkGRjt2Z zY@QJaQ4}F_kYDV7GQllDFHNtO&Hs*{q{;vT79d0id@`e648i%_oswAuBaa`GR!)sP znXBog7_F7vE|^kmmwIsoB}HHL%ncXhdfJ7LglXwe>HEKgg;2u6{ZO>z@^;Ov5&nEV zm8ue1D!sbD3I8-GW8t!1#$IF%=FThThe6{LuBMK|*ju7}fdCcW1_3jLTJ=Y`r56CZ ze?{*6GM}7Tjo7? zBrC+iK!_a`+o4vip&*>D5;4#xUrW2zo32Cv!N_W}ln+_# z+cM@5zcM+Kh2`OHXekHL)bA$^DJgQ|IqFsR5w(G4=0klF>enQSB209)QH~nmgjNIn zA8##w*iJa3z?gMuW~H%#wT2$hh+)Y^yM!i_C2vj_&Y z^;#X-6(g_528ffo?U!Ezv<<>lDiv9Wb=YZ_?iX#B;iOT=K~q(rK0A%POAZK<$YQ$1 zgjK3ZflI5HvcJ`AA%N)_v9U1unx57GS_o=@ttus?&GIJBpyc^t$!Pc&jim}+N$9;l zJbRsD!yIuzsK@xwk2~7kAze?`vTVQz#~-A-@fLp5&dWg5QN^9B_T3Qsb zQe#n|U?LW35biarMx)j~HFtuVF?2OXlMR-cOisB^$y^r4)zwrs^8%f6 zpaAMKBwuMTQ0i@+$@8jY5t@3~n;SQ7jyI$`|9dM9;0}KFwt-^5oDw?VX;N}04{dog zkbNE9UqA4Jtv=M8Wcv=S#v?@4(A{-cfOwSC3}!@D)3huf@<94Kd?Z0P6tov|&=4+N z8$uP_3tY5JsVk0eOXUQhbDG*%Wld8Rd!vwdF9$l>hNC0{1)z%IDv!6-9i*YB_5BhT z?&lV}RcIQkZ-P7Xx5pKa_1c*3P>Un+pAU`AY&t0ynh??|@+80^JQxNAyXI9(bp9no z*O=@Xy@hxY!iBuho+n4m$pty)wU^e{+#LtQ3_!$pGH#9xtz~K-?vI#iJtP}6aFG>h zkdtBFXg##}puv8ueiIk-1FxukDJ^o%g4p%w6*xE|A;DuQ47>mnuSj12;VfGUjBC>; z05=&`6gJR7K@q`&T(PVyD4Ex?_(F&J+yJ#ioSTc7nOf%`@i9xnUs_6g+Ut*w0Yq$6t{-ZK@xo@IZ;5*GqCR*a4W_^nI|N(*pj;7q$|)Wdp#|3 z{)!2h2w%^NnZJ9lY{+2Eq@AFICphdqWY+$^hLn7lLl)fDya|c;A(-I*Fb<3O_M4CD z&z?61Q84%J)8dZAvxru%dzr)G6*TCs1CJHm-A{jun&|q<$zW}@ai_<`CvsS@rT{;b zRVsd1G(^}?6G)y`vV}R)f}1v!T7CD}ht~4U6$mr~w17ObAV>IsK`2H20M#O-;e~`$ zSw7W^zU|~feq@{}A1MNXWClWVWz6jo552(N^8|){KR;oXxLE4>=sR5`wQ}y4C~nTy zgsu0-D0r|yy>)xbL5OOnmJC8UqWR$8Jt@tJTDVe?O+ zAqx2MlB8Ox>38@eeCf@saZ|(>RT|$p#>fhwcWbyAM~m* zNV>bm)L)CuR1`%eIKnNSiNDS}wYARWib-?OOyizP>bLK>juMw@my|zj z@MKxZBiwWa6<(bN3kj1k^HuC{DWi@kPm7xF&$V9sv)Duvf=1gV%%xO>*EsEJQ%!c0 zHhBiKN|-8kUe&@fwO-fX#VjRwttig)S?m@}7Aqe|;VvSF|HHdUi}Y&ldPX`%S%P>8 zo^K40cD*YM(wLntSV>kYH*Ufas92Msf83d}ila-F0McWDwVv80u=b8iB`tX9qug#P zG*lb3R7xoEIAYMT2x+l0n+B$jtH3$CaI2-4d;4z1#JQ+6=`quS)q=@Gs0F#UJ%{b7 z0X-6*8z8iZg3Yq&y=axbA0TD6VAbf+XEM z;)2Zx4HrRxOI0;dSE_9Jz|3Z2g}l(rlwV$#I0Pq#inol~KbaM!f;_3>*lXYGv$Uc@ zUP!DY4q=j^B+e9y23b-YZ;T@k9`a1NM%Sy^v8*)sx>1N=CEQuo)Q!6nf;Old8RS3E zYMEWJ!kM4{0bxJy@k?Ca$jb-{0I=+OoK#H9C5bRfN~8NZLvSwoq^9 zZ!%Mf`7q8j5o<9vY7^%9NY<*=iQ_5UCO>N1f?L5w#4$pbJIKavJ~96W=b>O+m~n|Y z`F+B;LFr?qr1c5{3#*oEdt*7{iZ9)MA>mS(prYDA5v2?QbasYD7yx0~gZVfUFC#wV zBhk&Sz3EybzTU)lVe!m=z6-Vz5oB7DpxV?P1OF|ueT}HM*b|q2dg!~Rjh_;A5d2Zn zOM?ykZkaS1#1H{!XKKPJqW@;A^!#ANh+RUm~vmb^$8;en&P>l_U9|h|v^wgXX{J8<6!kwK0~IB~Hfk)koI^;-)``Ia3Pj9nKjC19Nq2uTFiUt#))0 z7*AMPeQKSA;%yg+l3j3Md;jEDYmdOKOq&Rth{GXfivUc2S9Awd0qt)?;YyJKM6CG% z#KELjOl`mUJQd|uQVW}9KWx_$;i%ScHh*zvw>Q+0fig=G&lkIuY3vOr~Pac*8}6Ejhn`x-5w@HcSeb)$~`^;62u zoh{6-k0)2MK_sfP0+Z!eC$R;nIu*xpxME8s09&9$YmBR`NJ4tvm#RjD(8X*$PXJHm z>dSrplyE8LG1}aK9yg$##gTKYk$>4`zr2ln+ax4QCJ9*tHCJbgD{lO`{jFqul4e}& zsh0dF-FXMP*6Ja(*4{BsGBV^}PxWSI4f@rsm)P2qM~`P+0SJHlf74-|sCNVp1vN?G zn}(3ACtQg4EsEx3V6fmIKtbXxd10*8xABJ)@_unj3067f6UklbcAtHoAq;O&^f4cm zsi<*pOHo+o&vOK+J0Xde?=Vq}KmorO>8oZ%DQu66GoZ^#i<6FNVhqkWXNal)Y3|dt zu{sqJerc$Aeg*j@P!&=`67GM{PD0oJQp$X|`u8+G8_;W2QDvika1tJ8^-r&qXB$>t z((jYsQ4nyZ->Yfk@`FTPJ_&FMNfTbKw?wYcj5es!Tl*P)s6`*W6JAT`A^*!SKJSMtDJ_jC5vz^PRy#lIaFrr^wiHfd_Fz)+U*yA$wCs;p-E`xL zSEsu#v1KaWlB`V|8?3hlPvNeAwj+QHB_38g{OVi0mOf5?)<9A7n3WQFVM?U1Sa(s@ zF>d#;H0n!SCx=;vZzsb>ebhl2=n(Z5r=oW7&--MOEl5ydgD8|^SrKd0J)(?R`>Iv) z1u5d;1mftovj+NK5o{G*w~0=~1JRj3x8A2Gbih^CI!HDgz-SAazmbU^sXpBwx%Mh%Tidh!~%^`@)G2Pv=js?zm>eFB6+e z6s=&z54dr2;sU29>gyswB!x0}G}NC{8<6H+2vmX$^l(jJX?~otzFSlGy=-H-v0mbA1S~3V+cSU)OV*C+cI2@WR`s=AQ~r&l zP)kEw=D?k6uagjhsRkR6rUT3{eU@X&_uD(M zCRd)e%H6Hi&A%-!kr4%1SRy!}bWimAe;)=&ae9ZD0&$7(Y_0J{ZrtVYagnAa+iDt< z59GXyyL|iAnXPm5IzytHxj)XLR*x#XD&fD@EE;u96cB)$0*z-3W(&c|2horriK9o9 zVHyE1MWad3aV>B1k#YOO&`7Ny)1`n0C^x}@KM+hY**h1MiG#80dYbJO49Nx0S!uU!_# zx~{rQlJr4=E`NBezIS|783HVfWjeWZ=|T7o*4N$x=p#3&Gxtk=W2Sfi3=(WE7pKW`2sx%|D%71G%tT zx>KL_L`-a3NPA-SJge=QFO-39W{wHn{OFZJM_mA5Ew1k!;TmvE_|e$v1GyhcjDiJh zzEz5Zh@D;mos(mVXAIjLl$6Aw5s8!}Vqa)z#8i%@=E6Hfp}Ey;6^wqe-4J6QC{F7I zh${USOzJ1JNiE|7dy|}8XcBp1pE5;=extb$AE2#X&<%o)RPg~Oa9EQ)#afv8kLQIh zCYq7C$_NcY(bx4G1?C-C7(W#z(dZ`&$b6$AfR6!jOs^72=bmnKSRdHL3Z9WT+^ZuK zL_vvQAr?6)n&?@R`?xvo^-CoZP^4j{t)$P;L2b&-;y$Q4m)VKGN zy$aoMUdXzl;zXM-RswmZIS=L9bGG4(OOoXV5PY%av9{A?Es4C!S2qbR7^pO*7~sH| zQK-1FjgvB+G>?sh3*QgV*SYgjkc}Y|un{7OpczZ`iGPo`=_le6!k5gY&*zOc_Tg}b z_UbF}8moV|vV?1HB#RKg>f2TUyhndpD&+pf9p=QESh$SpWZrnzRN;-(nsfldsju#} z$Ubq)yw-AeQL8k$4}m?k;PpQ3Eba8V_z{d5yQTTr@`0r|rO^jTahpTLqHvT9J9c{J zpHWVyWSz6jO(&WXnedkQ(6h;lCLX~M4;UA188Ux-Rr_b@Q*?|kzg>`-vvp}9aG1yF zg>_|FdS7E}J3Q~tFiz&^<~?vn9Scv+|3LW#c&;EHHW3zrp&v>_I@h;`yIvKXVt)Oz zKv(61SxDDlFhGXE;fw}o`MAfpKbcm#oNZ*?+3ZQ4bWqVnDHP87>i3#%aI3Q zA)FNp+oyhVIFB1KKoXdFDAMtD+t>kRSiyD98a=$Ax#ZFIDkBAI=~S^m=*`{cJevc7 z2PHK7dmpSZcQYA#UyisyL&hl6;03_Mh1v8TFDa*qdHxo0^sGfH8d92p?w4uRI~bs#Q3!OQ068g%n4ET}Qdg*;>NR8Nl8unaX`r~rz$ zyz%4ok3nq5c8xIqG~qs;1{AjFye8I3=-pj5haPTeHuLzJjgE8?CPor-u4HZ(O0r@> zISHeN*F`$y^}8=Leqh@m_2vCvb!mi~;bJ1f}GoOIi$N%&|ub2%S9@&py(bG}~~ODTTvptk~xJ6*86@5`7T&4l^jm zQ~3duR(pp&ecNLi5|-6z@$eW%1u^NqOQY2O^LMzRt}|NtC_f_$16`@QH83+Fv7*7}ZA3Pe50c?FTDZy03wD_Hq@jb7ZG+3FS?(H5lhw z+?OffClktaloDI&rZ-$?uMukoO+@^A*wVxBtYs?%P3ZlUhg8+R@atZ>*g{Zr`I7H@ ziz)3|3&sU)vStk|m^~2niVQ@%g_t!@Ejq+%slHA6?8Ugo3E7RyfVgzj(&QUK(rQYe z%#@req#o9->v!v5)#41uYW4RlGZ7L%0+)In`jdgftESoxNP8WJUHIFVPF!G&LlJ~` zhT{s*Ejb}S)MrOv&abwxsx?Rd^97+#tkZK~0&ChTRZz=y8pP_^c7;r z4fYp6UJMHC-pXqRCunBIPn!eTs*-={^yRH=BVF*iw!E9Z*nFjGY+b`GE8Zbwcuee% zl=`mox6@O&#v47`E76Jt0bz7;)fsWJVV0vq8Hdw0x zOrkhRxW_|=0CUgiZqd>?7h%JWjBk!=Y!*nT&ZJ;Y-XA13NrJc;vprm+ zY^4E(sR;p77>wCs70ojZ9+?U1K@Yng8xR%_i)a#-W?`1CEwv$c*`4}mEPH4zu@8(qjO^Y$SO^8ZnaY$Y+Y2V}QBU1J4C7u`w3(Pb#4y{W zKNu|PD)WVS8$G43sS6%)lkFPP2GG*%^2>wnCnQ4Oy_iGDmyS^`Gnrw}YEzVt)3dLV z6LNim?hV_aiUiosLc;u+l+XTU9WU>JjYZU5>jy1$qB7?AjLRY37tNAzbmEtk*uzLE zF%#N9bx`Pc^cfg-=0BiDc%}ibR9l|zXE`}Z@W6nahVjg$un#mz>l@lP=;q!ZGiV`+ zj+dL!I3zX>*UcBYHzp_p7y)9R{iEhJhO9U>8M*f8t6*q|FeoWs^Qf?U=_`U-EbY*5 z;PtFCb1`D;O?R4-s@NZ-4o_}B?S+N2B~?UxWJm1R`n;%)X0H%%cl+<)vt4g+@dHIE z^*=gNhWFcnDx5wesaTTG?)Vo532L12_m(Tcv-s6JJ}LdnEAWSFW8 zQfy*ir^9rXdtj@zM6RSKo`y=fFn!}^_lkI~vU$F|yh|bfNWDpPzfKV)!F|KQK%n>0 z5tI0SGcgLDe}vxH${eNL(7}$7&Bf*?wvCKF;uFhZUN-y8&r8dOUAKn;pD1=%6f4@rx2c0}xWUA8_BtrO%2!k_w`yZ zh>r^a3Q89pBsU(bC3$@&bsc`=Ik}$9B~?|S4-~z&?rH9-?I&$0ZagFL4xww-qgH$s z>1F3^_&$8Os%8;mZ>Mxw<=pVA=^MW3Mf7oygD9a)q7{2Q&6YYprAb9^;r_c1sJ!rD zSIQ7n^5lTvc*mRD_X+%fD6P&r0|zK72d&3o1ir*?uyY{9g7*#nr*sViS5XGiAf_eiD zX}HK^PHEJWBFY-k&fE=|%yX%S8C@opbnH%Bd9`pdY$IFYZ|~;lhF_zjEYh9o1lO)> zX|_Hj_nlA73OM^nSckl5zh^CSNxThQJ!v1Un@QdFq(iBnG-sx`7V)jda7T|n#QNKA z*gHZ2*B`&HCA1-nAh=VCq!*Ht%>N{O7XW)CnhD>UoINw-nxWfa>cZ}FGQN9xvb>Hd z%oeBDhBVTFy9h0Eu&z-Z<3!#91`zg%dNRdYD|l|OeV44T zn4|{#)bI$G#E6dR54ZIM8#x@tBOy`%@~*}Vs|+Y{T^s2x>oqn;=n`L#Hg?px4e5pRkHuadQuHCl&((2=E&1UvtYG1 z63>?i?k}T1H|jsdb+U8l&&b+xR{G(s{>1cn7-NBKel*dnP@l#xkQq#$#{ci|UL`h=t zDaFZbGzEf`dyX^L5PKXkref9JRcI?et%Do`g=mfZzsQ&G*R0*7S)Z$6gfbgO{W!w3 zVAR`Rb&AC7p@!YejNK-SKwt1meg6XKn{CP6s{_}|ln}R7st)mX{^YutZxB+8dGwS2 zxN&UJ!cQcin0FYF^E(674vO0^JGdXY*!ZrF)JhmM+Rx5`SpHgX>SsMI401DAXV)v# zjmlU~0ovT#oAZDXvP~E?S9pxN$(%|=BAv|i4j6LI1n^X{wp(scpz@gsf^1*LAZEWD zqoMs39@f}$>n|!=fH&N_GA|^+ahYL#FJ8P`^Kc=%`NSZ+Jsy~stz2lLoOw+lPPuim znw`#)likIFfXXuZhVRcNEqDN#wv9{}xYeU>3b6EHc#_A5c}L+$2f8ixW}V|=fUi@W zH?$-DK3Eu!$8&>`hv{3BOWpnX(cX5-JeRmQ(YaHDqNlX206XW%8_E#vgIs~$984P4nW_Q(cEYsTpW_wR01(>-Q zk<%^)<4=s7pK1&W+DYDx0_^OJSM-{LnUdZHqE=JmY5C&xDXtzUDv3t{eDBh?RVR64t=3V$pa)^6+ZS8d4aJ9)BIV9&j z9n~R0KQ#NAVY{;Gjsm|AwD`lnQz)3wTCYGGXSXj&ECnOz>qqcs;;c-$xsBC&i~Yq2 z)L`*5nVofsr3{o5Y3^JLCzAfa!oJn-e!bcpI0J>Mh!*J0dH5h~DWwXO+5z<^0onI8 zYPwkrs~=Yzkbwn}Nnm%n(Uo}>$mbr*1!d_IX5}oxzeCURy=Dy}Dt$~LPb<=egycTU zdemxQF%~8f1s`PGD1j`SFSttLjMU+Jx#qdvj01e#C?(tsXs7WgkQYpL>Y&;q0}#0TK6&mxW|P?*L|s2KkN60wH7=9RsbT=!Kh|yk{A+R zENiYP{k%plD6uV`6l$*Y+cyzXqK7cx&}Z{EHI*M~#wt(kZkh~llZilK2+MU z?Ei7jXeFDMjpIG)RK;Y(94g!c@LI-s+rcnM~9Y_hFocU5al zm_K>t;oNW(<$>C&(Bi7LwN_6W9ULfGdE8<&f&lADfKQFAvQM*>2&Sc`{@EB+-)X~p zbp)Mv1RI;;_8>3@vd-I_%4~Sy5a|B4l-1SR8dUECVWdwDoO9R%l5~xm&CI?(1$txa zU83KJBJo_Kee80WzxV~W{BfF8bgH&{hhpq)N&oI1?Dq+&sK^SRk?%|~mhg8_Se3O} zS82=*0AaA|p|!S*Q3g0g9b!IOo{Tc&B9K9ebY!7k!*|C!PRM9XU&M!n4Bf}5wQ zQeA0#q&4Zky@uODIC*O?EbBb@hOh46LN_;!CvdH0;hP0+91`~Dn6Wap24e%sxhW*I zWp)$9%J!?hA@)g@KUwR88kkw%E`S(}+DueLUY^{D$aiR!o*CVxf2iR3*ex1FaPRpX zeC|`U@l;?mHE;lJ-=i*Y4abc#F%*7Tu2^W4zLfx<)K zuvM9?vo}hzEAfJ3F-$alCBYNQS0uXgqh&&YlHiXZMwRr6gPSN;6wf~Ni~fi&;ZLtZ zT%?cZ^4les60vCdh|4z#q7246+QVObj}~0;c_xbtL+7`X<8LqB&O+D0wPEc`8W%Wn zPOB%K_B;=<&6YD3&NdL$1^PjkI?%)#(^G}=nyF=~ORC2#)iQGeRGb}Gx>W3!-NNki zvR_VVlWzDdV&i^AzDql=LPA>jx=O4yq)}QLX_p0y;PK~P(v#n-!ImydpC!HT0YS#q zXTiXF`W}!Xmx9eq+gFhHO1&U8aiqs6f)_(-x8`;oz(BdwQ8mDZo2S#}C3Gu(8L?Dn zmDvNwMiv*Z-7}7)`68;1m@<}ZAJ-*qcivLtH-hOCVT+_l#BPl<)&l1f&6Th0psuJu zuB!CN8v_j(7y|qT$x5b zFC22v@xA0(RwB4woh3^3OPrG@w!Y*^Td^iWSsa(EcR!98EF4=BqLhk$zd=50wAoen zV5s?~umDk|)&;>oVt^}QX>;MO56|PZ-?b=x2p^}X`UUrlN%cT~_Agg;-rbsbu0$S^ z)D$`m>wgXjTyK`yF?>a-0Qi0|>^9jNsa55eFmxT~pcI$hcAIhLbF^{@-}5ecA5xev z$ayJ*l@UZwZyR}vq6B<(20Rxk-%OcrzhDnko2$E!gZ^~i1=n~QtbEZovMIdrjrGb; zU&1k+;t^MdmS-z_o~v&96t9mQLRR0(ppn(w=w2iX*2~8VwTU*a<$0qc>NtcgE^M2;PL*qQdr!;Vw* zQ9ywhLrQ?1>ux`nD9ScrRb>d0xU;dKa*CvStDa@_N&f*Rc0>>d$XhrtVhO zaW74K9)y;ZnmccC5Z~y{{FK*@cY9f?1f@lALJd0&W&2Ogy7>G347=<%To}XYMz!By z27KDAyJ^u2gX80~pwEpUd>2;pjSvK)6p&{4jKX`*uYBl4`X`RSYs_dq+`^mD-|DtfC`FrTC zIu=FdAe{voBg_M4@(eWVaqElmYowM^pd!Dd&*mL5p;P5D)dHE$zfrhTBGscYzABD}fKs#hql}{p8Z>$kn2QCE(idvb+@4po? z%6PlP-qxo7DcAsKcz@tKj}dz21Yjw&S_kZV33Yk_gLQ!cBr~e}_i?oLnJf)heuF8& ztLbX=MjVMcdytj!p9g32ljbA|XZUBA{?sS6bPJ_uFa_uU7r|EeB##T5hYVnT3AJeoad`Wm{!i*3CY@$oPvOPpM=wqYensC++ zR?%pk^db_%OVTyD)%SqtwDKHbP-JulvYpbM$K z=>g{8dcH@7N`m1Wwn1M$b&Bg@+Qsi}>e~P5{oS#O>IEAD?&MRQ;vF@FoE(Q0KtIZS zFUj)3*h2lwtf?s|FxN;ZgaZu<6Wxp#@(8i>u`xq>l8R&4Y?1ElDrJ=Gx0s>;5%;HF zyb#n`G}^g6rDM}Kxs);jW%ne~pkVKdSN2;rq~3KHJ-9At3d1g70FLW$gx=H4}@T2^WYJa&|qFYA9^NmQ`rw-)* z+Q$3x9u1U zx&;Q%u-nEjL#WhF$`)iP+uef=#1=vS8z1G}OFW2?9oHN7dm|3uLv~etnJL`0%nnwk zSJ8fHn9ZJJ*bF|rQ2oWWs?%?2J|Wi4Iz-*?BvjA ze%pFd-{W38>+v^)dRhQ&`Py+b4i^mzxhxyU>vc#6#uXkq^SHyHJ0l&Q|1lLh`3F1qiWHTEVzqfuk=REFDfZ9VP9y%g8*|^NBd`s=BOFYy4&Gb#7T*H5b>H??z*eRF3Bi@D^#`17j`;>TdE-OBb%W#K?mXTU z0PNGYx!TUm8s4aI@E7lzGL5)5>5byz6dL2+pX8Sm5Uh2wq=ZxKe+$LS4O`7;oAkyF z3;R^I+#oBNf)qgvIpn)o-*r3RuqLWYxbd>U))mmJBt(_1?{{5hPq8@1HG(L5WAN%R zl`_OVUGjcl3MIHUxhQHO09^=rWN*y<&TwlX1y^B=TGV z-q?u?RkQ3+3ZnA&y%UF1I1GH35e0t*m+=sjRmT`oEuN@_AmR5ba7C#R%~wc8b>0v$ zIQamGVnBP+6WMf-Z)aM_Py+7^EX-@0bvvk9oZ<2O&YX&hkyH6UU}%Gu|K2g8o?tFR zjZ!aL>=)NkU9L93=lHtVjgDPPf{d)at;u334qf_T{lr46N>f3yg<_xpdVtXtWR7CM z(1;Hnd6U|KeX<~{i?P%95)&66`Q_JZ|o%-zsOKY<$#oq5}z(>!7C zz5o$XJhtyuHeK-E+A8IY1P7_I7BC-zHZbUgACdpU((M3B+GEM7XHDTfd05>XncUpk z!JuYGnziBH@IIE^(5MJ5-9lL`A^RZyI-WL9OOm2R9?;(}w~OwD;90u5AZ_slz^ZPD zr&o$34C9RVY9?nc5iYo05;@mGJ)djk@_f;9=x36rvO&1zs@V~G|24FfB0c+h~zV)*9qUEO@T`31EePZO%k*y+!;`)=Y@ znxYaEmax*&VdOf0E#S;aY{c8pEe*aQS^+F)BS^%$#*Gj`$7Ge8M=Caa!Jlz0?FUUV zxk!&g|GUm`EYraoeERDjx&U?o!M*Ho%L~37DRIOR4j(_0M#%?1 zvD%dzkLQ3tkli2iKUB@WVIm!08e%~RWjN7=C=y^hp->~maTVlRLKN&vt4kUW-UTS3 zWm1(R{qsQGgoebMuJ(am@cg;NLL_ctcugUJTGbNNc;@L4j&(g`Tr55wxQ84BAd`-J z`}=0@+%cwiCn02_USiF_^cuX1W}>=AzZ5h;9G8ah9Z6)K&ZZ;``0*kA96K6S3|loT z{T=0FAIFnGURj%G$HCTmr0t5#;2h;3f=Imzou8pAhKT-c$mN~jpd|b`=$~H^!iKM& z5t%l7dH$1DE|#!1pU zk9KCN3X1#mimKKTm9e_V5$9cQme8|BhaTs_gA9Ma!2@PQw30BT@OwC}=x_9jE1z$; zDxgCzTI5Kge*^UBYJBK@)MX z)Pe$tC4!#5`0(&d&Rsic6&jI zH&WazNf;64RVTg*H1^{zMv7Axk|pz-ZE^%XrHvxH0RB=m^;P<;KaIn=F1rB%FOGmTH$vbLu`Ts`YV^s`Nvth?Jw+>uiG-G;qMA^>@dt z;#C%f{6Z6-eN`oiBDjRj*-o6u8pDv4m`dY`D}GMbtFr$s$1e%L@|U=$k=XHDmuYMu z4HVq%Uc=g70HgLf&VSr@RAvyx$<9o(vrVN?gz zTO?d4!8_I)Gz$SKQEs!!4_WCq3@DTgR5rJjYRN)j^Sr55&-pK6p~ZUB5dFT1(WtXtkX1WBQL?&0QW?(y^;HGYsp`ApZWz{ zrYHBGiFLVmgKU1HcBE*AW(x#abfya$m}Q8HtZ`&bV#Bi<>S(UjzO>f6sk%1<7(9@R z=V;AgUnmHnzvLUe5;#SGNRhsu5^o_ux)|L%0ps>sPC$5Si$VRO>HTnCgO7 zYdU?=Q7cQsq*&jBN4n+<$Q`FYA*?`Q_@hQN8ITX0*l+tj_N3$C8JSUe76If$VD`(q;WqL)u`Oy$gs7X!JF?A!azY4}InMXjEfmQhgt zl~o)8r|b8K4)wrR$0s3FRST9c@quJ5+?Lcu1MOmI5+L6Ha7*J>ekY_CI&1%}2rXz| z*n;r5L@41q>i9OJ+?-KdYk2wyPuT8-o;)M{T!#)ESj#r|k!2^uSUPa;3lW-E z;KUi?QK})8;rJRP(DE>iOh2z{5Mq2ABEO}~C`sqGS|1J(jM47$`9n2=S7JS6*f_=W z?M(`k{qdvhoma;mOs!njlizoN&tJTr*R>dh@pXT_6>aiFIduN8(KA7dL)guCk79g$ zefSw~>{@d?Pbm>kF7^ z1qA#HTUvGbb)0Gj`1cl3RFWtI{=ZffW|7$0o9%gQ&{q(17zdKBt6}tZWU0qZp z%>ag0#>(6)-<#G8O^`aRqLC1m4`xFNVh#r*6Kri9pOpP-?$SyZpEjz`|4 z56u75-c>fm5v*(6oe*Sk4H_I0Ah-p04H|552o{3F;vOWpI|SLqT@yTba1HM6cQ~i+ zZ@6FXr>(8pnys1ce!HLdk(L7tHL1o|%(dlh_D*&W&@)2(e<}qrJUlCT)AFyqH+-8%A3Eik#p-5b@eCo^oKqKEeLWfs^gFu6wb~G zr|<0K1)0PSiUK~IW!92C2bQ|fG4v%7g4+x|Qr*(c2!^@+7EIR#-`Bq>_Uh!_Tv(G^ z0C%d+W=r~ZM4nGYdWNdUZ+o!(a0-{$>KUf;Z@}8aX$GMAEJJxaCkLEdCn4N z<0@n-L#L~4onnob>t7BR&bQ){Y@b$vskm*!K;)1h%-j~)?jsc0P39fX_NI7?!}?P0 zc|CeiG@sfc-+=Saprs>uYnyMgAHi;FP5Mt#+0|2}SacC40h^2L;89b%<45$@I^ZTh z!a{oqYrYG5ap_S4Yk_+Cm}f3;Zgs4~12JMKW{*gwxIwI<P%kkmPz{mW^U0c2TPZC~47N}Kdb~YqmcWYs2LdrBc02B+}vVOJRk(vNT3V3yJY=`lU zmTdJ1j`gdwIaR1B8CCEcMr8?+l!UtcrqmFBwdDc{Sc;6N5@TZhMDvT{M74#rR2tFP zTPMz?dseb~Dw-B!$W9|N^lilcP)U_>RM4D`}wnj)I@eya+qJ;iW>7BX325{zFu z0k6UQ>`+a>&)dAD>JTix!KB%1!m?x-N?Qx_(`O$}LM z%(<~QdA6lJ4;t~*3?UoV$-5A9uwt4iD@|f2n!`1BH~8Y|P2RwsX;F>t;x8v}y!9}c z*{C-PqB`xL8vfmT?M`1+RGh93@QD4r7w?h?f$g%6buSsTfi#oyhP{ln)XDXF z#3AlBal=#cmbJN58-EDC$WdsrjO^_ABFS>a6{?-p$)N&`#0kXg8r$3FUsko z6{G|vfXcviwXqg+WQRVZ?~OFz&}5pd_FIOqTw#JzNStXSTHHyc^C=$Jt4_&aOjW^W zLYhBZEqA6AgKI|V9Mun_PUp!;%nU(~$Pu=f#pQ9+AP=S6Q8yLP+fcZu=wI9{bJZW9 z43b2I6z%_kt`!n~)$iIIu8ZH^ilzVJnY|+$JgEH9hMe{M%2ytjg}-u zPF+z{r0^@xu$Nu|Dw6UxEA?7w599sfFJdTB*sY3k7p%U{q5EyAEk@g?ghCtoMl@A_yv5X zIp5kHO!+Q-F2p1*K%Ei(uLt;WX@!QMOoM&}QqutG>e$uRftuo8% zSU&63^0?Hny!7&r`_+ieQ28FSO?aEr)L1$BK<9k{h;6)8BEdT4u{}k+>}Tncy>+pF zog@WSk$bJ=kHmhQt;D53lhP4{NM4IR#n1NrJXbG(Kl3h?s?9F^H~;l4#)Gb2RsqFY=Jme1V(1UV)z80{v%ib%Lkj?mbWHfsJMn2k$OBh;9xj&pd= zfBnY1WH&o_;5g~cgg0|Ui{+5Dn-IfltjKPefJW`x`dw_*h|?|N^RvSh#~{c90Zhb6 zR-RG8C6)`o1?E^cXk~S*Sr1?9g(*oY*=q0)lps6Df+UcL+^Z1IO!OS#ia$pI7x^&4E3|0v}v%=+u`T?QB6gG7E)CZ4AI6D4jv ztlBnljwpqv4$lX^mfza8kD;{eG7PjB_;89r}Ke$NR^Q7>(G2;AVGLEfkV* zd6bPp5o_BD?mcnN&i3&)*--XPFiJtMOKFz0oac>Pg6`u^hl9G(j7h9gimV1ocxhnp zVcFoCHcyF_v89)DUL{e=^el${HKPDT}bt6?66mCA=K2H6_Hay zAwT6oPUnh%HrH8I$ZTIIfH866h|h_6H-ImTgCTqoB=#YyCAsG)kPWFnrllVxp957x zV!!Fw)`zBEyihy_D!}N|P@4HCMh6Mx8&j*ZAu4|vA~q3fWZagXQ~8PTYKUF%`hr^X zS}>USg+htF8k#Y@ADhc)n#)2aXc~bOSBaT#BA|K@`!6nc(D2=Bi4&)+ZjYDp@Xq}l z;(um6FgnJdJ(y>{skT6s<3RD2xTJVd)9ZA(*YnbF;l?qE4!{Fd!5Pk{KucXAA3uWw zQy)}uK-#?7W9zF}FBQe{3j0OZP!h%n|>O0uUAHKj=9XC=J@UdLH>FTWlj0u{)1H!|`MmTVL@EOTe zJW|gwz?)x9sXF6+ni-O-4XaQvNVAt+X_fSmy`{a1we&Dk_~DLu^`N#nU?G1Q2K43& zDTTQ778HTM-tlU>%0>t>p2+1TS)HbeGr`7jJy;87G#?holOq}G>3S9g!`vP+SL!Me z-#Px)&stOD^fic^3yc@WRYalgVtS|@6ZgEP7Y-f@&N zbQUle_+i<@Sq(V4cF)t_zj^V?48p}0d14875OQGFzlLs|o}-tKEjVy9|OPYq+`kyy7l zYh5OQ&4to^z`!Ot8!MJUByF3i#Ph%bL- zxg6ZwGdW_ukk2O6eThh>vqpW`x^tqqEyT=3OGaSus>a=q@f;^3@zy1d$vucch+=sL zv9n~_?CbrsHf)!9(1>H6bggQEK3KaUmcv@`7bZSaJ!)FGa|fK&0y!m1b~#1e^bRGG zwDX|*Cd?1hVm)nv-@jxo`OWvZ!NQ~E(9T+?_QDNE^$y+z&;9M9i((y9X@N3ZzCCFX zN*wMFB-`D_Qc5^Kn)sN1xI?(nH08Zl?;wF5ZwFoB^K>>g9px`Wfj#L_q)wc@$T=f7 zm-tnERY1|%AwH$Du-pgU)XpGE=fRjE4+k`!@AQgaXX@CeJXMb%J7OjY;7vU&fdZWbj+H*!LPP&@wY&4-;7I zcVCSx3_u|e4wmb(w(*IaDot61AjTd+3n9%@uB{tmD?~`7#%n%qot<`2@D*{3bQx_e z+Ip3^{oREqz!b2yObyHcnK*%`9CcXXt4)oN=M9E*e4q7x-#PMG&ycnmp3g9SiIJc5 z@&k!%P2iXMqVI{JWsMlUNuMSi;IprtfTnI&S)xeT*fp_LA9M76M_1TwOF{Ad) zUEfQit(|L@XBZDHd99hSUbf}LBj|ASM;qOOMW+2<SugGz=TxC; zCZvin*^DzV;+J^6(a&s)LE*Pe{Qfx5Z|->4&LH=0LqAk6pX1*Vn#`}hk3AMqIuI%1 zNU8AjZTt_0rbBZ-+rH&bk2bDRxv~XC)l-gnSg_z5r?8sqS4-*w~;|Mx^>~6)t7f zqTW*fMUIcW`OiO>1*UHsF%>-#`F{UGE>P*%vZbJlyP#@(f3&x0{k}Mu9=`(N!=wW7 ztQC}p=jS)Y%k_I}y@H%fQzA}8p-KLO+HboLRkBn};?mws8#&?ib42q0LJ|+sN(X&D ztY0LK6Ff@}BAoCop#JbWAj*W#Q{q7Nd5`4NJ^0TLXv4K;9dgf%*}D}MsNg4n(E@tk zYGcDbrKuypC+LCD5LvTTB!Ze)1ffz)RMlxHUge4-+I>?)dRi%IX>u)lno91Icl}?^ zG1u6oc za90$P5Bs*w;&V(Zp>&7K^<%0h24N0y1oz+e!b;3a!K@K>LI=aNDA3Q3x;s9l-+})d zSb)NAkNd5w3MW$q&_k$oFY<->B%y#tjFT0FnCec9U=SAwL%~X^-ng%rC7?+*JMK0( zsPS+|ojkD!wWeR0{kW`$i`wYqoR2230%8`09y1uaQRG>;lr!R$Y;Ki=Llo#BQvQn; z^2NHTE=(wAk{1!ELktnV@T43=bCVh0^<<783i)@|2*|F77&VXDHSLeqiOx_}60{stLN??7Z`AHE8jDyv&rlT$T_YZ?) zcBG;la87i?NplM=O8)xe5fe`=-v+a6KKHNSa+l7U=6A z8!3`(jTNW(8vjM~?$PQgtcK;f%yRuBm+AHgCkrb76>B1TVk{NncL~*A(F+bzu}0j4 zJYj6f*N_3}EtUOL<7nyIlbFk24{TjO;oHeXV{)4>#QcgF?;~s6rJ%g8WJ$j@ar%v7 zL|#xx5vDQYLZD4#9Nn4t%Sb$fTbsLEYCXigs3za!t?uo|i1aG=NWw!j=gO^hKu!7Z zMP`m4XU7VnA(@zc2bZDdhm5NAM)`M48>U{QX#**EdmW~co<+o8NtG+jn#0@J^W=x~ zNk&x%H;uDk#6Eb5&4B!KRCSYGUV8&k`d*Vw7U*{9bsWg;#eTgc+MO#@{YvH=0e0JZ zO}B>#AL9)RZ;qzj3RG2LXk3WO`%UZ@%sZ?ZjIwGq@NMQMzJ1%Ce0d^X@qtC|Sjwfe zc@~$lwos4EXxd@+5WURQ>|r9 zGr3knE}!Gm!ko*OFCsUF>n{65IuiJcrxGn?)-axyrBc6{n%gZEnk}#o2=_y-dq4Jr zxdv>C+vw%O#AoKZZ?Cmbj83r*%d`*rwck);1v1i;t$OOIBU{_>z=W1AhrQH@OEOD2 z^^xlhGNVHLH`HX~uig(8sP)&}U)12O9&tYrGG^amlz%Fp>8}UHeBvoA0>21=CDq6So@G1n*N=^M0lBDw9fU(U8VI20l7?Gak>0 zqzmuoU<67_QM{z6qT-h7D;0PM7l&5%#jLFMUZPrWcVoiEtci3_S9SY(6?h$cN*#Q5f6t`-wvxCw(Qmi*7qK*G(z=QD#-v?rxjFtg=}ha zP3et#9l3u+;wbt6T~#-gJf+`-izuh<(ak^TsLA0d8cPE-A!r$ML0|h9-DMD@(=pJ9 zv;LMoMBihGi+y~@ME%>m^0^h*y`HfMFa3?qpv-GItD<7@&0|$EX0~|2^wVY}iP_$z z%V8&CRAQ&s9%}+R5{VcervkWt=rQqe&9Y?doxc#|PSI)Wp1j41tXpT+FQ>hEwAs{e zI8I-O>h!YQtAes;C@at)7|ON%W)r`yq0u+zeAP0TMd>$XIk0GZvkh?r0JQ$j5%pV6 z$v(}ju@9KiW(B6v^k5=MkNdbmmIFEFxxq88;>jg-nhdvr{+yvguxI#=Fo;Fod@deQ zi_6e6!k~q`@@GMZUP&g!FParbWZ*F+u5;lYZ5r}XV7!Jxztfk-0CKMANov61Ujl|P78cI>WKN>;3qn4>zKsKG|R*YT!Qet31@*7n;oLB@0l|g z9l&0x?mc;PF#VY+H{2^{qW|%xh@uUF&63XU#8nrMz9(4H|2aZFtpLrWcL8teg}P0X zZvi>o{4y2(hBUjsP4#|m6pyhLo;967-#ejwz^CPt=v#*s?V;;3YJa~G;Ui~ZWszD! zQ{F1e&lu)J8D(KqZZy=(v;y<4+kyEbN)?=fee&+eQSWc~^aHG|`R+{m%k-#`_vgn0 z;zXOGj8N5WWm2hS3Mf+IM8of}7f5lP0}_uL49*Xa&GqWjr84DWK4)>Fg+cgDDJ?80 zx>g0sSHWF`nXBv|u+1Wh>Un8gAId#IzvaNLRtdrB5|(_=i!!i)yxeG()=L`Ui3vG0 zhayZ=1fj96fKKFA^2L3<$Z&|?v1*0gk@i1jYbz=5tbWze7d}-Y%_iZuG+GpJ?ia-9~{s0q1$!7p|u+tu@aknAk zardKzPqzUqh*B|wDn~NasM%t79x=>zw=4V)MuE}SH|#bEFO%)Jv46>ZZgM{|g59{L zjy_$QHivH2%qu)aFJmTvP7m-3wo4OIRagB)%(Vtk>4D7U%Bd$of^ z2@w-HDoc$ud~>g8)hO`!NQCkE*ZZL8tMSpz)kJ+SiH$`Z0J!?X!8CmgoOaM-z0SYM z`ZQqs)tU|tN^;rO#8mFmsJYtH?| zY@a_+JfIlNH{$CH7EVDkb46HVwJ20ng%Zto!DleiSsC=)Xx0aq-aBhV{dqHQF-|BH zl*MmE%wKixv@SKUtX1PNHOdREQ+_Ktee_pqbLkYUK-9vJR~8NINHF|N>a!yq+luY} z!1el-fYqbEBlxQ*_X;kRio|&!vMb5yX^A+Ci)C_z%Z~tNOylMW_=x3^>>P~Ym*?V%f*m7)RN-UOi~_c>=-)e>e?@lOXLJnT+goU)M?SQn}~zGmrcS)kts?F zFI462H+M9n+lC3D?0Ev;m^%dcJxdlvjjG6?$miDTac6+ZOWb4QNjd z+g#qZUf}NrO7Vb$Enz%Jj+-;(92?JI-}duRJ1n2g%^EZW$(QuhSMOcMXz?h#2N1Q$ z3mq|wgVes~mgOVW^3IKYm)JB6Vb|5E39RpQ*STz#TJ3$|*v&-7)b>Y?PUeI}vs);= zg&R&E;0v5Ab=$zqNuc%I%%*o+qSe#@^YPL61BRQ{&uQ9H?U!D25|i4-VGdW;@-60+ zxynXoqd^aGD3uliAKE!qy}l_pvQp=0Iv9Dj(X_qna!PPBtZ zF%U!B3Ffe9^U%sQ>h>zaR+!0T>NR!QgcH6x&r4Y!&M%3DED}MQpWJiL`41`|SfRWG z`fg|gZ;q+0t*C2+Uf?M@etfB=Pb(SUpO3(y(-Z4*GK@ctSS56w$U(V z)SW>!`Mip5tSw7-lHRxyfNPaHfI&zUnQXS}9!f9i3UM&+bY0BeF9NoBpI$ z;!vBLUjMD6qltR9Jv{C2Blxi7n7ESQ6>SmZWpr{w)f|;*POC@eceX?_sKMS7@7sG- zei=$t@OdL0QaM+xm*Djlo`yBZ>4bX3PvwyiOnEHEnZLd6qBP} zy0<;>TuU`M#dV38IA{j+sM$!VzgEzJ!Bjv?{VeJa7xoJpqMv3poZ^&&kts;l=~BZV zPnSI=*LeSZ_Ne&br*s0dz#Rt=K=TshkdVXRHAGwelGTj&=Mxa&K+*T^3g>=Q3j>rP z7C0LHpU@HI7PK1nc|E=FTm4SJm|};PDep$&;X zO2^aQz|p21n7~GdEk*?ZH3Tb;UO^ZMsOI;W4{Q#x&Ic65-#NAmbyXYPW56p)hl+WB z`g`S~Kb*-ToG<+^dn^6QWbB=@Wz81@QCT;0cNko8 zeERfag@~$jxllwfArK-- z0?e{F01S@NxDF}--+^Di|NqDT*(1n#h44&;EybpDaD>SH@;()LX=SM@36sG80T`rB A&Hw-a literal 0 HcmV?d00001 diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index bacf6ffc6..3c70d770a 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -905,27 +905,27 @@ namespace ReSolve return resnorm / (norm_x * norm_A); } - const std::string SystemSolver::getFactorizationMethod() const + const std::string& SystemSolver::getFactorizationMethod() const { return factorization_method_; } - const std::string SystemSolver::getRefactorizationMethod() const + const std::string& SystemSolver::getRefactorizationMethod() const { return refactorization_method_; } - const std::string SystemSolver::getSolveMethod() const + const std::string& SystemSolver::getSolveMethod() const { return solve_method_; } - const std::string SystemSolver::getRefinementMethod() const + const std::string& SystemSolver::getRefinementMethod() const { return ir_method_; } - const std::string SystemSolver::getGramSchmidtMethod() const + const std::string& SystemSolver::getGramSchmidtMethod() const { return gs_method_; } diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index f660d6ad3..b1a13394f 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -78,11 +78,11 @@ namespace ReSolve real_type getNormOfScaledResiduals(vector_type* rhs, vector_type* x); // Get solver parameters - const std::string getFactorizationMethod() const; - const std::string getRefactorizationMethod() const; - const std::string getSolveMethod() const; - const std::string getRefinementMethod() const; - const std::string getGramSchmidtMethod() const; + const std::string& getFactorizationMethod() const; + const std::string& getRefactorizationMethod() const; + const std::string& getSolveMethod() const; + const std::string& getRefinementMethod() const; + const std::string& getGramSchmidtMethod() const; // Set solver parameters int setFactorizationMethod(std::string method); From a73deb8481e822a8b189b258a61884d4747a618f Mon Sep 17 00:00:00 2001 From: pelesh Date: Sat, 3 Oct 2026 19:51:00 -0400 Subject: [PATCH 12/16] Improve eror management. --- resolve/SystemSolver.cpp | 63 ++++++++++++++++++++++++++-------------- resolve/SystemSolver.hpp | 12 ++++---- 2 files changed, 48 insertions(+), 27 deletions(-) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 3c70d770a..ba6c43e7a 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -239,19 +239,14 @@ namespace ReSolve } // If we use iterative solver, we can set it up here - if (solve_method_ == "randgmres") - { - auto* rgmres = dynamic_cast(iterative_solver_.get()); - status += rgmres->setup(A_); - } - else if (solve_method_ == "fgmres") - { - auto* fgmres = dynamic_cast(iterative_solver_.get()); - status += fgmres->setup(A_); - } - else + if (solve_method_ == "randgmres" || solve_method_ == "fgmres") { - // do nothing + if (iterative_solver_ == nullptr) + { + out::error() << "Iterative solver not initialized!\n"; + return 1; + } + status += iterative_solver_->setup(A_); } return status; @@ -362,6 +357,9 @@ namespace ReSolve factorization_solver_->setup(A_); return factorization_solver_->analyze(); } + + out::error() << "Cannot analyze: factorization method '" << factorization_method_ + << "' does not provide a factorization solver.\n"; return 1; } @@ -372,6 +370,9 @@ namespace ReSolve is_solve_on_device_ = false; return factorization_solver_->factorize(); } + + out::error() << "Cannot factorize: factorization method '" << factorization_method_ + << "' does not provide a factorization solver.\n"; return 1; } @@ -437,7 +438,12 @@ namespace ReSolve { status += refactorization_solver_->setup(A_, L_, U_, P_, Q_); - LinSolverDirectCuSolverRf* Rf = dynamic_cast(refactorization_solver_.get()); + auto* Rf = dynamic_cast(refactorization_solver_.get()); + if (Rf == nullptr) + { + out::error() << "Refactorization solver is not a cuSolverRf instance!\n"; + return 1; + } Rf->setNumericalProperties(1e-14, 1e-1); is_solve_on_device_ = false; @@ -445,7 +451,12 @@ namespace ReSolve #ifdef RESOLVE_USE_CUDSS else if (refactorization_method_ == "cudssrf") { - LinSolverDirectCuDssRf* Rf = dynamic_cast(refactorization_solver_.get()); + auto* Rf = dynamic_cast(refactorization_solver_.get()); + if (Rf == nullptr) + { + out::error() << "Refactorization solver is not a cuDSS instance!\n"; + return 1; + } Rf->setNumericalProperties(1e-14, 1e-1); status += refactorization_solver_->setup(A_, L_, U_, P_, Q_); @@ -757,8 +768,11 @@ namespace ReSolve * * @param[in] method - string ID for the iterative refinement method * @param[in] gs_method - string ID for the orthogonalization method to be used + * + * @return int 0 if the requested method is now active, 1 if it was rejected + * or not recognized */ - void SystemSolver::setRefinementMethod(std::string method, std::string gs_method) + int SystemSolver::setRefinementMethod(std::string method, std::string gs_method) { // With an iterative solve method the Krylov solver belongs to the solve // path, not to iterative refinement, so leave it untouched. @@ -769,9 +783,11 @@ namespace ReSolve out::warning() << "Iterative refinement cannot be enabled together with an " << "iterative solve method ('randgmres' or 'fgmres'). " << "Keeping refinement method 'none'.\n"; + ir_method_ = "none"; + return 1; } ir_method_ = "none"; - return; + return 0; } // Direct solve path: the Krylov solver, if any, is the refinement solver. @@ -781,7 +797,7 @@ namespace ReSolve if (method == "none") { - return; + return 0; } gs_method_ = gs_method; @@ -789,11 +805,11 @@ namespace ReSolve if (method == "fgmres" && createIterativeSolver(method) == 0) { ir_method_ = method; + return 0; } - else - { - out::error() << "Iterative refinement method " << method << " not recognized.\n"; - } + + out::error() << "Iterative refinement method " << method << " not recognized.\n"; + return 1; } real_type SystemSolver::getVectorNorm(vector_type* rhs) @@ -958,6 +974,11 @@ namespace ReSolve { // TODO: Use cast here as a temporary solution; will be replaced by parameter setting framework auto* sol = dynamic_cast(iterative_solver_.get()); + if (sol == nullptr) + { + out::error() << "Iterative solver is not a randomized FGMRES instance!\n"; + return 1; + } sol->setSketchingMethod(tmp); } diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index b1a13394f..b19eceacc 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -85,12 +85,12 @@ namespace ReSolve const std::string& getGramSchmidtMethod() const; // Set solver parameters - int setFactorizationMethod(std::string method); - int setRefactorizationMethod(std::string method); - int setSolveMethod(std::string method); - void setRefinementMethod(std::string method, std::string gs = "cgs2"); - int setSketchingMethod(std::string method); - int setGramSchmidtMethod(std::string gs_method); + int setFactorizationMethod(std::string method); + int setRefactorizationMethod(std::string method); + int setSolveMethod(std::string method); + int setRefinementMethod(std::string method, std::string gs = "cgs2"); + int setSketchingMethod(std::string method); + int setGramSchmidtMethod(std::string gs_method); private: template From 4b88721ac43df847a833d6282227043d7b4b2354 Mon Sep 17 00:00:00 2001 From: pelesh Date: Sun, 4 Oct 2026 13:00:56 -0400 Subject: [PATCH 13/16] Fix potential memory leak in LinSolverDirectKLU class. --- resolve/LinSolverDirectKLU.cpp | 87 +++++++++++++--------------------- resolve/LinSolverDirectKLU.hpp | 1 + 2 files changed, 35 insertions(+), 53 deletions(-) diff --git a/resolve/LinSolverDirectKLU.cpp b/resolve/LinSolverDirectKLU.cpp index 5c2920b31..8ac1f5643 100644 --- a/resolve/LinSolverDirectKLU.cpp +++ b/resolve/LinSolverDirectKLU.cpp @@ -58,17 +58,7 @@ namespace ReSolve */ LinSolverDirectKLU::~LinSolverDirectKLU() { - if (factors_extracted_) - { - delete L_; - delete U_; - delete[] P_; - delete[] Q_; - L_ = nullptr; - U_ = nullptr; - P_ = nullptr; - Q_ = nullptr; - } + freeExtractedFactors(); klu_free_symbolic(&Symbolic_, &Common_); klu_free_numeric(&Numeric_, &Common_); } @@ -106,23 +96,11 @@ namespace ReSolve { klu_free_symbolic(&Symbolic_, &Common_); } - Symbolic_ = klu_analyze(A_->getNumRows(), + Symbolic_ = klu_analyze(A_->getNumRows(), A_->getRowData(memory::HOST), A_->getColData(memory::HOST), &Common_); - factors_extracted_ = false; - - if (L_ != nullptr) - { - delete L_; - L_ = nullptr; - } - - if (U_ != nullptr) - { - delete U_; - U_ = nullptr; - } + freeExtractedFactors(); if (Symbolic_ == nullptr) { @@ -151,19 +129,7 @@ namespace ReSolve Symbolic_, &Common_); - factors_extracted_ = false; - - if (L_ != nullptr) - { - delete L_; - L_ = nullptr; - } - - if (U_ != nullptr) - { - delete U_; - U_ = nullptr; - } + freeExtractedFactors(); if (Numeric_ == nullptr) { @@ -193,19 +159,7 @@ namespace ReSolve Numeric_, &Common_); - factors_extracted_ = false; - - if (L_ != nullptr) - { - delete L_; - L_ = nullptr; - } - - if (U_ != nullptr) - { - delete U_; - U_ = nullptr; - } + freeExtractedFactors(); if (!kluStatus) { @@ -349,8 +303,11 @@ namespace ReSolve { if (Numeric_ != nullptr) { - P_ = new index_type[A_->getNumRows()]; size_t nrows = static_cast(A_->getNumRows()); + if (P_ == nullptr) + { + P_ = new index_type[nrows]; + } std::memcpy(P_, Symbolic_->Q, nrows * sizeof(index_type)); // KLU's CSC Symbolic_->Q is the CSR P vector. // Only a symbolic factorization is needed to get Q, because there is only row pivoting for the numeric factorization. return P_; @@ -373,8 +330,11 @@ namespace ReSolve { if (Numeric_ != nullptr) { - Q_ = new index_type[A_->getNumRows()]; size_t nrows = static_cast(A_->getNumRows()); + if (Q_ == nullptr) + { + Q_ = new index_type[nrows]; + } std::memcpy(Q_, Numeric_->Pnum, nrows * sizeof(index_type)); // KLU's CSC Numeric_->Pnum is the CSR Q vector. // A numeric factorization is needed to get Pnum, because there is row pivoting for the numeric factorization. return Q_; @@ -595,4 +555,25 @@ namespace ReSolve params_list_["halt_if_singular"] = HALT_IF_SINGULAR; } + /** + * @brief Releases factors and permutation vectors extracted from KLU. + * + * L and U factors as well as copies of permutation vectors P and Q are + * allocated on demand by the getters and owned by this object. They are + * invalidated whenever a new symbolic or numeric factorization is + * computed and must be released then and in the destructor. + */ + void LinSolverDirectKLU::freeExtractedFactors() + { + delete L_; + delete U_; + delete[] P_; + delete[] Q_; + L_ = nullptr; + U_ = nullptr; + P_ = nullptr; + Q_ = nullptr; + factors_extracted_ = false; + } + } // namespace ReSolve diff --git a/resolve/LinSolverDirectKLU.hpp b/resolve/LinSolverDirectKLU.hpp index 69ac6e7c0..b9fb6020c 100644 --- a/resolve/LinSolverDirectKLU.hpp +++ b/resolve/LinSolverDirectKLU.hpp @@ -100,6 +100,7 @@ namespace ReSolve private: void initParamList(); + void freeExtractedFactors(); bool factors_extracted_{false}; klu_common Common_; // settings klu_symbolic* Symbolic_{nullptr}; From 4f1853d3db4c6a1ca7ac5b5a4521abe0d51f24af Mon Sep 17 00:00:00 2001 From: pelesh Date: Sun, 4 Oct 2026 18:51:13 -0400 Subject: [PATCH 14/16] Some code comments updates. --- resolve/SystemSolver.cpp | 107 +++++++++++++++++++++++++-------------- 1 file changed, 70 insertions(+), 37 deletions(-) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index ba6c43e7a..a279eef08 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -135,33 +135,6 @@ namespace ReSolve #endif } // namespace - /** - * @brief Shared tail of all constructors. - * - * Creates matrix and vector handlers for the given workspace, derives the - * memory space ID from the workspace type and instantiates solver - * components via `initialize()`. - * - * @tparam Workspace - one of LinAlgWorkspaceCpu, LinAlgWorkspaceCUDA, LinAlgWorkspaceHIP - * @param[in] workspace - pointer to the workspace (not owned) - */ - template - void SystemSolver::completeSetup(Workspace* workspace) - { - matrix_handler_.reset(new MatrixHandler(workspace)); - vector_handler_.reset(new VectorHandler(workspace)); - memspace_ = memorySpaceName(workspace); - - if (initialize() != 0) - { - out::error() << "SystemSolver initialization failed with factorization '" << factorization_method_ - << "', refactorization '" << refactorization_method_ - << "', solve '" << solve_method_ - << "', preconditioner '" << precondition_method_ - << "', iterative refinement '" << ir_method_ - << "'. Solver is not usable in this state.\n"; - } - } SystemSolver::SystemSolver(LinAlgWorkspaceCpu* workspace_cpu, std::string factor, @@ -407,9 +380,7 @@ namespace ReSolve * solver to run on GPU. * * @pre Factorization solver exists and provides access to L and U factors, - * as well as left and right permutation vectors P and Q. Since KLU is - * the only factorization solver available through ReSolve, the factors are - * expected in CSC format. + * as well as left and right permutation vectors P and Q. * * @return int 0 if successful, 1 if it fails */ @@ -812,6 +783,15 @@ namespace ReSolve return 1; } + /** + * @brief Compute L2 vector norm. + * + * @param[in] rhs - pointer to the vector. + * + * @pre `rhs` is not a multivector. + * @invariant `rhs` is unmodified by this method. + * @return L2 norm of vector `rhs`. + */ real_type SystemSolver::getVectorNorm(vector_type* rhs) { using namespace ReSolve::constants; @@ -819,7 +799,7 @@ namespace ReSolve if (memspace_ == "cpu") { norm_b = std::sqrt(vector_handler_->dot(rhs, rhs, memory::HOST)); -#if defined(RESOLVE_USE_HIP) || defined(RESOLVE_USE_CUDA) +#ifdef RESOLVE_USE_GPU } else if (memspace_ == "cuda" || memspace_ == "hip") { @@ -831,7 +811,7 @@ namespace ReSolve { norm_b = std::sqrt(vector_handler_->dot(rhs, rhs, memory::HOST)); } -#endif +#endif // RESOLVE_USE_GPU } else { @@ -841,6 +821,22 @@ namespace ReSolve return norm_b; } + /** + * @brief Compute relative residual norm. + * + * The method evaluates ||A x - b|| / ||b||, where ||...|| is L2 norm, + * for a linear system A x = b, where A is n x n matrix, x is unknown + * n x 1 vector, and b is known right-hand-side vector. When the system + * is solved for x, the relative residual norm provides a measure of solution + * consistency. + * + * @param[in] rhs - right-hand side vector of a linear system. + * @param[in] x - solution vector of a linear system. + * @invariant Matrix `A_` and vectors `x` and `rhs` are not modified by this + * method. + * + * @return Relative residual norm + */ real_type SystemSolver::getResidualNorm(vector_type* rhs, vector_type* x) { using namespace ReSolve::constants; @@ -852,7 +848,7 @@ namespace ReSolve { res_vector_->copyFromExternal(rhs, memory::HOST, memory::HOST); norm_b = std::sqrt(vector_handler_->dot(res_vector_.get(), res_vector_.get(), memory::HOST)); -#if defined(RESOLVE_USE_HIP) || defined(RESOLVE_USE_CUDA) +#ifdef RESOLVE_USE_GPU } else if (memspace_ == "cuda" || memspace_ == "hip") { @@ -866,10 +862,9 @@ namespace ReSolve res_vector_->copyFromExternal(rhs, memory::HOST, memory::DEVICE); res_vector_->syncData(memory::HOST); norm_b = std::sqrt(vector_handler_->dot(res_vector_.get(), res_vector_.get(), memory::HOST)); - // ms = memory::HOST; } ms = memory::DEVICE; -#endif +#endif // RESOLVE_USE_GPU } else { @@ -882,6 +877,16 @@ namespace ReSolve return resnorm / norm_b; } + /** + * @brief Norm of scaled residuals. + * + * @param[in] rhs - right-hand side vector of a linear system. + * @param[in] x - solution vector of a linear system. + * @invariant Matrix `A_` and vectors `x` and `rhs` are not modified by this + * method. + * + * @return Norm of scaled residuals. + */ real_type SystemSolver::getNormOfScaledResiduals(vector_type* rhs, vector_type* x) { using namespace ReSolve::constants; @@ -893,7 +898,7 @@ namespace ReSolve if (memspace_ == "cpu") { res_vector_->copyFromExternal(rhs, memory::HOST, memory::HOST); -#if defined(RESOLVE_USE_HIP) || defined(RESOLVE_USE_CUDA) +#ifdef RESOLVE_USE_GPU } else if (memspace_ == "cuda" || memspace_ == "hip") { @@ -906,7 +911,7 @@ namespace ReSolve res_vector_->copyFromExternal(rhs, memory::HOST, memory::DEVICE); } ms = memory::DEVICE; -#endif +#endif // RESOLVE_USE_GPU } else { @@ -1021,6 +1026,34 @@ namespace ReSolve // Private methods // + /** + * @brief Shared tail of all constructors. + * + * Creates matrix and vector handlers for the given workspace, derives the + * memory space ID from the workspace type and instantiates solver + * components via `initialize()`. + * + * @tparam Workspace - one of LinAlgWorkspaceCpu, LinAlgWorkspaceCUDA, LinAlgWorkspaceHIP + * @param[in] workspace - pointer to the workspace (not owned) + */ + template + void SystemSolver::completeSetup(Workspace* workspace) + { + matrix_handler_.reset(new MatrixHandler(workspace)); + vector_handler_.reset(new VectorHandler(workspace)); + memspace_ = memorySpaceName(workspace); + + if (initialize() != 0) + { + out::error() << "SystemSolver initialization failed with factorization '" << factorization_method_ + << "', refactorization '" << refactorization_method_ + << "', solve '" << solve_method_ + << "', preconditioner '" << precondition_method_ + << "', iterative refinement '" << ir_method_ + << "'. Solver is not usable in this state.\n"; + } + } + /** * @brief Enforces supported combinations of user-selected methods. * From 22932506f7b1c897aa5637f3021529bdd0f83bfe Mon Sep 17 00:00:00 2001 From: pelesh Date: Sun, 4 Oct 2026 22:51:44 +0000 Subject: [PATCH 15/16] Apply pre-commmit fixes --- resolve/SystemSolver.cpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index a279eef08..d84c4444c 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -135,7 +135,6 @@ namespace ReSolve #endif } // namespace - SystemSolver::SystemSolver(LinAlgWorkspaceCpu* workspace_cpu, std::string factor, std::string refactor, @@ -785,9 +784,9 @@ namespace ReSolve /** * @brief Compute L2 vector norm. - * + * * @param[in] rhs - pointer to the vector. - * + * * @pre `rhs` is not a multivector. * @invariant `rhs` is unmodified by this method. * @return L2 norm of vector `rhs`. @@ -823,18 +822,18 @@ namespace ReSolve /** * @brief Compute relative residual norm. - * + * * The method evaluates ||A x - b|| / ||b||, where ||...|| is L2 norm, * for a linear system A x = b, where A is n x n matrix, x is unknown * n x 1 vector, and b is known right-hand-side vector. When the system * is solved for x, the relative residual norm provides a measure of solution * consistency. - * + * * @param[in] rhs - right-hand side vector of a linear system. * @param[in] x - solution vector of a linear system. * @invariant Matrix `A_` and vectors `x` and `rhs` are not modified by this * method. - * + * * @return Relative residual norm */ real_type SystemSolver::getResidualNorm(vector_type* rhs, vector_type* x) @@ -879,12 +878,12 @@ namespace ReSolve /** * @brief Norm of scaled residuals. - * + * * @param[in] rhs - right-hand side vector of a linear system. * @param[in] x - solution vector of a linear system. * @invariant Matrix `A_` and vectors `x` and `rhs` are not modified by this * method. - * + * * @return Norm of scaled residuals. */ real_type SystemSolver::getNormOfScaledResiduals(vector_type* rhs, vector_type* x) From 17e442ba0621d7771a7c19c70c31801d6e383e2d Mon Sep 17 00:00:00 2001 From: Shaked Regev <35384901+shakedregev@users.noreply.github.com> Date: Mon, 5 Oct 2026 10:55:21 -0400 Subject: [PATCH 16/16] [skip ci] Clarify comments Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com> --- resolve/LinSolverDirectKLU.cpp | 8 ++++---- resolve/SystemSolver.cpp | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/resolve/LinSolverDirectKLU.cpp b/resolve/LinSolverDirectKLU.cpp index 8ac1f5643..2aff43fe1 100644 --- a/resolve/LinSolverDirectKLU.cpp +++ b/resolve/LinSolverDirectKLU.cpp @@ -558,10 +558,10 @@ namespace ReSolve /** * @brief Releases factors and permutation vectors extracted from KLU. * - * L and U factors as well as copies of permutation vectors P and Q are - * allocated on demand by the getters and owned by this object. They are - * invalidated whenever a new symbolic or numeric factorization is - * computed and must be released then and in the destructor. + * L and U factors and copies of permutation vectors P and Q are + * allocated on demand by the getters and owned by this object. The old factors + * are invalidated when a new symbolic or numeric factorization is + * computed and then must be released in the destructor. */ void LinSolverDirectKLU::freeExtractedFactors() { diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index d84c4444c..b5b5b381d 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -44,16 +44,16 @@ namespace ReSolve // Create a shortcut name for Logger static class using out = io::Logger; - // Helpers below map user-facing string IDs to solver enums. They are kept + // The helpers below map user-facing string IDs to solver enums. They are kept // in an anonymous namespace (still inside ReSolve) so they have internal - // linkage: they are implementation details of SystemSolver, are not part + // linkage. This means they are implementation details of SystemSolver, are not part // of the library's exported symbols, and cannot collide with same-named // functions in other translation units. If these mappings become useful // elsewhere, move them to the classes that own the enums (GramSchmidt and // LinSolverIterativeRandFGMRES) as public static methods. namespace { - /// Maps string ID to the sketching method enum; warns and defaults to count sketch. + /// Maps string ID to the sketching method enum, warns and defaults to count sketch. LinSolverIterativeRandFGMRES::SketchingMethod sketchingMethodFromString(const std::string& method) { if (method == "count") @@ -69,7 +69,7 @@ namespace ReSolve return LinSolverIterativeRandFGMRES::cs; } - /// Maps string ID to the Gram-Schmidt variant enum; warns and defaults to CGS2. + /// Maps string ID to the Gram-Schmidt variant enum, warns and defaults to CGS2. GramSchmidt::GSVariant gsVariantFromString(const std::string& variant) { if (variant == "cgs2")