From f599072f7b77dd6b28398956efa52e1a8981c566 Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Wed, 23 Sep 2026 22:24:25 -0400 Subject: [PATCH 1/7] For HyKKTSolver and SCCG, don't print convergence iteration count and final error by default. Log them with verbosity level misc. --- resolve/hykkt/HyKKTSolver.cpp | 10 ++++------ .../sccg/SchurComplementConjugateGradient.cpp | 7 ++++--- tests/unit/hykkt/HykktSCCGTests.hpp | 3 +++ tests/unit/hykkt/HykktSolverTests.hpp | 16 +++++++++++++++- 4 files changed, 26 insertions(+), 10 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index d83a2fc22..d1316b407 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -9,6 +9,8 @@ #include #include +#include + namespace ReSolve { using namespace constants; @@ -239,10 +241,6 @@ namespace ReSolve { // Assume all matrix blocks and RHS blocks are already set - std::cout << "H size: " << H_->getNumRows() << " " << H_->getNumColumns() << " " << H_->getNnz() << " \n"; - std::cout << "J size: " << J_->getNumRows() << " " << J_->getNumColumns() << " " << J_->getNnz() << " \n"; - std::cout << "D_s nnz = " << D_s_->getNnz() << "\n"; - if (!allocated_) { r_x_perm_ = new vector::Vector(n_x_); @@ -650,11 +648,11 @@ namespace ReSolve if (norm_r_x_sq > 0) { norm_res /= sqrt(norm_r_x_sq); - printf("||Ax-b||/||b|| = %32.32g\n\n", norm_res); + ReSolve::io::Logger::misc() << "||Ax-b||/||b|| = " << std::setprecision(16) << norm_res << std::endl; } else { - printf("||Ax-b|| = %32.32g\n\n", norm_res); + ReSolve::io::Logger::misc() << "||Ax-b|| = " << std::setprecision(16) << norm_res << std::endl; } allocated_ = true; diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index 791f8cd9f..64648a3b8 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -1,6 +1,7 @@ #include "SchurComplementConjugateGradient.hpp" #include +#include #include @@ -145,7 +146,7 @@ namespace ReSolve gamma_i1_ = vector_handler_->dot(r_, r_, memspace_); if (sqrt(gamma_i1_) < tol_) { - printf("Convergence occured at iteration %d\n", i); + ReSolve::io::Logger::misc() << "Convergence occured at iteration " << std::to_string(i) << std::endl; break; } matrix_handler_->matvec(J_tr_, r_, y_, &ONE, &ZERO, memspace_); @@ -157,10 +158,10 @@ namespace ReSolve alpha_ = gamma_i_ / (delta_ - beta_ * gamma_i_ / alpha_); } - printf("Conjugate gradient error is %32.32g \n", sqrt(gamma_i1_)); + ReSolve::io::Logger::misc() << "Conjugate gradient error is " << std::setprecision(16) << sqrt(gamma_i1_) << std::endl; if (i == itmax_) { - printf("No CG convergence in %d iterations\n", itmax_); + ReSolve::io::Logger::misc() << "No CG convergence in " << std::to_string(itmax_) << " iterations" << std::endl; return 1; } diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index f3958b310..d09d787c7 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -97,7 +98,9 @@ namespace ReSolve sccg.addMatrixInfo(J, J_tr); sccg.addVectorInfo(x_0, b); sccg.setup(); + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values int converged_n = sccg.solve(); // 0 if converged, 1 if not + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); TestStatus status; std::string testname(__func__); diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index b71232a91..d953de8f7 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -132,8 +133,11 @@ namespace ReSolve hykktSolver.setLHSPointers(x, s, y, y_d); hykktSolver.setGamma(gamma); hykktSolver.addHandlers(&matrixHandler_, &vectorHandler_); - + + + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type error = hykktSolver.solve(); + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); TestStatus status; std::string testname(__func__); @@ -185,7 +189,9 @@ namespace ReSolve // Change gamma to verify the cached SpGEMM coefficient is refreshed. hykktSolver.setGamma(gamma * 1.1); + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type second_error = hykktSolver.solve(); + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); status *= validateResult(second_error, tol); // Check that a zero RHS doesn't result in NaNs. @@ -193,7 +199,9 @@ namespace ReSolve r_s->setToZero(memspace_); r_y->setToZero(memspace_); r_yd->setToZero(memspace_); + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type zero_rhs_error = hykktSolver.solve(); + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); status *= validateResult(zero_rhs_error, tol); // Check that the solver raises an error when trying to change J_d @@ -213,10 +221,16 @@ namespace ReSolve no_jd_solver.setGamma(gamma); no_jd_solver.addHandlers(&matrixHandler_, &vectorHandler_); + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type no_jd_error = no_jd_solver.solve(); + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); + status *= validateResult(no_jd_error, tol); + + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type no_jd_reuse_error = no_jd_solver.solve(); + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); status *= validateResult(no_jd_reuse_error, tol); delete D_s_reuse; From 10f56f7992093a5247193f098cd47732bdb69e4a Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Thu, 24 Sep 2026 02:36:37 +0000 Subject: [PATCH 2/7] Apply pre-commmit fixes --- resolve/hykkt/HyKKTSolver.cpp | 4 ++-- tests/unit/hykkt/HykktSCCGTests.hpp | 2 +- tests/unit/hykkt/HykktSolverTests.hpp | 4 +--- 3 files changed, 4 insertions(+), 6 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index d1316b407..790a6177c 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -6,11 +6,11 @@ #include "HyKKTSolver.hpp" +#include + #include #include -#include - namespace ReSolve { using namespace constants; diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index d09d787c7..c515466f0 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -99,7 +99,7 @@ namespace ReSolve sccg.addVectorInfo(x_0, b); sccg.setup(); ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values - int converged_n = sccg.solve(); // 0 if converged, 1 if not + int converged_n = sccg.solve(); // 0 if converged, 1 if not ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); TestStatus status; diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index d953de8f7..5bcde7354 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -133,8 +133,7 @@ namespace ReSolve hykktSolver.setLHSPointers(x, s, y, y_d); hykktSolver.setGamma(gamma); hykktSolver.addHandlers(&matrixHandler_, &vectorHandler_); - - + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type error = hykktSolver.solve(); ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); @@ -227,7 +226,6 @@ namespace ReSolve status *= validateResult(no_jd_error, tol); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type no_jd_reuse_error = no_jd_solver.solve(); ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); From a8ca37b334cb4dc94d03f134b4f789ad0e8bdb70 Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Thu, 24 Sep 2026 11:12:30 -0400 Subject: [PATCH 3/7] Revert "For HyKKTSolver and SCCG, don't print convergence iteration count and final error by default. Log them with verbosity level misc." This reverts commit f599072f7b77dd6b28398956efa52e1a8981c566. --- resolve/hykkt/HyKKTSolver.cpp | 8 ++++++-- .../hykkt/sccg/SchurComplementConjugateGradient.cpp | 7 +++---- tests/unit/hykkt/HykktSCCGTests.hpp | 5 +---- tests/unit/hykkt/HykktSolverTests.hpp | 12 +----------- 4 files changed, 11 insertions(+), 21 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index 790a6177c..e3ec89bfa 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -241,6 +241,10 @@ namespace ReSolve { // Assume all matrix blocks and RHS blocks are already set + std::cout << "H size: " << H_->getNumRows() << " " << H_->getNumColumns() << " " << H_->getNnz() << " \n"; + std::cout << "J size: " << J_->getNumRows() << " " << J_->getNumColumns() << " " << J_->getNnz() << " \n"; + std::cout << "D_s nnz = " << D_s_->getNnz() << "\n"; + if (!allocated_) { r_x_perm_ = new vector::Vector(n_x_); @@ -648,11 +652,11 @@ namespace ReSolve if (norm_r_x_sq > 0) { norm_res /= sqrt(norm_r_x_sq); - ReSolve::io::Logger::misc() << "||Ax-b||/||b|| = " << std::setprecision(16) << norm_res << std::endl; + printf("||Ax-b||/||b|| = %32.32g\n\n", norm_res); } else { - ReSolve::io::Logger::misc() << "||Ax-b|| = " << std::setprecision(16) << norm_res << std::endl; + printf("||Ax-b|| = %32.32g\n\n", norm_res); } allocated_ = true; diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index 64648a3b8..791f8cd9f 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -1,7 +1,6 @@ #include "SchurComplementConjugateGradient.hpp" #include -#include #include @@ -146,7 +145,7 @@ namespace ReSolve gamma_i1_ = vector_handler_->dot(r_, r_, memspace_); if (sqrt(gamma_i1_) < tol_) { - ReSolve::io::Logger::misc() << "Convergence occured at iteration " << std::to_string(i) << std::endl; + printf("Convergence occured at iteration %d\n", i); break; } matrix_handler_->matvec(J_tr_, r_, y_, &ONE, &ZERO, memspace_); @@ -158,10 +157,10 @@ namespace ReSolve alpha_ = gamma_i_ / (delta_ - beta_ * gamma_i_ / alpha_); } - ReSolve::io::Logger::misc() << "Conjugate gradient error is " << std::setprecision(16) << sqrt(gamma_i1_) << std::endl; + printf("Conjugate gradient error is %32.32g \n", sqrt(gamma_i1_)); if (i == itmax_) { - ReSolve::io::Logger::misc() << "No CG convergence in " << std::to_string(itmax_) << " iterations" << std::endl; + printf("No CG convergence in %d iterations\n", itmax_); return 1; } diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index c515466f0..f3958b310 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -12,7 +12,6 @@ #include #include #include -#include #include #include @@ -98,9 +97,7 @@ namespace ReSolve sccg.addMatrixInfo(J, J_tr); sccg.addVectorInfo(x_0, b); sccg.setup(); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values - int converged_n = sccg.solve(); // 0 if converged, 1 if not - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); + int converged_n = sccg.solve(); // 0 if converged, 1 if not TestStatus status; std::string testname(__func__); diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index 5bcde7354..436bd2494 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -14,7 +14,6 @@ #include #include #include -#include #include #include @@ -134,9 +133,7 @@ namespace ReSolve hykktSolver.setGamma(gamma); hykktSolver.addHandlers(&matrixHandler_, &vectorHandler_); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type error = hykktSolver.solve(); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); TestStatus status; std::string testname(__func__); @@ -188,9 +185,7 @@ namespace ReSolve // Change gamma to verify the cached SpGEMM coefficient is refreshed. hykktSolver.setGamma(gamma * 1.1); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type second_error = hykktSolver.solve(); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); status *= validateResult(second_error, tol); // Check that a zero RHS doesn't result in NaNs. @@ -198,9 +193,7 @@ namespace ReSolve r_s->setToZero(memspace_); r_y->setToZero(memspace_); r_yd->setToZero(memspace_); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type zero_rhs_error = hykktSolver.solve(); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); status *= validateResult(zero_rhs_error, tol); // Check that the solver raises an error when trying to change J_d @@ -220,15 +213,12 @@ namespace ReSolve no_jd_solver.setGamma(gamma); no_jd_solver.addHandlers(&matrixHandler_, &vectorHandler_); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type no_jd_error = no_jd_solver.solve(); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); - status *= validateResult(no_jd_error, tol); + ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type no_jd_reuse_error = no_jd_solver.solve(); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::WARNINGS); status *= validateResult(no_jd_reuse_error, tol); delete D_s_reuse; From 08ad6eff30ae2898aa69da794ca2e80ea8cefe7e Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Thu, 24 Sep 2026 11:50:13 -0400 Subject: [PATCH 4/7] Removed all prints in solvers. Return relevant values and print inside the tests --- resolve/hykkt/HyKKTSolver.cpp | 34 ++----------------- .../sccg/SchurComplementConjugateGradient.cpp | 13 ++++--- tests/unit/hykkt/HykktSCCGTests.hpp | 14 ++++---- tests/unit/hykkt/HykktSolverTests.hpp | 18 ++++++++-- 4 files changed, 35 insertions(+), 44 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index e3ec89bfa..f2b918e7b 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -241,10 +241,6 @@ namespace ReSolve { // Assume all matrix blocks and RHS blocks are already set - std::cout << "H size: " << H_->getNumRows() << " " << H_->getNumColumns() << " " << H_->getNnz() << " \n"; - std::cout << "J size: " << J_->getNumRows() << " " << J_->getNumColumns() << " " << J_->getNnz() << " \n"; - std::cout << "D_s nnz = " << D_s_->getNnz() << "\n"; - if (!allocated_) { r_x_perm_ = new vector::Vector(n_x_); @@ -609,32 +605,18 @@ namespace ReSolve } /** - * @brief calculates the error of Ax - b + * @brief calculates the error of Ax - b (||Ax-b||, without normalization) * * @pre solution properly recovered using recoverSolution() * - * @param[out] norm_res - Error of Ax - b + * @return - Error of Ax - b */ real_type hykkt::HyKKTSolver::checkError() { - // Start of block, calculate error of Ax-b // Calculate error in r_x - real_type norm_r_x_sq = 0; - real_type norm_rs_sq = 0; - real_type norm_r_y_sq = 0; - real_type norm_r_yd_sq = 0; real_type norm_resx_sq = 0; real_type norm_resy_sq = 0; - // This will aggregate the squared norms of the residual and rhs - // Note that by construction the residuals of r_s and r_yd are 0 - norm_r_x_sq = vectorHandler_->dot(r_x_, r_x_, memspace_); - norm_rs_sq = vectorHandler_->dot(r_s_, r_s_, memspace_); - norm_r_y_sq = vectorHandler_->dot(r_y_copy_, r_y_copy_, memspace_); - norm_r_yd_sq = vectorHandler_->dot(r_yd_, r_yd_, memspace_); - - norm_r_x_sq += norm_rs_sq + norm_r_y_sq + norm_r_yd_sq; - matrixHandler_->matvec(H_, x_, r_x_, &MINUS_ONE, &ONE, memspace_); if (J_d_flag_) { @@ -648,19 +630,9 @@ namespace ReSolve norm_resy_sq = vectorHandler_->dot(r_y_copy_, r_y_copy_, memspace_); norm_resx_sq += norm_resy_sq; - real_type norm_res = sqrt(norm_resx_sq); - if (norm_r_x_sq > 0) - { - norm_res /= sqrt(norm_r_x_sq); - printf("||Ax-b||/||b|| = %32.32g\n\n", norm_res); - } - else - { - printf("||Ax-b|| = %32.32g\n\n", norm_res); - } allocated_ = true; - return norm_res; + return sqrt(norm_resx_sq); } } // namespace ReSolve diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index 791f8cd9f..0431ada15 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -101,6 +101,12 @@ namespace ReSolve w_->allocate(memspace_); } + /** + * @brief Solves the system + * @return The number of iterations required to converge, or -1 if solver failed to converge + * @pre x_0_ contains the initial guess + * @post x_0_ contains the solution + */ int SchurComplementConjugateGradient::solve() { using namespace constants; @@ -145,7 +151,6 @@ namespace ReSolve gamma_i1_ = vector_handler_->dot(r_, r_, memspace_); if (sqrt(gamma_i1_) < tol_) { - printf("Convergence occured at iteration %d\n", i); break; } matrix_handler_->matvec(J_tr_, r_, y_, &ONE, &ZERO, memspace_); @@ -157,14 +162,12 @@ namespace ReSolve alpha_ = gamma_i_ / (delta_ - beta_ * gamma_i_ / alpha_); } - printf("Conjugate gradient error is %32.32g \n", sqrt(gamma_i1_)); if (i == itmax_) { - printf("No CG convergence in %d iterations\n", itmax_); - return 1; + return -1; } - return 0; + return i + 1; // Add 1 to return the total number of iterations elapsed instead of the iteration index } } // namespace hykkt diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index f3958b310..a135b63f6 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -97,20 +97,22 @@ namespace ReSolve sccg.addMatrixInfo(J, J_tr); sccg.addVectorInfo(x_0, b); sccg.setup(); - int converged_n = sccg.solve(); // 0 if converged, 1 if not + int iters = sccg.solve(); // -1 if solver failed to converge + printf("Convergence occured after %d iteration(s) for non-trivial test system.\n", iters); TestStatus status; std::string testname(__func__); testname += " n=" + std::to_string(n) + ", m=" + std::to_string(m) + ", nnz =" + std::to_string(nnz); - status *= validateResult(x_0, converged_n); + status *= validateResult(x_0, iters); // A zero initial residual is already converged and must not enter // conjugate-gradient divisions with zero numerator and denominator. x_0->setToZero(memspace_); b->setToZero(memspace_); sccg.addVectorInfo(x_0, b); - int zero_residual_converged_n = sccg.solve(); - status *= (zero_residual_converged_n == 0); + int zero_residual_iters = sccg.solve(); + printf("Convergence occured after %d iteration(s) for zero-residual (trivial) test system.\n", zero_residual_iters); + status *= (zero_residual_iters == 0); status *= (vector_handler_.dot(x_0, x_0, memspace_) <= sccg_tol); delete H; @@ -144,9 +146,9 @@ namespace ReSolve * @brief Validate the SCCG result. * @param[in] x_0 Pointer to the output x_0 vector. */ - bool validateResult(vector::Vector* x_0, int converged_n) + bool validateResult(vector::Vector* x_0, int iters) { - if (converged_n != 0) + if (iters == -1) { return false; } diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index 436bd2494..81bde5c2a 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -116,6 +116,20 @@ namespace ReSolve r_yd->syncData(memory::DEVICE); } + // Get RHS norm + real_type norm_r_x_sq = 0; + real_type norm_r_s_sq = 0; + real_type norm_r_y_sq = 0; + real_type norm_r_yd_sq = 0; + + // This will aggregate the squared norms of the residual and rhs + // Note that by construction the residuals of r_s and r_yd are 0 + norm_r_x_sq = vectorHandler_.dot(r_x, r_x, memspace_); + norm_r_s_sq = vectorHandler_.dot(r_s, r_s, memspace_); + norm_r_y_sq = vectorHandler_.dot(r_y, r_y, memspace_); + norm_r_yd_sq = vectorHandler_.dot(r_yd, r_yd, memspace_); + real_type norm_rhs = sqrt(norm_r_x_sq + norm_r_s_sq + norm_r_y_sq + norm_r_yd_sq); + // LHS vector blocks vector::Vector* x = new vector::Vector(n_x); vector::Vector* s = new vector::Vector(m_d); @@ -133,7 +147,7 @@ namespace ReSolve hykktSolver.setGamma(gamma); hykktSolver.addHandlers(&matrixHandler_, &vectorHandler_); - real_type error = hykktSolver.solve(); + real_type error = hykktSolver.solve() / norm_rhs; TestStatus status; std::string testname(__func__); @@ -185,7 +199,7 @@ namespace ReSolve // Change gamma to verify the cached SpGEMM coefficient is refreshed. hykktSolver.setGamma(gamma * 1.1); - real_type second_error = hykktSolver.solve(); + real_type second_error = hykktSolver.solve() / norm_rhs; status *= validateResult(second_error, tol); // Check that a zero RHS doesn't result in NaNs. From 187b844e38f9465b595e0a65d0b33c1d94fd3d6c Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Thu, 24 Sep 2026 15:53:03 +0000 Subject: [PATCH 5/7] Apply pre-commmit fixes --- tests/unit/hykkt/HykktSolverTests.hpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index 81bde5c2a..d514a3e7d 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -118,16 +118,16 @@ namespace ReSolve // Get RHS norm real_type norm_r_x_sq = 0; - real_type norm_r_s_sq = 0; + real_type norm_r_s_sq = 0; real_type norm_r_y_sq = 0; real_type norm_r_yd_sq = 0; - + // This will aggregate the squared norms of the residual and rhs // Note that by construction the residuals of r_s and r_yd are 0 - norm_r_x_sq = vectorHandler_.dot(r_x, r_x, memspace_); - norm_r_s_sq = vectorHandler_.dot(r_s, r_s, memspace_); - norm_r_y_sq = vectorHandler_.dot(r_y, r_y, memspace_); - norm_r_yd_sq = vectorHandler_.dot(r_yd, r_yd, memspace_); + norm_r_x_sq = vectorHandler_.dot(r_x, r_x, memspace_); + norm_r_s_sq = vectorHandler_.dot(r_s, r_s, memspace_); + norm_r_y_sq = vectorHandler_.dot(r_y, r_y, memspace_); + norm_r_yd_sq = vectorHandler_.dot(r_yd, r_yd, memspace_); real_type norm_rhs = sqrt(norm_r_x_sq + norm_r_s_sq + norm_r_y_sq + norm_r_yd_sq); // LHS vector blocks @@ -230,7 +230,6 @@ namespace ReSolve real_type no_jd_error = no_jd_solver.solve(); status *= validateResult(no_jd_error, tol); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values real_type no_jd_reuse_error = no_jd_solver.solve(); status *= validateResult(no_jd_reuse_error, tol); From 9e626cba79a8d774b45c3943ae34a5020f8e7797 Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Fri, 25 Sep 2026 15:36:20 -0400 Subject: [PATCH 6/7] Add getters for error and itcount information. Print them only if tests fail. --- resolve/hykkt/HyKKTSolver.cpp | 8 ++-- resolve/hykkt/HyKKTSolver.hpp | 4 +- .../sccg/SchurComplementConjugateGradient.cpp | 37 ++++++++++++++++--- .../sccg/SchurComplementConjugateGradient.hpp | 5 +++ tests/unit/hykkt/HykktSCCGTests.hpp | 30 +++++++++------ tests/unit/hykkt/HykktSolverTests.hpp | 24 ++++++++---- 6 files changed, 79 insertions(+), 29 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index f2b918e7b..b742e7a4e 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -171,11 +171,11 @@ namespace ReSolve * * @pre matrix files have been loaded into the solver * - * @param[out] - Error of Ax - b + * @param[out] - Status code. 0 for successful solve * * @post solution to given KKT system is computed using Hykkt */ - real_type hykkt::HyKKTSolver::solve() + int hykkt::HyKKTSolver::solve() { setupParameters(); @@ -225,7 +225,7 @@ namespace ReSolve computeConjugateGradient(); recoverSolution(); - return checkError(); + return 0; } /** @@ -607,7 +607,7 @@ namespace ReSolve /** * @brief calculates the error of Ax - b (||Ax-b||, without normalization) * - * @pre solution properly recovered using recoverSolution() + * @pre solve() is called * * @return - Error of Ax - b */ diff --git a/resolve/hykkt/HyKKTSolver.hpp b/resolve/hykkt/HyKKTSolver.hpp index a1c93a07e..06f5a73d4 100644 --- a/resolve/hykkt/HyKKTSolver.hpp +++ b/resolve/hykkt/HyKKTSolver.hpp @@ -49,7 +49,8 @@ namespace ReSolve void addHandlers(MatrixHandler* matrixHandler, VectorHandler* vectorHandler); - real_type solve(); + int solve(); + real_type checkError(); private: // Intermediate steps of solving the system @@ -68,7 +69,6 @@ namespace ReSolve void setupConjugateGradient(); void computeConjugateGradient(); void recoverSolution(); - real_type checkError(); static constexpr int ruiz_its_ = 2; static constexpr double cholesky_tol_ = 1e-12; diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index 0431ada15..730df804a 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -103,9 +103,9 @@ namespace ReSolve /** * @brief Solves the system - * @return The number of iterations required to converge, or -1 if solver failed to converge + * @return 0 if the system converged. 1 if it did not converge * @pre x_0_ contains the initial guess - * @post x_0_ contains the solution + * @post x_0_ contains the solution, itcount_ contains the number of iterations elapsed until convergence */ int SchurComplementConjugateGradient::solve() { @@ -130,6 +130,8 @@ namespace ReSolve if (sqrt(gamma_i_) < tol_) { gamma_i1_ = gamma_i_; + error_ = sqrt(gamma_i_); + itcount_ = 0; return 0; } @@ -164,11 +166,36 @@ namespace ReSolve if (i == itmax_) { - return -1; + return 1; } - return i + 1; // Add 1 to return the total number of iterations elapsed instead of the iteration index + error_ = sqrt(gamma_i1_); + itcount_ = i + 1; // Add 1 to return the total number of iterations elapsed instead of the iteration index + return 0; + } + + /** + * @brief get the number of iterations elapsed until solver convergence + * + * @pre solve() is called + * + * @return - Iteration count + */ + int SchurComplementConjugateGradient::checkItcount() + { + return itcount_; + } + + /** + * @brief calculates the error of the solution + * + * @pre solve() is called + * + * @return - Error of Ax - b + */ + real_type SchurComplementConjugateGradient::checkError() + { + return error_; } - } // namespace hykkt } // namespace ReSolve diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp index bf0700983..391e3777b 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp @@ -56,6 +56,8 @@ namespace ReSolve */ void setup(); int solve(); + int checkItcount(); + real_type checkError(); private: index_type n_; // Dimension of outer system @@ -63,6 +65,9 @@ namespace ReSolve int itmax_ = 100; // Maximum iterations for conjugate gradient double tol_ = 1e-12; // Solver tolerance for Schur + int itcount_; // Number of iterations elapsed until convergence + real_type error_; // Error of the solution + CholeskySolver* choleskySolver_{nullptr}; // Cholesky factorization on 1,1 block MatrixHandler* matrix_handler_{nullptr}; ///< Backend-specific matrix handler. diff --git a/tests/unit/hykkt/HykktSCCGTests.hpp b/tests/unit/hykkt/HykktSCCGTests.hpp index a135b63f6..d182e1b8e 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -97,22 +97,34 @@ namespace ReSolve sccg.addMatrixInfo(J, J_tr); sccg.addVectorInfo(x_0, b); sccg.setup(); - int iters = sccg.solve(); // -1 if solver failed to converge - printf("Convergence occured after %d iteration(s) for non-trivial test system.\n", iters); + bool converged = (sccg.solve() == 0); + if (!converged) + { + real_type error = sccg.checkError(); + printf("Convergence failed to converge for non-trivial test system, with error %32.32g.\n", error); + } TestStatus status; std::string testname(__func__); testname += " n=" + std::to_string(n) + ", m=" + std::to_string(m) + ", nnz =" + std::to_string(nnz); - status *= validateResult(x_0, iters); + status *= converged; + if (converged) + { + status *= validateResult(x_0); + } // A zero initial residual is already converged and must not enter // conjugate-gradient divisions with zero numerator and denominator. x_0->setToZero(memspace_); b->setToZero(memspace_); sccg.addVectorInfo(x_0, b); - int zero_residual_iters = sccg.solve(); - printf("Convergence occured after %d iteration(s) for zero-residual (trivial) test system.\n", zero_residual_iters); - status *= (zero_residual_iters == 0); + bool zero_residual_converged = (sccg.solve() == 0); + if (!zero_residual_converged) + { + real_type error = sccg.checkError(); + printf("Convergence failed to converge for zero-initial-residual (trivial) test system, with error %32.32g.\n", error); + } + status *= zero_residual_converged; status *= (vector_handler_.dot(x_0, x_0, memspace_) <= sccg_tol); delete H; @@ -146,12 +158,8 @@ namespace ReSolve * @brief Validate the SCCG result. * @param[in] x_0 Pointer to the output x_0 vector. */ - bool validateResult(vector::Vector* x_0, int iters) + bool validateResult(vector::Vector* x_0) { - if (iters == -1) - { - return false; - } if (memspace_ == memory::DEVICE) { x_0->syncData(memory::HOST); diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index d514a3e7d..63e59f58d 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -147,7 +147,8 @@ namespace ReSolve hykktSolver.setGamma(gamma); hykktSolver.addHandlers(&matrixHandler_, &vectorHandler_); - real_type error = hykktSolver.solve() / norm_rhs; + hykktSolver.solve(); + real_type error = hykktSolver.checkError() / norm_rhs; TestStatus status; std::string testname(__func__); @@ -199,7 +200,8 @@ namespace ReSolve // Change gamma to verify the cached SpGEMM coefficient is refreshed. hykktSolver.setGamma(gamma * 1.1); - real_type second_error = hykktSolver.solve() / norm_rhs; + hykktSolver.solve(); + real_type second_error = hykktSolver.checkError() / norm_rhs; status *= validateResult(second_error, tol); // Check that a zero RHS doesn't result in NaNs. @@ -207,7 +209,8 @@ namespace ReSolve r_s->setToZero(memspace_); r_y->setToZero(memspace_); r_yd->setToZero(memspace_); - real_type zero_rhs_error = hykktSolver.solve(); + hykktSolver.solve(); + real_type zero_rhs_error = hykktSolver.checkError(); status *= validateResult(zero_rhs_error, tol); // Check that the solver raises an error when trying to change J_d @@ -227,11 +230,12 @@ namespace ReSolve no_jd_solver.setGamma(gamma); no_jd_solver.addHandlers(&matrixHandler_, &vectorHandler_); - real_type no_jd_error = no_jd_solver.solve(); + no_jd_solver.solve(); + real_type no_jd_error = no_jd_solver.checkError(); status *= validateResult(no_jd_error, tol); - ReSolve::io::Logger::setVerbosity(ReSolve::io::Logger::EVERYTHING); // Enable printing for solver convergence & error values - real_type no_jd_reuse_error = no_jd_solver.solve(); + no_jd_solver.solve(); + real_type no_jd_reuse_error = no_jd_solver.checkError(); status *= validateResult(no_jd_reuse_error, tol); delete D_s_reuse; @@ -265,7 +269,13 @@ namespace ReSolve */ bool validateResult(real_type error, real_type tol) { - return error < tol; + if (error >= tol) + { + printf("Solve failed! Error is %32.32g.\n", error); + return false; + } + + return true; } }; // class HykktSolverTests } // namespace tests From b3c8d6d7a6fb646a7da2768d10fc88e92c171e2e Mon Sep 17 00:00:00 2001 From: andrewxu319 Date: Fri, 25 Sep 2026 19:38:07 +0000 Subject: [PATCH 7/7] Apply pre-commmit fixes --- resolve/hykkt/HyKKTSolver.hpp | 32 +++++++++---------- .../sccg/SchurComplementConjugateGradient.cpp | 10 +++--- .../sccg/SchurComplementConjugateGradient.hpp | 10 +++--- tests/unit/hykkt/HykktSolverTests.hpp | 2 +- 4 files changed, 27 insertions(+), 27 deletions(-) diff --git a/resolve/hykkt/HyKKTSolver.hpp b/resolve/hykkt/HyKKTSolver.hpp index 06f5a73d4..142642be6 100644 --- a/resolve/hykkt/HyKKTSolver.hpp +++ b/resolve/hykkt/HyKKTSolver.hpp @@ -49,26 +49,26 @@ namespace ReSolve void addHandlers(MatrixHandler* matrixHandler, VectorHandler* vectorHandler); - int solve(); + int solve(); real_type checkError(); private: // Intermediate steps of solving the system - void setupParameters(); - void setupSpGEMMHtilde(); - void computeSpGEMMHtilde(); - void setupSolutionCheck(); - void setupRuizScaling(); - void computeRuizScaling(); - void setupSpGEMMHgamma(); - void computeSpGEMMHgamma(); - void setupPermutation(); - void applyPermutation(); - void setupHgammaFactorization(); // Uses Cholesky - void computeHgammaFactorization(); - void setupConjugateGradient(); - void computeConjugateGradient(); - void recoverSolution(); + void setupParameters(); + void setupSpGEMMHtilde(); + void computeSpGEMMHtilde(); + void setupSolutionCheck(); + void setupRuizScaling(); + void computeRuizScaling(); + void setupSpGEMMHgamma(); + void computeSpGEMMHgamma(); + void setupPermutation(); + void applyPermutation(); + void setupHgammaFactorization(); // Uses Cholesky + void computeHgammaFactorization(); + void setupConjugateGradient(); + void computeConjugateGradient(); + void recoverSolution(); static constexpr int ruiz_its_ = 2; static constexpr double cholesky_tol_ = 1e-12; diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp index 730df804a..04fdfd223 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp @@ -130,8 +130,8 @@ namespace ReSolve if (sqrt(gamma_i_) < tol_) { gamma_i1_ = gamma_i_; - error_ = sqrt(gamma_i_); - itcount_ = 0; + error_ = sqrt(gamma_i_); + itcount_ = 0; return 0; } @@ -169,11 +169,11 @@ namespace ReSolve return 1; } - error_ = sqrt(gamma_i1_); + error_ = sqrt(gamma_i1_); itcount_ = i + 1; // Add 1 to return the total number of iterations elapsed instead of the iteration index return 0; } - + /** * @brief get the number of iterations elapsed until solver convergence * @@ -185,7 +185,7 @@ namespace ReSolve { return itcount_; } - + /** * @brief calculates the error of the solution * diff --git a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp index 391e3777b..89ebe7e56 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp @@ -54,9 +54,9 @@ namespace ReSolve * * @pre Must be called exactly once before the first call to solve(). */ - void setup(); - int solve(); - int checkItcount(); + void setup(); + int solve(); + int checkItcount(); real_type checkError(); private: @@ -65,8 +65,8 @@ namespace ReSolve int itmax_ = 100; // Maximum iterations for conjugate gradient double tol_ = 1e-12; // Solver tolerance for Schur - int itcount_; // Number of iterations elapsed until convergence - real_type error_; // Error of the solution + int itcount_; // Number of iterations elapsed until convergence + real_type error_; // Error of the solution CholeskySolver* choleskySolver_{nullptr}; // Cholesky factorization on 1,1 block diff --git a/tests/unit/hykkt/HykktSolverTests.hpp b/tests/unit/hykkt/HykktSolverTests.hpp index 63e59f58d..823f52453 100644 --- a/tests/unit/hykkt/HykktSolverTests.hpp +++ b/tests/unit/hykkt/HykktSolverTests.hpp @@ -234,7 +234,7 @@ namespace ReSolve real_type no_jd_error = no_jd_solver.checkError(); status *= validateResult(no_jd_error, tol); - no_jd_solver.solve(); + no_jd_solver.solve(); real_type no_jd_reuse_error = no_jd_solver.checkError(); status *= validateResult(no_jd_reuse_error, tol);