Andrew/sccg and matrix class

CHANGELOG.md

@@ -2,7 +2,11 @@
 
 ## HyKKT Release changes
 
-- Added classes and tests for permutation, Ruiz scaling, Cholesky factorization, and matrix multiplication and addition.
+- Added classes and tests for permutation, Ruiz scaling, Cholesky factorization, Schur complement conjugate gradient and matrix multiplication and addition.
+
+- Changed random number generation int tests to be C++ style and fixed-seed, to avoid random failures.
+
+- Added Schur Complement Conjugate Gradient class.
 
 ## Changes to Re::Solve since release 0.99.2
 
@@ -22,6 +26,9 @@
 
 - Improved coding guidelines for developers on floating point conventions.
 
+- Added `SchurComplementConjugateGradient` with CPU implementation.
+
+- Created `Matrix`, a parent class for all matrix classes.
 - Added Developer Guide documentation for vector and matrix classes and handlers.
 
 ## Changes to Re::Solve in release 0.99.2

resolve/Common.hpp

@@ -14,11 +14,12 @@ namespace ReSolve
 
   namespace constants
   {
-    constexpr real_type ZERO      = 0.0;
-    constexpr real_type ONE       = 1.0;
-    constexpr real_type TWO       = 2.0;
-    constexpr real_type HALF      = 0.5;
-    constexpr real_type MINUS_ONE = -1.0;
+    constexpr real_type  ZERO      = 0.0;
+    constexpr real_type  ONE       = 1.0;
+    constexpr real_type  TWO       = 2.0;
+    constexpr real_type  HALF      = 0.5;
+    constexpr real_type  MINUS_ONE = -1.0;
+    constexpr index_type SEED      = 12345;
 
     constexpr real_type MACHINE_EPSILON = std::numeric_limits<real_type>::epsilon();
   } // namespace constants

resolve/hykkt/CMakeLists.txt

@@ -2,3 +2,4 @@ add_subdirectory(permutation)
 add_subdirectory(ruiz)
 add_subdirectory(cholesky)
 add_subdirectory(spgemm)
+add_subdirectory(sccg)

resolve/hykkt/cholesky/CholeskySolverHip.cpp

@@ -60,6 +60,7 @@ namespace ReSolve
       }
       A_chol_ = convertToCholmod(A);
       A_      = A;
+      mem_.deviceSynchronize();
     }
 
     void CholeskySolverHip::symbolicAnalysis()
@@ -70,6 +71,7 @@ namespace ReSolve
       {
         out::error() << "Cholesky symbolic analysis failed with status: " << Common_.status << "\n";
       }
+      mem_.deviceSynchronize();
     }
 
     /**
@@ -156,6 +158,7 @@ namespace ReSolve
           out::error() << "Refactorization step failed with status: " << status << "\n";
         }
         L_->setUpdated(memory::DEVICE);
+        mem_.deviceSynchronize();
       }
     }
 

resolve/hykkt/sccg/CMakeLists.txt

@@ -0,0 +1,40 @@
+# Add source files for the SCCG module
+set(SCCG_SRC SchurComplementConjugateGradient.cpp)
+set(SCCG_HEADER_INSTALL SchurComplementConjugateGradient.hpp)
+
+# Build shared library ReSolve::hykkt
+add_library(resolve_hykkt_sccg SHARED ${SCCG_SRC})
+
+target_link_libraries(
+  resolve_hykkt_sccg
+  PRIVATE resolve_hykkt resolve_hykkt_chol resolve_logger resolve_vector
+          resolve_matrix ${suitesparse_cholmod}
+)
+
+# Link to CUDA ReSolve backend if CUDA is support enabled
+
+target_include_directories(resolve_hykkt_sccg PUBLIC ${SUITESPARSE_INCLUDE_DIR})
+
+if(RESOLVE_USE_CUDA)
+  target_sources(resolve_hykkt_sccg PRIVATE ${Matrix_CUDASDK_SRC})
+  target_link_libraries(resolve_hykkt_sccg PUBLIC resolve_backend_cuda)
+endif()
+
+# Link to HIP ReSolve backend if HIP is support enabled
+
+if(RESOLVE_USE_HIP)
+  target_sources(resolve_hykkt_sccg PRIVATE ${Matrix_ROCM_SRC})
+  target_link_libraries(resolve_hykkt_sccg PUBLIC resolve_backend_hip)
+endif()
+
+# Link to dummy device backend if GPU support is not enabled
+if(NOT RESOLVE_USE_GPU)
+  target_link_libraries(resolve_hykkt_sccg PUBLIC resolve_backend_cpu)
+endif()
+
+target_include_directories(
+  resolve_hykkt_sccg INTERFACE $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}>
+                               $<INSTALL_INTERFACE:include>
+)
+
+install(FILES ${SCCG_HEADER_INSTALL} DESTINATION include/resolve/hykkt)

resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp

@@ -0,0 +1,153 @@
+#include "SchurComplementConjugateGradient.hpp"
+
+#include <cmath>
+
+#include <resolve/Common.hpp>
+
+namespace ReSolve
+{
+  namespace hykkt
+  {
+    /** Constructor for SchurComplementConjugateGradient.
+     *  @param n[in] - Dimension of outer system.
+     *  @param m[in] - Dimension of inner system.
+     *  @param choleskySolver[in] - Factorization of Hgamma to use for direct solve.
+     *  @param memspace[in] - Memory space of incoming data and for computation.
+     *  @param matrix_handler[in] - Matrix handler for the selected backend.
+     *  @param vector_handler[in] - Vector handler for the selected backend.
+     */
+    SchurComplementConjugateGradient::SchurComplementConjugateGradient(
+        index_type          n,
+        index_type          m,
+        CholeskySolver*     choleskySolver,
+        memory::MemorySpace memspace,
+        MatrixHandler&      matrix_handler,
+        VectorHandler&      vector_handler)
+      : n_(n),
+        m_(m),
+        matrix_handler_(matrix_handler),
+        vector_handler_(vector_handler),
+        choleskySolver_(choleskySolver),
+        y_(m_),
+        z_(m_),
+        r_(n_),
+        p_(n_),
+        s_(n_),
+        w_(n_),
+        memspace_(memspace)
+    {
+      ;
+    }
+
+    /**
+     * @brief Loads or reloads matrix pointers to the solver
+     * @param[in] jc - Pointer to the JC matrix in CSR format.
+     * @param[in] jc_tr - Pointer to the transposed JC matrix in CSR format.
+     */
+    void SchurComplementConjugateGradient::addMatrixInfo(matrix::Csr* jc, matrix::Csr* jc_tr)
+    {
+      jc_    = jc;
+      jc_tr_ = jc_tr;
+    }
+
+    /**
+     * @brief Loads or reloads vector pointers to the solver
+     * @param[in] x0 - Pointer to the left-hand side vector.
+     * @param[in] b - Pointer to the right-hand side vector.
+     */
+    void SchurComplementConjugateGradient::addVectorInfo(vector::Vector* x0, vector::Vector* b)
+    {
+      x0_ = x0;
+      b_  = b;
+    }
+
+    /**
+     * @brief Reloads pointer to the Cholesky solver
+     * @param[in] choleskySolver - Factorization of Hgamma to use for direct solve.
+     */
+    void SchurComplementConjugateGradient::updateCholeskySolver(CholeskySolver* choleskySolver)
+    {
+      choleskySolver_ = choleskySolver;
+    }
+
+    void SchurComplementConjugateGradient::setSolverTolerance(double tol)
+    {
+      tol_ = tol;
+    }
+
+    void SchurComplementConjugateGradient::setSolverItmax(int itmax)
+    {
+      itmax_ = itmax;
+    }
+
+    void SchurComplementConjugateGradient::setup()
+    {
+      y_.allocate(memspace_);
+      z_.allocate(memspace_);
+      r_.allocate(memspace_);
+      p_.allocate(memspace_);
+      s_.allocate(memspace_);
+      w_.allocate(memspace_);
+
+      y_.setToZero(memspace_);
+      z_.setToZero(memspace_);
+      r_.setToZero(memspace_);
+      p_.setToZero(memspace_);
+      s_.setToZero(memspace_);
+      w_.setToZero(memspace_);
+
+      beta_ = 0;
+    }
+
+    int SchurComplementConjugateGradient::solve()
+    {
+      using namespace constants;
+
+      matrix_handler_.matvec(jc_tr_, x0_, &y_, &ONE, &ZERO, memspace_);
+      choleskySolver_->solve(&z_, &y_);
+      matrix_handler_.matvec(jc_, &z_, &r_, &MINUS_ONE, &ONE, memspace_);
+      gam_i_ = vector_handler_.dot(&r_, &r_, memspace_);
+
+      matrix_handler_.matvec(jc_tr_, &r_, &y_, &ONE, &ZERO, memspace_);
+      choleskySolver_->solve(&z_, &y_);
+      matrix_handler_.matvec(jc_, &z_, &w_, &ONE, &ZERO, memspace_);
+      delta_ = vector_handler_.dot(&w_, &r_, memspace_);
+      alpha_ = gam_i_ / delta_;
+
+      int i;
+      for (i = 0; i < itmax_; i++)
+      {
+        vector_handler_.scal(beta_, &p_, memspace_);
+        vector_handler_.axpy(ONE, &r_, &p_, memspace_);
+        vector_handler_.scal(beta_, &s_, memspace_);
+        vector_handler_.axpy(ONE, &w_, &s_, memspace_);
+        vector_handler_.axpy(alpha_, &p_, x0_, memspace_);
+        minalpha_ = -alpha_;
+        vector_handler_.axpy(minalpha_, &s_, &r_, memspace_);
+        gam_i1_ = vector_handler_.dot(&r_, &r_, memspace_);
+        if (sqrt(gam_i1_) < tol_)
+        {
+          printf("Convergence occured at iteration %d\n", i);
+          break;
+        }
+        matrix_handler_.matvec(jc_tr_, &r_, &y_, &ONE, &ZERO, memspace_);
+        choleskySolver_->solve(&z_, &y_);
+        matrix_handler_.matvec(jc_, &z_, &w_, &ONE, &ZERO, memspace_);
+        delta_ = vector_handler_.dot(&w_, &r_, memspace_);
+        beta_  = gam_i1_ / gam_i_;
+        gam_i_ = gam_i1_;
+        alpha_ = gam_i_ / (delta_ - beta_ * gam_i_ / alpha_);
+      }
+
+      printf("Error is %32.32g \n", sqrt(gam_i1_));
+      if (i == itmax_)
+      {
+        printf("No CG convergence in %d iterations\n", itmax_);
+        return 1;
+      }
+
+      return 0;
+    }
+
+  } // namespace hykkt
+} // namespace ReSolve

resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp

@@ -0,0 +1,85 @@
+/**
+ * @file SchurComplementConjugateGradient.hpp
+ * @brief Schur complement conjugate gradient solver for HyKKT.
+ */
+
+#pragma once
+
+#include <resolve/Common.hpp>
+#include <resolve/MemoryUtils.hpp>
+#include <resolve/hykkt/cholesky/CholeskySolver.hpp>
+#include <resolve/matrix/Csr.hpp>
+#include <resolve/matrix/MatrixHandler.hpp>
+#include <resolve/vector/Vector.hpp>
+#include <resolve/vector/VectorHandler.hpp>
+
+namespace ReSolve
+{
+  using index_type = ReSolve::index_type;
+  using real_type  = ReSolve::real_type;
+
+  namespace hykkt
+  {
+    class SchurComplementConjugateGradient
+    {
+    public:
+      /**
+       * @brief Constructor for SchurComplementConjugateGradient.
+       *
+       * The solver uses caller-provided matrix and vector handlers so the same solver can be run with CPU, CUDA, or HIP backends.
+       *
+       * @param[in] n Dimension of outer system.
+       * @param[in] m Dimension of inner system.
+       * @param[in] choleskySolver Factorization of Hgamma to use for direct solves.
+       * @param[in] memspace Memory space of incoming data and for computation.
+       * @param[in] matrix_handler Matrix handler for the selected backend.
+       * @param[in] vector_handler Vector handler for the selected backend.
+       */
+      SchurComplementConjugateGradient(index_type n, index_type m, CholeskySolver* choleskySolver, memory::MemorySpace memspace, MatrixHandler& matrix_handler, VectorHandler& vector_handler);
+
+      void addMatrixInfo(matrix::Csr* jc, matrix::Csr* jc_tr);
+      void addVectorInfo(vector::Vector* x0, vector::Vector* b);
+      void updateCholeskySolver(CholeskySolver* choleskySolver);
+      void setSolverTolerance(double tol);
+      void setSolverItmax(int itmax);
+
+      void setup();
+      int  solve();
+
+    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
+
+      MatrixHandler& matrix_handler_; ///< Backend-specific matrix handler.
+      VectorHandler& vector_handler_; ///< Backend-specific vector handler.
+
+      CholeskySolver* choleskySolver_; // Cholesky factorization on 1,1 block
+
+      matrix::Csr* jc_;
+      matrix::Csr* jc_tr_;
+
+      vector::Vector* x0_; // LHS of entire system
+      vector::Vector* b_;  // RHS of entire system
+
+      // scalars used for conjugate gradient
+      double beta_;
+      double delta_;
+      double alpha_;
+      double minalpha_;
+      double gam_i_;
+      double gam_i1_;
+
+      // Vectors used for conjugate gradient
+      vector::Vector y_; // Internal RHS of system
+      vector::Vector z_; // Internal LHS of system
+      vector::Vector r_; // Residual
+      vector::Vector p_;
+      vector::Vector s_;
+      vector::Vector w_;
+
+      memory::MemorySpace memspace_;
+    }; // class SchurComplementConjugateGradient
+  } // namespace hykkt
+} // namespace ReSolve