diff --git a/resolve/hykkt/HyKKTSolver.cpp b/resolve/hykkt/HyKKTSolver.cpp index d83a2fc22..b742e7a4e 100644 --- a/resolve/hykkt/HyKKTSolver.cpp +++ b/resolve/hykkt/HyKKTSolver.cpp @@ -6,6 +6,8 @@ #include "HyKKTSolver.hpp" +#include + #include #include @@ -169,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(); @@ -223,7 +225,7 @@ namespace ReSolve computeConjugateGradient(); recoverSolution(); - return checkError(); + return 0; } /** @@ -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_); @@ -607,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() + * @pre solve() is called * - * @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_) { @@ -646,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/HyKKTSolver.hpp b/resolve/hykkt/HyKKTSolver.hpp index a1c93a07e..142642be6 100644 --- a/resolve/hykkt/HyKKTSolver.hpp +++ b/resolve/hykkt/HyKKTSolver.hpp @@ -49,26 +49,26 @@ namespace ReSolve void addHandlers(MatrixHandler* matrixHandler, VectorHandler* vectorHandler); - real_type 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(); - real_type checkError(); + 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 791f8cd9f..04fdfd223 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 0 if the system converged. 1 if it did not converge + * @pre x_0_ contains the initial guess + * @post x_0_ contains the solution, itcount_ contains the number of iterations elapsed until convergence + */ int SchurComplementConjugateGradient::solve() { using namespace constants; @@ -124,6 +130,8 @@ namespace ReSolve if (sqrt(gamma_i_) < tol_) { gamma_i1_ = gamma_i_; + error_ = sqrt(gamma_i_); + itcount_ = 0; return 0; } @@ -145,7 +153,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,15 +164,38 @@ 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; } + 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..89ebe7e56 100644 --- a/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp +++ b/resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp @@ -54,8 +54,10 @@ namespace ReSolve * * @pre Must be called exactly once before the first call to solve(). */ - void setup(); - int solve(); + 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 f3958b310..d182e1b8e 100644 --- a/tests/unit/hykkt/HykktSCCGTests.hpp +++ b/tests/unit/hykkt/HykktSCCGTests.hpp @@ -97,20 +97,34 @@ 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 + 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, converged_n); + 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_converged_n = sccg.solve(); - status *= (zero_residual_converged_n == 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; @@ -144,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 converged_n) + bool validateResult(vector::Vector* x_0) { - if (converged_n != 0) - { - 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 b71232a91..823f52453 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,8 @@ namespace ReSolve hykktSolver.setGamma(gamma); hykktSolver.addHandlers(&matrixHandler_, &vectorHandler_); - real_type error = hykktSolver.solve(); + hykktSolver.solve(); + real_type error = hykktSolver.checkError() / norm_rhs; TestStatus status; std::string testname(__func__); @@ -185,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(); + 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. @@ -193,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 @@ -213,10 +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); - 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; @@ -250,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