Disable debug info printing for HykktSolver

resolve/hykkt/HyKKTSolver.cpp

@@ -6,6 +6,8 @@
 
 #include "HyKKTSolver.hpp"
 
+#include <iomanip>
+
 #include <resolve/matrix/io.hpp>
 #include <resolve/utilities/logger/Logger.hpp>
 
@@ -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

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;

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

resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp

@@ -54,15 +54,20 @@ 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
       index_type m_;             // Dimension of inner system
       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.

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);

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,15 +200,17 @@ 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.
         r_x->setToZero(memspace_);
         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