diff --git a/CHANGELOG.md b/CHANGELOG.md index ad503b2e..99259977 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -80,4 +80,6 @@ It is seamless from the user perspective and fixed many bugs. 16. Updated MatrixHandler::addConst to return integer error codes instead of void. -17. Added a preconditioner interface class so users can define thier own preconditioners. +17. Added a preconditioner interface class so users can define their own preconditioners. + +18. Added left preconditioning support for GMRES and a user-defined preconditioner class. diff --git a/examples/ExampleHelper.hpp b/examples/ExampleHelper.hpp index 3d569371..50f5716b 100644 --- a/examples/ExampleHelper.hpp +++ b/examples/ExampleHelper.hpp @@ -141,12 +141,12 @@ namespace ReSolve } res_->copyFromExternal(r_, memspace_, memspace_); - real_type norm = computeResidualNorm(*A_, *x_, *res_, memspace_); - real_type rnorm = norm2(*r_, memspace_); + real_type norm = computeResidualNorm(*A_, *x_, *res_, memspace_); + real_type rhs_norm = norm2(*r_, memspace_); std::cout << "\t2-Norm of the residual: " << std::scientific << std::setprecision(16) - << norm / rnorm << "\n"; + << norm / rhs_norm << "\n"; } /// Summary of direct solve diff --git a/examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp b/examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp index 4d67d49e..70e4298d 100644 --- a/examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp +++ b/examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp @@ -186,8 +186,8 @@ int main(int argc, char* argv[]) << status << std::endl; vec_rhs->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); status = Rf->solve(vec_rhs, vec_x); - ReSolve::PreconditionerLU precond_lu(Rf); - FGMRES->setPreconditioner(&precond_lu); + ReSolve::PreconditionerLU preconditioner(Rf); + FGMRES->setPreconditioner(&preconditioner); } // if (i%2!=0) vec_x->setToZero(ReSolve::memory::DEVICE); real_type norm_x = vector_handler->dot(vec_x, vec_x, ReSolve::memory::DEVICE); diff --git a/examples/gpuRefactor.cpp b/examples/gpuRefactor.cpp index ae6c366a..17516fcd 100644 --- a/examples/gpuRefactor.cpp +++ b/examples/gpuRefactor.cpp @@ -297,8 +297,8 @@ int gpuRefactor(int argc, char* argv[]) { // Setup iterative refinement FGMRES.resetMatrix(A); - ReSolve::PreconditionerLU precond_lu(&Rf); - FGMRES.setPreconditioner(&precond_lu); + ReSolve::PreconditionerLU preconditioner(&Rf); + FGMRES.setPreconditioner(&preconditioner); // If refactorization produced finite solution do iterative refinement if (std::isfinite(helper.getNormRelativeResidual())) diff --git a/examples/kluFactor.cpp b/examples/kluFactor.cpp index 615ea504..0e9f6b15 100644 --- a/examples/kluFactor.cpp +++ b/examples/kluFactor.cpp @@ -191,8 +191,8 @@ int main(int argc, char* argv[]) { // Setup iterative refinement FGMRES.setup(A); - ReSolve::PreconditionerLU precond_lu(&KLU); - FGMRES.setPreconditioner(&precond_lu); + ReSolve::PreconditionerLU preconditioner(&KLU); + FGMRES.setPreconditioner(&preconditioner); // If refactorization produced finite solution do iterative refinement if (std::isfinite(helper.getNormRelativeResidual())) diff --git a/examples/kluRefactor.cpp b/examples/kluRefactor.cpp index 5f68efde..c080b93d 100644 --- a/examples/kluRefactor.cpp +++ b/examples/kluRefactor.cpp @@ -194,8 +194,8 @@ int main(int argc, char* argv[]) { // Setup iterative refinement FGMRES.setup(A); - ReSolve::PreconditionerLU precond_lu(KLU); - FGMRES.setPreconditioner(&precond_lu); + ReSolve::PreconditionerLU preconditioner(KLU); + FGMRES.setPreconditioner(&preconditioner); // If refactorization produced finite solution do iterative refinement if (std::isfinite(helper.getNormRelativeResidual())) diff --git a/examples/randGmres.cpp b/examples/randGmres.cpp index f8ec86a5..0a46d8da 100644 --- a/examples/randGmres.cpp +++ b/examples/randGmres.cpp @@ -116,7 +116,7 @@ int runGmresExample(int argc, char* argv[]) GramSchmidt GS(&vector_handler, GramSchmidt::CGS2); - precon_type Precond(&workspace); + precon_type precondition_solver(&workspace); LinSolverIterativeRandFGMRES FGMRES(&matrix_handler, &vector_handler, LinSolverIterativeRandFGMRES::cs, @@ -170,15 +170,16 @@ int runGmresExample(int argc, char* argv[]) matrix_handler.setValuesChanged(true, memspace); - Precond.setup(A); FGMRES.setRestart(150); FGMRES.setMaxit(2500); FGMRES.setTol(1e-12); FGMRES.setup(A); - FGMRES.resetMatrix(A); - ReSolve::PreconditionerLU precond_lu(&Precond); - FGMRES.setPreconditioner(&precond_lu); + + ReSolve::PreconditionerLU preconditioner(&precondition_solver); + preconditioner.setup(A); + + FGMRES.setPreconditioner(&preconditioner); FGMRES.setFlexible(1); FGMRES.solve(vec_rhs, vec_x); diff --git a/examples/sysGmres.cpp b/examples/sysGmres.cpp index faecbd01..538510ed 100644 --- a/examples/sysGmres.cpp +++ b/examples/sysGmres.cpp @@ -35,6 +35,7 @@ void printHelpInfo() std::cout << "\t-g \tGram-Schmidt method: cgs1, cgs2, or mgs (default 'cgs2').\n"; std::cout << "\t-s \tSketching method: count or fwht (default 'count')\n"; std::cout << "\t-x \tEnable flexible: yes or no (default 'yes')\n\n"; + std::cout << "\t-p \tPreconditioner side: left or right (default 'right')\n\n"; } // @@ -56,7 +57,8 @@ static int sysGmres(int argc, char* argv[]); static void processInputs(std::string& method, std::string& gs, std::string& sketch, - std::string& flexible); + std::string& flexible, + std::string& side); /// Main function selects example to be run int main(int argc, char* argv[]) @@ -164,7 +166,10 @@ int sysGmres(int argc, char* argv[]) opt = options.getParamFromKey("-x"); std::string flexible = opt ? (*opt).second : "yes"; - processInputs(method, gs, sketch, flexible); + opt = options.getParamFromKey("-p"); + std::string side = opt ? (*opt).second : "right"; + + processInputs(method, gs, sketch, flexible, side); std::cout << "Matrix file: " << matrix_pathname << "\n" << "RHS file: " << rhs_pathname << "\n"; @@ -247,7 +252,7 @@ int sysGmres(int argc, char* argv[]) // Set up the preconditioner if (return_code == 0) { - status = solver.preconditionerSetup(); + status = solver.preconditionerSetup(side); std::cout << "solver.preconditionerSetup returned status: " << status << "\n"; if (status != 0) { @@ -281,7 +286,7 @@ int sysGmres(int argc, char* argv[]) return return_code; } -void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& flexible) +void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& flexible, std::string& side) { if (method == "randgmres") { @@ -314,4 +319,10 @@ void processInputs(std::string& method, std::string& gs, std::string& sketch, st std::cout << "Setting flexible to the default (yes).\n\n"; flexible = "yes"; } + + if ((side != "left") && (side != "right")) + { + std::cout << "Preconditioning side " << side << " not recognized.\n"; + std::cout << "Setting preconditioning side to the default (right).\n\n"; + } } diff --git a/resolve/CMakeLists.txt b/resolve/CMakeLists.txt index 628963c8..d54644ca 100644 --- a/resolve/CMakeLists.txt +++ b/resolve/CMakeLists.txt @@ -21,6 +21,7 @@ set(ReSolve_SRC SystemSolver.cpp Preconditioner.cpp PreconditionerLU.cpp + PreconditionerUserMatrix.cpp ) set(ReSolve_KLU_SRC LinSolverDirectKLU.cpp) @@ -52,6 +53,7 @@ set(ReSolve_HEADER_INSTALL MemoryUtils.hpp Preconditioner.hpp PreconditionerLU.hpp + PreconditionerUserMatrix.hpp ) set(ReSolve_KLU_HEADER_INSTALL LinSolverDirectKLU.hpp) diff --git a/resolve/LinSolverIterativeFGMRES.cpp b/resolve/LinSolverIterativeFGMRES.cpp index 67efc53e..144ba74e 100644 --- a/resolve/LinSolverIterativeFGMRES.cpp +++ b/resolve/LinSolverIterativeFGMRES.cpp @@ -103,10 +103,39 @@ namespace ReSolve return 0; } + /** + * @brief Solve linear system A*x = rhs + * + * Implements restarted GMRES with optional flexible (FGMRES) variant. + * + * Flexible GMRES allows the preconditioner to vary periteration and + * uses right preconditioning. Standard GMRES supports both left and + * right preconditioning. + * + * Left preconditioning solves M^{-1}Ax = M^{-1}b and checks convergence + * with ||M^{-1}(b - Ax)||. Right preconditioning solves AM^{-1}(Mx) = b + * and checks convergence with ||b - Ax||. Both report the true relative + * residual ||b - Ax||/||b||. + * + * @param rhs - right hand side vector + * @param x - solution vector + * @return int - zero if successful, error code otherwise + * + * @invariant rhs vector is unchanged. + * @post x is overwritten with the solution to the linear system. + */ int LinSolverIterativeFGMRES::solve(vector_type* rhs, vector_type* x) { using namespace constants; + // Flexible GMRES only supports right preconditioning + if (flexible_ && preconditioner_->getSide() == Preconditioner::Side::LEFT) + { + out::warning() << "Flexible GMRES does not support left preconditioning. " + << "Switching to right preconditioning.\n"; + preconditioner_->setSide(Preconditioner::Side::RIGHT); + } + // io::Logger::setVerbosity(io::Logger::EVERYTHING); int outer_flag = 1; @@ -117,34 +146,73 @@ namespace ReSolve int k = 0; int k1 = 0; - real_type t = 0.0; - real_type rnorm = 0.0; - real_type bnorm = 0.0; + real_type t = 0.0; + real_type res_norm = 0.0; // Residual norm used used for convergence + real_type rhs_norm = 0.0; // Right-hand side norm used for convergence + real_type true_res_norm = 0.0; // True (unpreconditioned) residual norm ||b - Ax|| for reporting + real_type true_rhs_norm = 0.0; // True (unpreconditioned) right-hand side norm ||b|| for reporting real_type tolrel; vector_type vec_v(n_); vector_type vec_z(n_); - // V[0] = b-A*x_0 - // debug + + // Compute initial residual norm. + // V[0] = ||b - A*x0|| for right preconditioning + // V[0] = ||M^{-1}{b - A*x0}|| for left preconditioning + vec_Z_->setToZero(memspace_); vec_V_->setToZero(memspace_); rhs->copyToExternal(vec_V_->getData(memspace_), 0, memspace_, memspace_); matrix_handler_->matvec(A_, x, vec_V_, &MINUS_ONE, &ONE, memspace_); - rnorm = 0.0; - bnorm = vector_handler_->dot(rhs, rhs, memspace_); - rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_); - // rnorm = ||V_1|| - rnorm = std::sqrt(rnorm); - bnorm = std::sqrt(bnorm); + + vec_v.setData(vec_V_->getData(0, memspace_), memspace_); + vec_z.setData(vec_Z_->getData(0, memspace_), memspace_); + + // True residual norm ||b - A*x0|| + true_res_norm = vector_handler_->dot(&vec_v, &vec_v, memspace_); + true_res_norm = std::sqrt(true_res_norm); + + // True right-hand side norm ||b|| + true_rhs_norm = vector_handler_->dot(rhs, rhs, memspace_); + true_rhs_norm = std::sqrt(true_rhs_norm); + + switch (preconditioner_->getSide()) + { + case Preconditioner::Side::RIGHT: + // Right preconditioning uses true norms for convergence + res_norm = true_res_norm; + rhs_norm = true_rhs_norm; + break; + case Preconditioner::Side::LEFT: + // Left preconditioning uses preconditioned norms for convergence + // Left-preconditioned residual norm ||M^{-1}*(b-A*x0)|| + preconditioner_->apply(&vec_v, &vec_z); + vec_v.copyFromExternal(&vec_z, memspace_, memspace_); + res_norm = vector_handler_->dot(vec_V_, vec_V_, memspace_); + res_norm = std::sqrt(res_norm); + + // Left-preconditioned right-hand side norm ||M^{-1}*b|| + preconditioner_->apply(rhs, &vec_z); + rhs_norm = vector_handler_->dot(&vec_z, &vec_z, memspace_); + rhs_norm = std::sqrt(rhs_norm); + break; + default: + out::error() << "Unknown preconditioner side.\n"; + return 1; + } + io::Logger::misc() << "it 0: norm of residual " << std::scientific << std::setprecision(16) - << rnorm << " Norm of rhs: " << bnorm << "\n"; - initial_residual_norm_ = rnorm / bnorm; // relative residual norm + << res_norm << " Norm of rhs: " << rhs_norm << "\n"; + + // Report the true initial relative residual norm + initial_residual_norm_ = true_res_norm / true_rhs_norm; + while (outer_flag) { if (it == 0) { - tolrel = tol_ * rnorm; + tolrel = tol_ * res_norm; if (std::abs(tolrel) < MACHINE_EPSILON) { tolrel = MACHINE_EPSILON; @@ -155,30 +223,30 @@ namespace ReSolve switch (conv_cond_) { case 0: - exit_cond = ((std::abs(rnorm - ZERO) <= MACHINE_EPSILON)); + exit_cond = ((std::abs(res_norm - ZERO) <= MACHINE_EPSILON)); break; case 1: - exit_cond = ((std::abs(rnorm - ZERO) <= MACHINE_EPSILON) || (rnorm < tol_)); + exit_cond = ((std::abs(res_norm - ZERO) <= MACHINE_EPSILON) || (res_norm < tol_)); break; case 2: - exit_cond = ((std::abs(rnorm - ZERO) <= MACHINE_EPSILON) || (rnorm < (tol_ * bnorm))); + exit_cond = ((std::abs(res_norm - ZERO) <= MACHINE_EPSILON) || (res_norm < (tol_ * rhs_norm))); break; } if (exit_cond) { outer_flag = 0; - final_residual_norm_ = rnorm; - initial_residual_norm_ = rnorm; + final_residual_norm_ = res_norm; + initial_residual_norm_ = res_norm; total_iters_ = 0; break; } // normalize first vector - t = 1.0 / rnorm; + t = 1.0 / res_norm; vector_handler_->scal(t, vec_V_, memspace_); // initialize norm history - h_rs_[0] = rnorm; + h_rs_[0] = res_norm; i = -1; notconv = 1; @@ -187,7 +255,6 @@ namespace ReSolve i++; it++; - // Z_i = (LU)^{-1}*V_i vec_v.setData(vec_V_->getData(i, memspace_), memspace_); if (flexible_) { @@ -197,14 +264,35 @@ namespace ReSolve { vec_z.setData(vec_Z_->getData(0, memspace_), memspace_); } - this->precV(&vec_v, &vec_z); - mem_.deviceSynchronize(); - - // V_{i+1}=A*Z_i - vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_); + // Expand the Krylov subspace. + // + // New basis vector: + // V[i+1] = A*M^{-1}*V[i] (right preconditioning) + // V[i+1] = M^{-1}*A*V[i] (left preconditioning) - matrix_handler_->matvec(A_, &vec_z, &vec_v, &ONE, &ZERO, memspace_); + switch (preconditioner_->getSide()) + { + case Preconditioner::Side::RIGHT: + // Compute vec_z = M^{-1}*V[i], then V[i+1] = A*vec_z + preconditioner_->apply(&vec_v, &vec_z); + mem_.deviceSynchronize(); + + vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_); + matrix_handler_->matvec(A_, &vec_z, &vec_v, &ONE, &ZERO, memspace_); + break; + case Preconditioner::Side::LEFT: + // Compute vec_z = A*V[i], then V[i+1] = M^{-1}*vec_z + matrix_handler_->matvec(A_, &vec_v, &vec_z, &ONE, &ZERO, memspace_); + + vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_); + preconditioner_->apply(&vec_z, &vec_v); + mem_.deviceSynchronize(); + break; + default: + out::error() << "Unknown preconditioner side.\n"; + return 1; + } // orthogonalize V[i+1], form a column of h_H_ @@ -238,12 +326,12 @@ namespace ReSolve h_H_[(i) * (restart_ + 1) + (i + 1)] = h_c_[i] * Hii1 - h_s_[i] * Hii; // residual norm estimate - rnorm = std::abs(h_rs_[i + 1]); + res_norm = std::abs(h_rs_[i + 1]); io::Logger::misc() << "it: " << it << " --> norm of the residual " << std::scientific << std::setprecision(16) - << rnorm << "\n"; + << res_norm << "\n"; // check convergence - if (i + 1 >= restart_ || rnorm <= tolrel || it >= maxit_) + if (i + 1 >= restart_ || res_norm <= tolrel || it >= maxit_) { notconv = 0; } @@ -251,7 +339,7 @@ namespace ReSolve io::Logger::misc() << "End of cycle, ESTIMATED norm of residual " << std::scientific << std::setprecision(16) - << rnorm << "\n"; + << res_norm << "\n"; // solve tri system h_rs_[i] = h_rs_[i] / h_H_[i * (restart_ + 1) + i]; for (int ii = 2; ii <= i + 1; ii++) @@ -266,7 +354,11 @@ namespace ReSolve h_rs_[k] = t / h_H_[k * (restart_ + 1) + k]; } - // get solution + // Update the approximate solution x using h_rs_. + // Flexible GMRES uses the preconditioned basis Z[j] directly. + // Standard GMRES first forms vec_z from V[j], then applies M^{-1} + // only for right preconditioning. + if (flexible_) { for (j = 0; j <= i; j++) @@ -277,6 +369,7 @@ namespace ReSolve } else { + // Accumulate the correction vec_z = sum_j h_rs_[j] * V[j] vec_Z_->setToZero(memspace_); vec_z.setData(vec_Z_->getData(0, memspace_), memspace_); for (j = 0; j <= i; j++) @@ -284,37 +377,61 @@ namespace ReSolve vec_v.setData(vec_V_->getData(j, memspace_), memspace_); vector_handler_->axpy(h_rs_[j], &vec_v, &vec_z, memspace_); } - // now multiply d_Z by precon - - vec_v.setData(vec_V_->getData(memspace_), memspace_); - this->precV(&vec_z, &vec_v); - // and add to x - vector_handler_->axpy(ONE, &vec_v, x, memspace_); + // Apply the correction to x based on preconditioning side + switch (preconditioner_->getSide()) + { + case Preconditioner::Side::RIGHT: + // Right preconditioning: x += M^{-1} * vec_z + preconditioner_->apply(&vec_z, &vec_v); + vector_handler_->axpy(ONE, &vec_v, x, memspace_); + break; + case Preconditioner::Side::LEFT: + // Left preconditioning: x += vec_z + vector_handler_->axpy(ONE, &vec_z, x, memspace_); + break; + default: + out::error() << "Unknown preconditioner side.\n"; + return 1; + } } /* test solution */ - if (rnorm <= tolrel || it >= maxit_) + if (res_norm <= tolrel || it >= maxit_) { - // rnorm_aux = rnorm; + // res_norm_aux = res_norm; outer_flag = 0; } rhs->copyToExternal(vec_V_->getData(memspace_), 0, memspace_, memspace_); matrix_handler_->matvec(A_, x, vec_V_, &MINUS_ONE, &ONE, memspace_); - rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_); - // rnorm = ||V_1|| - rnorm = std::sqrt(rnorm); + + vec_v.setData(vec_V_->getData(0, memspace_), memspace_); + true_res_norm = vector_handler_->dot(&vec_v, &vec_v, memspace_); + true_res_norm = std::sqrt(true_res_norm); + + // Left-preconditioned GMRES applies M^{-1} to the residual + if (preconditioner_->getSide() == Preconditioner::Side::LEFT) + { + preconditioner_->apply(&vec_v, &vec_z); + vec_v.copyFromExternal(&vec_z, memspace_, memspace_); + } + res_norm = vector_handler_->dot(vec_V_, vec_V_, memspace_); + + // res_norm = ||V_1|| + res_norm = std::sqrt(res_norm); if (!outer_flag) { - final_residual_norm_ = rnorm / bnorm; // relative residual norm + // Report the true relative residual norm + final_residual_norm_ = true_res_norm / true_rhs_norm; total_iters_ = it; io::Logger::misc() << "End of cycle, COMPUTED norm of residual " << std::scientific << std::setprecision(16) - << rnorm << "\n"; + << res_norm << "\n"; } } // outer while + return 0; } @@ -583,11 +700,6 @@ namespace ReSolve return 0; } - void LinSolverIterativeFGMRES::precV(vector_type* rhs, vector_type* x) - { - preconditioner_->apply(rhs, x); - } - void LinSolverIterativeFGMRES::setMemorySpace() { bool is_matrix_handler_cuda = matrix_handler_->getIsCudaEnabled(); diff --git a/resolve/LinSolverIterativeFGMRES.hpp b/resolve/LinSolverIterativeFGMRES.hpp index de9ebc8d..a35ca22a 100644 --- a/resolve/LinSolverIterativeFGMRES.hpp +++ b/resolve/LinSolverIterativeFGMRES.hpp @@ -91,7 +91,6 @@ namespace ReSolve int freeSolverData(); void setMemorySpace(); void initParamList(); - void precV(vector_type* rhs, vector_type* x); ///< Apply preconditioner memory::MemorySpace memspace_; diff --git a/resolve/LinSolverIterativeRandFGMRES.cpp b/resolve/LinSolverIterativeRandFGMRES.cpp index 6651f289..3fdd1316 100644 --- a/resolve/LinSolverIterativeRandFGMRES.cpp +++ b/resolve/LinSolverIterativeRandFGMRES.cpp @@ -126,6 +126,18 @@ namespace ReSolve /** * @brief Solve linear system A * x = rhs * + * Implements restarted randomized GMRES with optional flexible + * variant. A sketching operator is used to orthogonalize in a reduced space + * + * Flexible randomized GMRES allows the preconditioner to vary per + * iteration and uses right preconditioning. Standard randomized GMRES + * supports both left and right preconditioning. + * + * Left preconditioning solves M^{-1}Ax = M^{-1}b and checks convergence + * with ||M^{-1}(b - Ax)||. Right preconditioning solves AM^{-1}(Mx) = b + * and checks convergence with ||b - Ax||. Both report the true relative + * residual ||b - Ax||/||b||. + * * @param rhs - right hand side vector * @param x - solution vector * @return int - zero if successful, error code otherwise @@ -137,6 +149,14 @@ namespace ReSolve { using namespace constants; + // Flexible GMRES only supports right preconditioning + if (flexible_ && preconditioner_->getSide() == Preconditioner::Side::LEFT) + { + out::warning() << "Flexible GMRES does not support left preconditioning. " + << "Switching to right preconditioning.\n"; + preconditioner_->setSide(Preconditioner::Side::RIGHT); + } + // io::Logger::setVerbosity(io::Logger::EVERYTHING); int outer_flag = 1; @@ -147,15 +167,20 @@ namespace ReSolve int k; int k1; - real_type t; - real_type rnorm; - real_type bnorm; + real_type t = 0.0; + real_type res_norm = 0.0; // Residual norm used used for convergence + real_type rhs_norm = 0.0; // Right-hand side norm used for convergence + real_type true_res_norm = 0.0; // True (unpreconditioned) residual norm ||b - Ax|| for reporting + real_type true_rhs_norm = 0.0; // True (unpreconditioned) right-hand side norm ||b|| for reporting real_type tolrel; vector_type vec_v(n_); vector_type vec_z(n_); vector_type vec_s(k_rand_); - // V[0] = b-A*x_0 - // debug + + // Compute initial residual norm. + // V[0] = ||b - A*x0|| for right preconditioning + // V[0] = ||M^{-1}{b - A*x0}|| for left preconditioning + vec_Z_->setToZero(memspace_); vec_V_->setToZero(memspace_); @@ -163,6 +188,23 @@ namespace ReSolve matrix_handler_->matvec(A_, x, vec_V_, &MINUS_ONE, &ONE, memspace_); vec_v.setData(vec_V_->getData(0, memspace_), memspace_); + vec_z.setData(vec_Z_->getData(0, memspace_), memspace_); + + // True residual norm ||b - A*x0|| + true_res_norm = vector_handler_->dot(&vec_v, &vec_v, memspace_); + true_res_norm = std::sqrt(true_res_norm); + + // True right-hand side norm ||b|| + true_rhs_norm = vector_handler_->dot(rhs, rhs, memspace_); + true_rhs_norm = std::sqrt(true_rhs_norm); + + // Left preconditioning uses preconditioned norms for convergence + if (preconditioner_->getSide() == Preconditioner::LEFT) + { + preconditioner_->apply(&vec_v, &vec_z); + vec_v.copyFromExternal(&vec_z, memspace_, memspace_); + } + vec_s.setData(vec_S_->getData(0, memspace_), memspace_); sketching_handler_->Theta(&vec_v, &vec_s); @@ -172,22 +214,33 @@ namespace ReSolve vector_handler_->scal(one_over_k_, &vec_s, memspace_); } mem_.deviceSynchronize(); + res_norm = vector_handler_->dot(&vec_s, &vec_s, memspace_); + rhs_norm = vector_handler_->dot(rhs, rhs, memspace_); + + // Left-preconditioned right-hand side norm ||M^{-1}*b|| + if (!flexible_ && preconditioner_->getSide() == Preconditioner::LEFT) + { + vec_v.setData(rhs->getData(memspace_), memspace_); + preconditioner_->apply(&vec_v, &vec_z); + rhs_norm = vector_handler_->dot(&vec_z, &vec_z, memspace_); + vec_v.setData(vec_V_->getData(0, memspace_), memspace_); + } + + res_norm = std::sqrt(res_norm); // res_norm = ||S_0|| + rhs_norm = std::sqrt(rhs_norm); - rnorm = 0.0; - bnorm = vector_handler_->dot(rhs, rhs, memspace_); - rnorm = vector_handler_->dot(&vec_s, &vec_s, memspace_); - rnorm = std::sqrt(rnorm); // rnorm = ||V_1|| - bnorm = std::sqrt(bnorm); io::Logger::misc() << "it 0: norm of residual " << std::scientific << std::setprecision(16) - << rnorm << " Norm of rhs: " << bnorm << "\n"; + << res_norm << " Norm of rhs: " << rhs_norm << "\n"; + + // Report the true initial relative residual norm + initial_residual_norm_ = true_res_norm / true_rhs_norm; - initial_residual_norm_ = rnorm / bnorm; // compute relative residual norm while (outer_flag) { if (it == 0) { - tolrel = tol_ * rnorm; + tolrel = tol_ * res_norm; if (std::abs(tolrel) < MACHINE_EPSILON) { tolrel = MACHINE_EPSILON; @@ -198,34 +251,34 @@ namespace ReSolve switch (conv_cond_) { case 0: - exit_cond = ((std::abs(rnorm - ZERO) <= MACHINE_EPSILON)); + exit_cond = ((std::abs(res_norm - ZERO) <= MACHINE_EPSILON)); break; case 1: - exit_cond = ((std::abs(rnorm - ZERO) <= MACHINE_EPSILON) || (rnorm < tol_)); + exit_cond = ((std::abs(res_norm - ZERO) <= MACHINE_EPSILON) || (res_norm < tol_)); break; case 2: - exit_cond = ((std::abs(rnorm - ZERO) <= MACHINE_EPSILON) || (rnorm < (tol_ * bnorm))); + exit_cond = ((std::abs(res_norm - ZERO) <= MACHINE_EPSILON) || (res_norm < (tol_ * rhs_norm))); break; } if (exit_cond) { outer_flag = 0; - final_residual_norm_ = rnorm; - initial_residual_norm_ = rnorm; + final_residual_norm_ = res_norm; + initial_residual_norm_ = res_norm; total_iters_ = 0; break; } // normalize first vector - t = 1.0 / rnorm; + t = 1.0 / res_norm; vector_handler_->scal(t, vec_V_, memspace_); vector_handler_->scal(t, vec_S_, memspace_); mem_.deviceSynchronize(); // initialize norm history - h_rs_[0] = rnorm; + h_rs_[0] = res_norm; i = -1; notconv = 1; @@ -234,7 +287,6 @@ namespace ReSolve i++; it++; - // Z_i = (LU)^{-1}*V_i vec_v.setData(vec_V_->getData(i, memspace_), memspace_); if (flexible_) { @@ -244,14 +296,35 @@ namespace ReSolve { vec_z.setData(vec_Z_->getData(0, memspace_), memspace_); } - this->precV(&vec_v, &vec_z); - - mem_.deviceSynchronize(); - // V_{i+1}=A*Z_i - vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_); + // Expand the Krylov subspace. + // + // New basis vector: + // V[i+1] = A*M^{-1}*V[i] (right preconditioning) + // V[i+1] = M^{-1}*A*V[i] (left preconditioning) - matrix_handler_->matvec(A_, &vec_z, &vec_v, &ONE, &ZERO, memspace_); + switch (preconditioner_->getSide()) + { + case Preconditioner::Side::RIGHT: + // Compute vec_z = M^{-1}*V[i], then V[i+1] = A*vec_z + preconditioner_->apply(&vec_v, &vec_z); + mem_.deviceSynchronize(); + + vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_); + matrix_handler_->matvec(A_, &vec_z, &vec_v, &ONE, &ZERO, memspace_); + break; + case Preconditioner::Side::LEFT: + // Compute vec_z = A*V[i], then V[i+1] = M^{-1}*vec_z + matrix_handler_->matvec(A_, &vec_v, &vec_z, &ONE, &ZERO, memspace_); + + vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_); + preconditioner_->apply(&vec_z, &vec_v); + mem_.deviceSynchronize(); + break; + default: + out::error() << "Unknown preconditioner side.\n"; + return 1; + } // orthogonalize V[i+1], form a column of h_H_ // this is where it differs from normal solver GS @@ -305,13 +378,13 @@ namespace ReSolve h_H_[(i) * (restart_ + 1) + (i + 1)] = h_c_[i] * Hii1 - h_s_[i] * Hii; // residual norm estimate - rnorm = std::abs(h_rs_[i + 1]); + res_norm = std::abs(h_rs_[i + 1]); io::Logger::misc() << "it: " << it << " --> norm of the residual " << std::scientific << std::setprecision(16) - << rnorm << "\n"; + << res_norm << "\n"; // check convergence - if (i + 1 >= restart_ || rnorm <= tolrel || it >= maxit_) + if (i + 1 >= restart_ || res_norm <= tolrel || it >= maxit_) { notconv = 0; } @@ -319,7 +392,7 @@ namespace ReSolve io::Logger::misc() << "End of cycle, ESTIMATED norm of residual " << std::scientific << std::setprecision(16) - << rnorm << "\n"; + << res_norm << "\n"; // solve tri system h_rs_[i] = h_rs_[i] / h_H_[i * (restart_ + 1) + i]; for (int ii = 2; ii <= i + 1; ii++) @@ -334,7 +407,11 @@ namespace ReSolve h_rs_[k] = t / h_H_[k * (restart_ + 1) + k]; } - // get solution + // Update the approximate solution x using h_rs_. + // Flexible GMRES uses the preconditioned basis Z[j] directly. + // Standard GMRES first forms vec_z from V[j], then applies M^{-1} + // only for right preconditioning. + if (flexible_) { for (j = 0; j <= i; j++) @@ -345,23 +422,34 @@ namespace ReSolve } else { - vec_Z_->setToZero(0, memspace_); + // Accumulate the correction vec_z = sum_j h_rs_[j] * V[j] + vec_Z_->setToZero(memspace_); vec_z.setData(vec_Z_->getData(0, memspace_), memspace_); for (j = 0; j <= i; j++) { vec_v.setData(vec_V_->getData(j, memspace_), memspace_); vector_handler_->axpy(h_rs_[j], &vec_v, &vec_z, memspace_); } - // now multiply d_Z by precon - - vec_v.setData(vec_V_->getData(memspace_), memspace_); - this->precV(&vec_z, &vec_v); - // and add to x - vector_handler_->axpy(ONE, &vec_v, x, memspace_); + // Apply the correction to x based on preconditioning side + switch (preconditioner_->getSide()) + { + case Preconditioner::Side::RIGHT: + // Right preconditioning: x += M^{-1} * vec_z + preconditioner_->apply(&vec_z, &vec_v); + vector_handler_->axpy(ONE, &vec_v, x, memspace_); + break; + case Preconditioner::Side::LEFT: + // Left preconditioning: x += vec_z + vector_handler_->axpy(ONE, &vec_z, x, memspace_); + break; + default: + out::error() << "Unknown preconditioner side.\n"; + return 1; + } } /* test solution */ - if (rnorm <= tolrel || it >= maxit_) + if (res_norm <= tolrel || it >= maxit_) { outer_flag = 0; } @@ -370,6 +458,13 @@ namespace ReSolve matrix_handler_->matvec(A_, x, vec_V_, &MINUS_ONE, &ONE, memspace_); if (outer_flag) { + // Left-preconditioned GMRES applies M^{-1} to the residual before sketching + if (!flexible_ && preconditioner_->getSide() == Preconditioner::LEFT) + { + vec_v.setData(vec_V_->getData(0, memspace_), memspace_); + preconditioner_->apply(&vec_v, &vec_z); + vec_v.copyFromExternal(&vec_z, memspace_, memspace_); + } sketching_handler_->reset(); @@ -385,22 +480,23 @@ namespace ReSolve vector_handler_->scal(one_over_k_, &vec_s, memspace_); } mem_.deviceSynchronize(); - rnorm = vector_handler_->dot(vec_S_, vec_S_, memspace_); - // rnorm = ||S_0|| - rnorm = std::sqrt(rnorm); + res_norm = vector_handler_->dot(vec_S_, vec_S_, memspace_); + // res_norm = ||S_0|| + res_norm = std::sqrt(res_norm); } if (!outer_flag) { - rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_); - // rnorm = ||V_0|| - rnorm = std::sqrt(rnorm); + // Compute the true residual norm = ||b - A*x|| + true_res_norm = vector_handler_->dot(vec_V_, vec_V_, memspace_); + true_res_norm = std::sqrt(true_res_norm); io::Logger::misc() << "End of cycle, COMPUTED norm of residual " << std::scientific << std::setprecision(16) - << rnorm << "\n"; + << true_res_norm << "\n"; - final_residual_norm_ = rnorm / bnorm; // relative residual norm + // Report the true relative residual norm + final_residual_norm_ = true_res_norm / true_rhs_norm; total_iters_ = it; } } // outer while @@ -773,11 +869,6 @@ namespace ReSolve return 0; } - void LinSolverIterativeRandFGMRES::precV(vector_type* rhs, vector_type* x) - { - preconditioner_->apply(rhs, x); - } - /** * @brief Set memory space and device tape based on how MatrixHandler * and VectorHandler are configured. diff --git a/resolve/LinSolverIterativeRandFGMRES.hpp b/resolve/LinSolverIterativeRandFGMRES.hpp index e2da24b8..2098201b 100644 --- a/resolve/LinSolverIterativeRandFGMRES.hpp +++ b/resolve/LinSolverIterativeRandFGMRES.hpp @@ -108,7 +108,6 @@ namespace ReSolve int freeSketchingData(); void setMemorySpace(); void initParamList(); - void precV(vector_type* rhs, vector_type* x); ///< Apply preconditioner memory::MemorySpace memspace_; diff --git a/resolve/Preconditioner.cpp b/resolve/Preconditioner.cpp index 4d85556d..6f3886fc 100644 --- a/resolve/Preconditioner.cpp +++ b/resolve/Preconditioner.cpp @@ -22,4 +22,21 @@ namespace ReSolve return 1; } + Preconditioner::Side Preconditioner::getSide() const + { + return side_; + } + + /** + * @brief Set the preconditioning side + * + * @param[in] side - side of preconditioning + * @return 0 if successful + */ + int Preconditioner::setSide(Side side) + { + side_ = side; + return 0; + } + } // namespace ReSolve diff --git a/resolve/Preconditioner.hpp b/resolve/Preconditioner.hpp index 4ed21cda..e484c97f 100644 --- a/resolve/Preconditioner.hpp +++ b/resolve/Preconditioner.hpp @@ -29,11 +29,23 @@ namespace ReSolve using vector_type = vector::Vector; using matrix_type = matrix::Sparse; + enum Side + { + LEFT = 0, + RIGHT + }; + Preconditioner(); virtual ~Preconditioner(); - virtual int setup(matrix_type* A) = 0; virtual int apply(vector_type* rhs, vector_type* x) = 0; + virtual int setup(matrix_type* A) = 0; virtual int reset(matrix_type* /* A */); + + Side getSide() const; // Gets the preconditioning side + int setSide(Side side); // Sets the preconditioning side + + private: + Side side_{Side::RIGHT}; // Right preconditioning by default }; } // namespace ReSolve diff --git a/resolve/PreconditionerLU.cpp b/resolve/PreconditionerLU.cpp index fd6ec3c7..0c1d2960 100644 --- a/resolve/PreconditionerLU.cpp +++ b/resolve/PreconditionerLU.cpp @@ -1,7 +1,7 @@ /** * @file PreconditionerLU.cpp * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) - * @brief Declaration of preconditioner ILU0 class. + * @brief Implementation of preconditioner LU class. * */ @@ -29,7 +29,7 @@ namespace ReSolve } /** - * @brief Sets up the preconditioner with the given matrix + * @brief Sets up the lu solver with the given matrix * * @param[in] A - System matrix to set up the preconditioner with * diff --git a/resolve/PreconditionerLU.hpp b/resolve/PreconditionerLU.hpp index 14ede135..136a5e92 100644 --- a/resolve/PreconditionerLU.hpp +++ b/resolve/PreconditionerLU.hpp @@ -1,7 +1,7 @@ /** - * @file PreconditionerILU0.cpp + * @file PreconditionerLU.hpp * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) - * @brief Declaration of preconditioner ILU0 class. + * @brief Declaration of preconditioner LU class. * */ diff --git a/resolve/PreconditionerUserMatrix.cpp b/resolve/PreconditionerUserMatrix.cpp new file mode 100644 index 00000000..4c309f20 --- /dev/null +++ b/resolve/PreconditionerUserMatrix.cpp @@ -0,0 +1,97 @@ +/** + * @file PreconditionerUserMatrix.cpp + * @author Jeffery Zhang (jefferyz@vt.edu) + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Implementation of preconditioner class with user matrix + * + */ + +#include "PreconditionerUserMatrix.hpp" + +namespace ReSolve +{ + /** + * @brief Constructor for PreconditionerUserMatrix. + * + * @param[in] matrix_handler - Pointer to the matrix handler + */ + PreconditionerUserMatrix::PreconditionerUserMatrix(MatrixHandler* matrix_handler) + { + matrix_handler_ = matrix_handler; + setMemorySpace(); + } + + /** + * @brief Destructor for PreconditionerUserMatrix + */ + PreconditionerUserMatrix::~PreconditionerUserMatrix() + { + } + + /** + * @brief This preconditioner does not depend on A. + */ + int PreconditionerUserMatrix::setup(matrix_type* /* A */) + { + return 0; + } + + /** + * @brief Set the explicit preconditioner matrix B. + * + * @param[in] B - Pointer to the preconditioning matrix + */ + int PreconditionerUserMatrix::setPrecMatrix(matrix_type* B) + { + B_ = B; + return 0; + } + + /** + * @brief Get the preconditioner matrix + * + * Necessary for some calculations + * + */ + matrix::Sparse* PreconditionerUserMatrix::getPrecMatrix() + { + return B_; + } + + /** + * @brief Applies the preconditioner depending on the side + * + * @param[in] rhs - Right-hand-side vector + * @param[in] x - Solution vector + * + * @return int 0 if successful, 1 if fails + */ + int PreconditionerUserMatrix::apply(vector_type* rhs, vector_type* x) + { + using namespace constants; + + if (matrix_handler_ == nullptr || B_ == nullptr) + { + return 1; + } + + matrix_handler_->matvec(B_, rhs, x, &ONE, &ZERO, memspace_); + return 0; + } + + void PreconditionerUserMatrix::setMemorySpace() + { + bool is_matrix_handler_cuda = matrix_handler_->getIsCudaEnabled(); + bool is_matrix_handler_hip = matrix_handler_->getIsHipEnabled(); + + if (is_matrix_handler_cuda || is_matrix_handler_hip) + { + memspace_ = memory::DEVICE; + } + else + { + memspace_ = memory::HOST; + } + } + +} // namespace ReSolve diff --git a/resolve/PreconditionerUserMatrix.hpp b/resolve/PreconditionerUserMatrix.hpp new file mode 100644 index 00000000..edf5bc95 --- /dev/null +++ b/resolve/PreconditionerUserMatrix.hpp @@ -0,0 +1,60 @@ +/** + * @file PreconditionerUserMatrix.hpp + * @author Jeffery Zhang (jefferyz@vt.edu) + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Declaration of preconditioner class with user matrix + * + */ + +#pragma once + +#include + +#include "Common.hpp" +#include +#include +#include + +namespace ReSolve +{ + + namespace matrix + { + class Sparse; + } // namespace matrix + + namespace vector + { + class Vector; + } // namespace vector + + // Forward declaration of MatrixHandler class + class MatrixHandler; + + /** + * @brief Class allows for a user specified matrix to be used for preconditioning + * + */ + class PreconditionerUserMatrix : public Preconditioner + { + public: + using vector_type = vector::Vector; + using matrix_type = matrix::Sparse; + + PreconditionerUserMatrix(MatrixHandler* matrix_handler); + ~PreconditionerUserMatrix(); + + int setup(matrix_type* A) override; + int setPrecMatrix(matrix_type* B); + int apply(vector_type* rhs, vector_type* x) override; + + matrix_type* getPrecMatrix(); + + private: + void setMemorySpace(); + + matrix_type* B_{nullptr}; + MatrixHandler* matrix_handler_{nullptr}; + memory::MemorySpace memspace_; + }; +} // namespace ReSolve diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index f7a4acf8..6b4cbc0c 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -562,18 +562,46 @@ namespace ReSolve * * @return int 0 if successful, 1 if it fails */ - int SystemSolver::preconditionerSetup() + int SystemSolver::preconditionerSetup(std::string side) { int status = 0; - if (precondition_method_ == "ilu0") + + if (preconditioner_ == nullptr) { - status += preconditioner_->setup(A_); - if (memspace_ != "cpu") - { - is_solve_on_device_ = true; - } - status += iterativeSolver_->setPreconditioner(preconditioner_); + out::error() << "Preconditioner not initialized!\n"; + status += 1; + } + + if (iterativeSolver_ == nullptr) + { + out::error() << "Iterative solver not initialized!\n"; + status += 1; + } + + Preconditioner::Side prec_side; + if (side == "left") + { + prec_side = Preconditioner::LEFT; + } + else if (side == "right") + { + prec_side = Preconditioner::RIGHT; } + else + { + out::warning() << "Preconditioning side " << " not recognized.\n"; + out::warning() << "Using default preconditioning side (right).\n"; + prec_side = Preconditioner::RIGHT; + } + + status += preconditioner_->setSide(prec_side); + status += preconditioner_->setup(A_); + + if (memspace_ != "cpu") + { + is_solve_on_device_ = true; + } + status += iterativeSolver_->setPreconditioner(preconditioner_); return status; } @@ -591,11 +619,15 @@ namespace ReSolve { int status = 0; A_ = A; - if (precondition_method_ == "ilu0") + + if (preconditioner_ == nullptr) { - status += preconditioner_->reset(A); + out::error() << "Preconditioner not initialized!\n"; + status += 1; } + status += preconditioner_->reset(A); + return status; } @@ -623,6 +655,11 @@ namespace ReSolve return *iterativeSolver_; } + Preconditioner& SystemSolver::getPreconditioner() + { + return *preconditioner_; + } + void SystemSolver::setFactorizationMethod(std::string method) { factorizationMethod_ = method; diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp index ca0c313d..67a3307d 100644 --- a/resolve/SystemSolver.hpp +++ b/resolve/SystemSolver.hpp @@ -1,3 +1,10 @@ +#pragma once + +#include + +#include +#include + // this is to solve the system, can call different linear solvers if necessary namespace ReSolve { @@ -10,7 +17,6 @@ namespace ReSolve class LinAlgWorkspaceCpu; class MatrixHandler; class VectorHandler; - class Preconditioner; namespace vector { @@ -55,7 +61,7 @@ namespace ReSolve int factorize(); // numeric part int refactorize(); int refactorizationSetup(); - int preconditionerSetup(); + int preconditionerSetup(std::string side); int resetPreconditioner(matrix_type* A); int solve(vector_type* rhs, vector_type* x); // for direct and iterative int refine(vector_type* rhs, vector_type* x); // for iterative refinement @@ -66,6 +72,7 @@ namespace ReSolve LinSolverDirect& getFactorizationSolver(); LinSolverDirect& getRefactorizationSolver(); LinSolverIterative& getIterativeSolver(); + Preconditioner& getPreconditioner(); real_type getVectorNorm(vector_type* rhs); real_type getResidualNorm(vector_type* rhs, vector_type* x); diff --git a/tests/functionality/CMakeLists.txt b/tests/functionality/CMakeLists.txt index ab8c2413..c36e18a2 100644 --- a/tests/functionality/CMakeLists.txt +++ b/tests/functionality/CMakeLists.txt @@ -183,6 +183,56 @@ add_test(NAME sys_gmres_mgspm_test "-g" "mgs_pm" ) +# Left-preconditioned Krylov solvers tests (GMRES) +add_test(NAME sys_rand_count_gmres_left_cgs2_test + COMMAND $ "-x" "no" "-i" + "randgmres" "-g" "cgs2" "-s" "count" "-p" "left" +) +add_test(NAME sys_rand_count_gmres_left_mgs_test + COMMAND $ "-x" "no" "-i" + "randgmres" "-g" "mgs" "-s" "count" "-p" "left" +) +add_test(NAME sys_rand_count_gmres_left_mgs2sync_test + COMMAND $ "-x" "no" "-i" + "randgmres" "-g" "mgs_two_sync" "-s" "count" "-p" "left" +) +add_test(NAME sys_rand_count_gmres_left_mgspm_test + COMMAND $ "-x" "no" "-i" + "randgmres" "-g" "mgs_pm" "-s" "count" "-p" "left" +) +add_test(NAME sys_rand_fwht_gmres_left_cgs2_test + COMMAND $ "-x" "no" "-i" + "randgmres" "-g" "cgs2" "-s" "fwht" "-p" "left" +) +add_test(NAME sys_rand_fwht_gmres_left_mgs_test + COMMAND $ "-x" "no" "-i" + "randgmres" "-g" "mgs" "-s" "fwht" "-p" "left" +) +add_test(NAME sys_rand_fwht_gmres_left_mgs2sync_test + COMMAND $ "-x" "no" "-i" + "randgmres" "-g" "mgs_two_sync" "-s" "fwht" "-p" "left" +) +add_test(NAME sys_rand_fwht_gmres_left_mgspm_test + COMMAND $ "-x" "no" "-i" + "randgmres" "-g" "mgs_pm" "-s" "fwht" "-p" "left" +) +add_test(NAME sys_gmres_left_cgs2_test + COMMAND $ "-x" "no" "-i" "fgmres" + "-g" "cgs2" "-p" "left" +) +add_test(NAME sys_gmres_left_mgs_test + COMMAND $ "-x" "no" "-i" "fgmres" + "-g" "mgs" "-p" "left" +) +add_test(NAME sys_gmres_left_mgs2sync_test + COMMAND $ "-x" "no" "-i" "fgmres" + "-g" "mgs_two_sync" "-p" "left" +) +add_test(NAME sys_gmres_left_mgspm_test + COMMAND $ "-x" "no" "-i" "fgmres" + "-g" "mgs_pm" "-p" "left" +) + # Randomized GMRES solvers tested without SystemSolver class add_test(NAME rand_gmres_test COMMAND $) diff --git a/tests/functionality/testKlu.cpp b/tests/functionality/testKlu.cpp index af537357..4df39afd 100644 --- a/tests/functionality/testKlu.cpp +++ b/tests/functionality/testKlu.cpp @@ -151,8 +151,8 @@ int runTest(int argc, char* argv[], std::string& solver_name) status = FGMRES.setup(A); error_sum += status; - ReSolve::PreconditionerLU precond_lu(&KLU); - status = FGMRES.setPreconditioner(&precond_lu); + ReSolve::PreconditionerLU precond(&KLU); + status = FGMRES.setPreconditioner(&precond); error_sum += status; status = FGMRES.solve(&vec_rhs, &vec_x); error_sum += status; @@ -199,8 +199,8 @@ int runTest(int argc, char* argv[], std::string& solver_name) if (is_ir) { FGMRES.resetMatrix(A); - ReSolve::PreconditionerLU precond_lu(&KLU); - status = FGMRES.setPreconditioner(&precond_lu); + ReSolve::PreconditionerLU precond(&KLU); + status = FGMRES.setPreconditioner(&precond); error_sum += status; status = FGMRES.solve(&vec_rhs, &vec_x); error_sum += status; diff --git a/tests/functionality/testRandGmres.cpp b/tests/functionality/testRandGmres.cpp index ca9adb46..e667bfa6 100644 --- a/tests/functionality/testRandGmres.cpp +++ b/tests/functionality/testRandGmres.cpp @@ -133,8 +133,8 @@ int runTest(int argc, char* argv[]) FGMRES.setRestart(200); FGMRES.setSketchingMethod(LinSolverIterativeRandFGMRES::cs); - PreconditionerLU precond_lu(&ILU); - status = FGMRES.setPreconditioner(&precond_lu); + PreconditionerLU precond(&ILU); + status = FGMRES.setPreconditioner(&precond); error_sum += status; FGMRES.setFlexible(true); diff --git a/tests/functionality/testRefactor.cpp b/tests/functionality/testRefactor.cpp index a765ea7b..7d72a401 100644 --- a/tests/functionality/testRefactor.cpp +++ b/tests/functionality/testRefactor.cpp @@ -191,8 +191,8 @@ int runTest(int argc, char* argv[], std::string& solver_name) status = FGMRES.setup(A); error_sum += status; - ReSolve::PreconditionerLU precond_lu(&Rf); - status = FGMRES.setPreconditioner(&precond_lu); + ReSolve::PreconditionerLU precond(&Rf); + status = FGMRES.setPreconditioner(&precond); error_sum += status; status = FGMRES.solve(&vec_rhs, &vec_x); @@ -241,8 +241,8 @@ int runTest(int argc, char* argv[], std::string& solver_name) if (is_ir) { FGMRES.resetMatrix(A); - ReSolve::PreconditionerLU precond_lu(&Rf); - status = FGMRES.setPreconditioner(&precond_lu); + ReSolve::PreconditionerLU precond(&Rf); + status = FGMRES.setPreconditioner(&precond); error_sum += status; status = FGMRES.solve(&vec_rhs, &vec_x); diff --git a/tests/functionality/testSysGmres.cpp b/tests/functionality/testSysGmres.cpp index 277c1842..a7b5d162 100644 --- a/tests/functionality/testSysGmres.cpp +++ b/tests/functionality/testSysGmres.cpp @@ -16,6 +16,7 @@ #include "TestHelper.hpp" #include #include +#include #include #include #include @@ -41,7 +42,7 @@ template static int test(int argc, char* argv[]); /// Checks if inputs are valid, otherwise sets defaults -static void processInputs(std::string& method, std::string& gs, std::string& sketch); +static void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& side); /// Creates string with test description static std::string headerInfo(const std::string& method, @@ -102,7 +103,10 @@ int test(int argc, char* argv[]) opt = options.getParamFromKey("-x"); std::string flexible = opt ? (*opt).second : "yes"; - processInputs(method, gs, sketch); + opt = options.getParamFromKey("-p"); + std::string side = opt ? (*opt).second : "right"; + + processInputs(method, gs, sketch, side); // Create workspace and initialize its handles. workspace_type workspace; @@ -166,7 +170,7 @@ int test(int argc, char* argv[]) solver.getIterativeSolver().setCliParam("restart", "200"); // Set preconditioner (default in this case ILU0) - status = solver.preconditionerSetup(); + status = solver.preconditionerSetup(side); error_sum += status; // Solve system @@ -196,7 +200,7 @@ int test(int argc, char* argv[]) // Definitions of helper functions // -void processInputs(std::string& method, std::string& gs, std::string& sketch) +void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& side) { if (method == "randgmres") { @@ -214,6 +218,7 @@ void processInputs(std::string& method, std::string& gs, std::string& sketch) std::cout << "Setting iterative solver method to the default (FGMRES).\n\n"; method = "fgmres"; } + if (gs != "cgs1" && gs != "cgs2" && gs != "mgs" && gs != "mgs_two_sync" && gs != "mgs_pm") { @@ -221,6 +226,13 @@ void processInputs(std::string& method, std::string& gs, std::string& sketch) std::cout << "Setting orthogonalization to the default (CGS2).\n\n"; gs = "cgs2"; } + + if ((side != "left") && (side != "right")) + { + std::cout << "Preconditioning side " << side << " not recognized.\n"; + std::cout << "Setting preconditioning side to the default (right).\n\n"; + side = "right"; + } } std::string headerInfo(const std::string& method, diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 1f4df955..a2526457 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -8,6 +8,7 @@ add_subdirectory(matrix) add_subdirectory(vector) +add_subdirectory(preconditioner) add_subdirectory(utilities) add_subdirectory(memory) add_subdirectory(params) diff --git a/tests/unit/preconditioner/CMakeLists.txt b/tests/unit/preconditioner/CMakeLists.txt new file mode 100644 index 00000000..63df670e --- /dev/null +++ b/tests/unit/preconditioner/CMakeLists.txt @@ -0,0 +1,27 @@ +#[[ + +@brief Build ReSolve preconditioner unit tests + +]] + +add_executable( + runPreconditionerUserMatrixTests.exe runPreconditionerUserMatrixTests.cpp +) +add_executable(runPreconditionerLUTests.exe runPreconditionerLUTests.cpp) + +target_link_libraries(runPreconditionerUserMatrixTests.exe PRIVATE ReSolve) +target_link_libraries(runPreconditionerLUTests.exe PRIVATE ReSolve) + +install(TARGETS runPreconditionerUserMatrixTests.exe + RUNTIME DESTINATION bin/resolve/tests/unit +) +install(TARGETS runPreconditionerLUTests.exe + RUNTIME DESTINATION bin/resolve/tests/unit +) + +add_test(NAME preconditioner_user_matrix_test + COMMAND $ +) +add_test(NAME preconditioner_lu_test + COMMAND $ +) diff --git a/tests/unit/preconditioner/PreconditionerLUTests.hpp b/tests/unit/preconditioner/PreconditionerLUTests.hpp new file mode 100644 index 00000000..7f415e81 --- /dev/null +++ b/tests/unit/preconditioner/PreconditionerLUTests.hpp @@ -0,0 +1,130 @@ +/** + * @file PreconditionerLUTests.hpp + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Tests for PreconditionerILU class + * + */ + +#pragma once + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace ReSolve +{ + namespace tests + { + /** + * @brief Unit tests for PreconditionerLU. + * + */ + class PreconditionerLUTests : public TestBase + { + public: + PreconditionerLUTests(memory::MemorySpace memspace, LinSolverDirect* lu_solver) + { + memspace_ = memspace; + lu_solver_ = lu_solver; + } + + virtual ~PreconditionerLUTests() + { + } + + /** + * @brief Solve a 3x3 diagonal system through PreconditionerLU + * + */ + TestOutcome solve() + { + TestStatus status; + std::string testname(__func__); + + const index_type n = 3; + const index_type nnz = 3; + + matrix::Csr* A = new matrix::Csr(n, n, nnz); + index_type row_data[4] = {0, 1, 2, 3}; + index_type col_data[3] = {0, 1, 2}; + real_type val_data[3] = {4.0, 5.0, 6.0}; + + A->copyFromExternal(row_data, col_data, val_data, memory::HOST, memory::HOST); + + if (memspace_ == memory::DEVICE) + { + A->syncData(memory::DEVICE); + } + + PreconditionerLU precond(lu_solver_); + precond.setup(A); + + real_type rhs_data[3] = {4.0, 10.0, 18.0}; + vector::Vector* rhs = new vector::Vector(n); + rhs->copyFromExternal(rhs_data, memory::HOST, memspace_); + + vector::Vector* x = new vector::Vector(n); + x->allocate(memspace_); + + precond.apply(rhs, x); + + if (memspace_ == memory::DEVICE) + { + x->syncData(memory::HOST); + } + + const real_type expected[3] = {1.0, 2.0, 3.0}; + const real_type tol = 1e-12; + for (index_type i = 0; i < n; ++i) + { + if (std::fabs(x->getData(memory::HOST)[i] - expected[i]) > tol) + { + status *= false; + std::cout << testname << ": x[" << i << "] = " << x->getData(memory::HOST)[i] + << ", expected " << expected[i] << "\n"; + } + } + + delete A; + delete rhs; + delete x; + + return status.report(testname.c_str()); + } + + TestOutcome checkSide() + { + TestStatus status; + std::string testname(__func__); + + PreconditionerLU precond(lu_solver_); + + if (precond.getSide() != Preconditioner::Side::RIGHT) + { + status *= false; + std::cout << "Default side should be Side::RIGHT\n"; + } + + precond.setSide(Preconditioner::Side::LEFT); + if (precond.getSide() != Preconditioner::Side::LEFT) + { + status *= false; + std::cout << "After setSide(Side::LEFT), side was not updated\n"; + } + + return status.report(testname.c_str()); + } + + private: + memory::MemorySpace memspace_; + LinSolverDirect* lu_solver_; + }; + + } // namespace tests +} // namespace ReSolve diff --git a/tests/unit/preconditioner/PreconditionerUserMatrixTests.hpp b/tests/unit/preconditioner/PreconditionerUserMatrixTests.hpp new file mode 100644 index 00000000..e79267a6 --- /dev/null +++ b/tests/unit/preconditioner/PreconditionerUserMatrixTests.hpp @@ -0,0 +1,182 @@ +/** + * @file PreconditionerUserMatrixTests.hpp + * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp) + * @brief Tests for PreconditionerUserMatrix class. + * + */ + +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace ReSolve +{ + namespace tests + { + /** + * @class Unit tests for PreconditionerUserMatrix. + * + */ + class PreconditionerUserMatrixTests : public TestBase + { + public: + PreconditionerUserMatrixTests(memory::MemorySpace memspace, MatrixHandler& handler) + : memspace_(memspace), + handler_(handler) + { + } + + virtual ~PreconditionerUserMatrixTests() + { + } + + /** + * @brief Test default side is right and setSide works for both enum values. + * + */ + TestOutcome checkSide() + { + TestStatus status; + std::string testname(__func__); + + PreconditionerUserMatrix precond(&handler_); + + if (precond.getSide() != Preconditioner::Side::RIGHT) + { + status *= false; + std::cout << testname << ": default side is not Side::RIGHT\n"; + } + + if (precond.setSide(Preconditioner::Side::LEFT) != 0 + || precond.getSide() != Preconditioner::Side::LEFT) + { + status *= false; + std::cout << testname << ": setSide(Side::LEFT) failed\n"; + } + + if (precond.setSide(Preconditioner::Side::RIGHT) != 0 + || precond.getSide() != Preconditioner::Side::RIGHT) + { + status *= false; + std::cout << testname << ": setSide(Side::RIGHT) failed\n"; + } + + return status.report(testname.c_str()); + } + + /** + * @brief Test setPrecMatrix(B) stores B and getPrecMatrix() retrieves it. + */ + TestOutcome setPrecMatrix() + { + TestStatus status; + status = true; + + PreconditionerUserMatrix precond(&handler_); + matrix::Csr B(3, 3, 3); + + int ret = precond.setPrecMatrix(&B); + if (ret != 0) + { + std::cout << "setMatrix(B) returned " << ret << ", expected 0\n"; + status *= false; + } + if (precond.getPrecMatrix() != &B) + { + std::cout << "getPrecMatrix() did not return the matrix passed to setPrecMatrix()\n"; + status *= false; + } + + return status.report(__func__); + } + + /** + * @brief Test apply() performs correct matvec with B. + * + * B = diag(2, 3, 4), x = [1, 1, 1] => y = [2, 3, 4]. + */ + TestOutcome apply() + { + TestStatus status; + status = true; + + const index_type n = 3; + const index_type nnz = 3; + matrix::Csr* B = new matrix::Csr(n, n, nnz); + + index_type row_data[4] = {0, 1, 2, 3}; + index_type col_data[3] = {0, 1, 2}; + real_type val_data[3] = {2.0, 3.0, 4.0}; + + B->copyFromExternal(row_data, col_data, val_data, memory::HOST, memory::HOST); + + if (memspace_ == memory::DEVICE) + { + B->allocateMatrixData(memspace_); + B->syncData(memspace_); + } + + vector::Vector x(n); + x.allocate(memory::HOST); + x.allocate(memspace_); + + for (index_type i = 0; i < n; ++i) + { + x.getData(memory::HOST)[i] = 1.0; + } + + x.setDataUpdated(memory::HOST); + + if (memspace_ == memory::DEVICE) + { + x.syncData(memspace_); + } + + vector::Vector y(n); + y.allocate(memory::HOST); + y.allocate(memspace_); + y.setToZero(memspace_); + + PreconditionerUserMatrix precond(&handler_); + precond.setPrecMatrix(B); + + int ret = precond.apply(&x, &y); + + if (ret != 0) + { + std::cout << "apply() returned " << ret << ", expected 0\n"; + status *= false; + } + + if (memspace_ == memory::DEVICE) + { + y.syncData(memory::HOST); + } + + const real_type expected[3] = {2.0, 3.0, 4.0}; + for (index_type i = 0; i < n; ++i) + { + if (!isEqual(y.getData(memory::HOST)[i], expected[i])) + { + std::cout << "y[" << i << "] = " << y.getData(memory::HOST)[i] + << ", expected " << expected[i] << "\n"; + status *= false; + } + } + + delete B; + return status.report(__func__); + } + + private: + memory::MemorySpace memspace_{memory::HOST}; + ReSolve::MatrixHandler& handler_; + }; + + } // namespace tests +} // namespace ReSolve diff --git a/tests/unit/preconditioner/runPreconditionerLUTests.cpp b/tests/unit/preconditioner/runPreconditionerLUTests.cpp new file mode 100644 index 00000000..f39f739c --- /dev/null +++ b/tests/unit/preconditioner/runPreconditionerLUTests.cpp @@ -0,0 +1,64 @@ +/** + * @file runPreconditionerLUTests.cpp + * @brief Tests for PreconditionerLU class. + * + */ + +#include +#include + +#include +#include + +#ifdef RESOLVE_USE_CUDA +#include +#endif + +#ifdef RESOLVE_USE_HIP +#include +#endif + +#include "PreconditionerLUTests.hpp" + +/** + * @brief Run PreconditionerLU tests with a given backend. + * + */ +template +void runTests(const std::string& backend, + ReSolve::memory::MemorySpace memspace, + ReSolve::tests::TestingResults& result) +{ + std::cout << "Running PreconditionerLU tests on " << backend << ":\n"; + + WorkspaceType workspace; + workspace.initializeHandles(); + ILUSolverType ilu_solver(&workspace); + + ReSolve::tests::PreconditionerLUTests test(memspace, &ilu_solver); + + result += test.checkSide(); + result += test.solve(); + + std::cout << "\n"; +} + +int main(int, char**) +{ + ReSolve::tests::TestingResults result; + + runTests("CPU", ReSolve::memory::HOST, result); + +#ifdef RESOLVE_USE_CUDA + runTests("CUDA", ReSolve::memory::DEVICE, result); +#endif + +#ifdef RESOLVE_USE_HIP + runTests("HIP", ReSolve::memory::DEVICE, result); +#endif + + return result.summary(); +} diff --git a/tests/unit/preconditioner/runPreconditionerUserMatrixTests.cpp b/tests/unit/preconditioner/runPreconditionerUserMatrixTests.cpp new file mode 100644 index 00000000..3d8f1c83 --- /dev/null +++ b/tests/unit/preconditioner/runPreconditionerUserMatrixTests.cpp @@ -0,0 +1,53 @@ +/** + * @file runPreconditionerUserMatrixTests.cpp + * @brief Tests for PreconditionerUserMatrix class. + * + */ + +#include +#include + +#include "PreconditionerUserMatrixTests.hpp" +#include +#include + +/** + * @brief Run PreconditionerUserMatrix tests with a given backend. + * + */ +template +void runExplicitMatrixTests(const std::string& backend, + ReSolve::memory::MemorySpace memspace, + ReSolve::tests::TestingResults& result) +{ + std::cout << "Running PreconditionerUserMatrix tests on " << backend << ":\n"; + + WorkspaceType workspace; + workspace.initializeHandles(); + ReSolve::MatrixHandler handler(&workspace); + + ReSolve::tests::PreconditionerUserMatrixTests test(memspace, handler); + + result += test.checkSide(); + result += test.setPrecMatrix(); + result += test.apply(); + + std::cout << "\n"; +} + +int main(int, char**) +{ + ReSolve::tests::TestingResults result; + + runExplicitMatrixTests("CPU", ReSolve::memory::HOST, result); + +#ifdef RESOLVE_USE_CUDA + runExplicitMatrixTests("CUDA", ReSolve::memory::DEVICE, result); +#endif + +#ifdef RESOLVE_USE_HIP + runExplicitMatrixTests("HIP", ReSolve::memory::DEVICE, result); +#endif + + return result.summary(); +}