Make SCCG test path CPU, CUDA, and HIP capable

CHANGELOG.md

@@ -4,6 +4,10 @@
 
 - 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
 
 - Added cmake-format.

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/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

@@ -13,9 +13,7 @@ target_link_libraries(
 
 # Link to CUDA ReSolve backend if CUDA is support enabled
 
-target_include_directories(
-  resolve_hykkt_spgemm PUBLIC ${SUITESPARSE_INCLUDE_DIR}
-)
+target_include_directories(resolve_hykkt_sccg PUBLIC ${SUITESPARSE_INCLUDE_DIR})
 
 if(RESOLVE_USE_CUDA)
   target_sources(resolve_hykkt_sccg PRIVATE ${Matrix_CUDASDK_SRC})

resolve/hykkt/sccg/SchurComplementConjugateGradient.cpp

@@ -13,26 +13,28 @@ namespace ReSolve
      *  @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)
+        memory::MemorySpace memspace,
+        MatrixHandler&      matrix_handler,
+        VectorHandler&      vector_handler)
       : n_(n),
         m_(m),
+        matrix_handler_(matrix_handler),
+        vector_handler_(vector_handler),
         choleskySolver_(choleskySolver),
-        memspace_(memspace),
-        mem_(),
-        workspace_(),
-        matrixhandler_(&workspace_),
-        vectorhandler_(&workspace_),
         y_(m_),
         z_(m_),
         r_(n_),
         p_(n_),
         s_(n_),
-        w_(n_)
+        w_(n_),
+        memspace_(memspace)
     {
       ;
     }
@@ -101,37 +103,37 @@ namespace ReSolve
     {
       using namespace constants;
 
-      matrixhandler_.matvec(jc_tr_, x0_, &y_, &ONE, &ZERO, memspace_);
+      matrix_handler_.matvec(jc_tr_, x0_, &y_, &ONE, &ZERO, memspace_);
       choleskySolver_->solve(&z_, &y_);
-      matrixhandler_.matvec(jc_, &z_, &r_, &MINUS_ONE, &ONE, memspace_);
-      gam_i_ = vectorhandler_.dot(&r_, &r_, memspace_);
+      matrix_handler_.matvec(jc_, &z_, &r_, &MINUS_ONE, &ONE, memspace_);
+      gam_i_ = vector_handler_.dot(&r_, &r_, memspace_);
 
-      matrixhandler_.matvec(jc_tr_, &r_, &y_, &ONE, &ZERO, memspace_);
+      matrix_handler_.matvec(jc_tr_, &r_, &y_, &ONE, &ZERO, memspace_);
       choleskySolver_->solve(&z_, &y_);
-      matrixhandler_.matvec(jc_, &z_, &w_, &ONE, &ZERO, memspace_);
-      delta_ = vectorhandler_.dot(&w_, &r_, memspace_);
+      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++)
       {
-        vectorhandler_.scal(beta_, &p_, memspace_);
-        vectorhandler_.axpy(ONE, &r_, &p_, memspace_);
-        vectorhandler_.scal(beta_, &s_, memspace_);
-        vectorhandler_.axpy(ONE, &w_, &s_, memspace_);
-        vectorhandler_.axpy(alpha_, &p_, x0_, memspace_);
+        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_;
-        vectorhandler_.axpy(minalpha_, &s_, &r_, memspace_);
-        gam_i1_ = vectorhandler_.dot(&r_, &r_, memspace_);
+        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;
         }
-        matrixhandler_.matvec(jc_tr_, &r_, &y_, &ONE, &ZERO, memspace_);
+        matrix_handler_.matvec(jc_tr_, &r_, &y_, &ONE, &ZERO, memspace_);
         choleskySolver_->solve(&z_, &y_);
-        matrixhandler_.matvec(jc_, &z_, &w_, &ONE, &ZERO, memspace_);
-        delta_ = vectorhandler_.dot(&w_, &r_, memspace_);
+        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_);

resolve/hykkt/sccg/SchurComplementConjugateGradient.hpp

@@ -1,6 +1,6 @@
 /**
  * @file SchurComplementConjugateGradient.hpp
- * @brief Currently CPU-only.
+ * @brief Schur complement conjugate gradient solver for HyKKT.
  */
 
 #pragma once
@@ -12,7 +12,6 @@
 #include <resolve/matrix/MatrixHandler.hpp>
 #include <resolve/vector/Vector.hpp>
 #include <resolve/vector/VectorHandler.hpp>
-#include <resolve/workspace/LinAlgWorkspace.hpp>
 
 namespace ReSolve
 {
@@ -24,7 +23,19 @@ namespace ReSolve
     class SchurComplementConjugateGradient
     {
     public:
-      SchurComplementConjugateGradient(index_type n, index_type m, CholeskySolver* choleskySolver, memory::MemorySpace memspace);
+      /**
+       * @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);
@@ -41,15 +52,8 @@ namespace ReSolve
       int        itmax_ = 100;   // Maximum iterations for conjugate gradient
       double     tol_   = 1e-12; // Solver tolerance for Schur
 
-#ifdef RESOLVE_USE_CUDA
-      LinAlgWorkspaceCUDA workspace_;
-#elif defined(RESOLVE_USE_HIP)
-      LinAlgWorkspaceHIP workspace_;
-#else
-      LinAlgWorkspaceCpu workspace_;
-#endif
-      MatrixHandler matrixhandler_;
-      VectorHandler vectorhandler_;
+      MatrixHandler& matrix_handler_; ///< Backend-specific matrix handler.
+      VectorHandler& vector_handler_; ///< Backend-specific vector handler.
 
       CholeskySolver* choleskySolver_; // Cholesky factorization on 1,1 block
 
@@ -75,7 +79,6 @@ namespace ReSolve
       vector::Vector s_;
       vector::Vector w_;
 
-      MemoryHandler       mem_;
       memory::MemorySpace memspace_;
     }; // class SchurComplementConjugateGradient
   } // namespace hykkt

resolve/matrix/MatrixHandlerCuda.cpp

@@ -76,17 +76,23 @@ namespace ReSolve
     // result = alpha *A*x + beta * result
     cusparseStatus_t     status;
     cusparseDnVecDescr_t vecx = workspace_->getVecX();
-    cusparseCreateDnVec(&vecx, A->getNumRows(), vec_x->getData(memory::DEVICE), CUDA_R_64F);
+    cusparseCreateDnVec(&vecx, A->getNumColumns(), vec_x->getData(memory::DEVICE), CUDA_R_64F);
 
     cusparseDnVecDescr_t vecAx = workspace_->getVecY();
     cusparseCreateDnVec(&vecAx, A->getNumRows(), vec_result->getData(memory::DEVICE), CUDA_R_64F);
 
-    cusparseSpMatDescr_t matA = workspace_->getSpmvMatrixDescriptor();
-
-    void*            buffer_spmv     = workspace_->getSpmvBuffer();
     cusparseHandle_t handle_cusparse = workspace_->getCusparseHandle();
-    if (values_changed_)
+    // Rebuild cached SpMV setup if the matrix object or dimensions changed.
+    bool matrix_changed = (matrix_for_matvec_ != A) || (matvec_num_rows_ != A->getNumRows()) || (matvec_num_cols_ != A->getNumColumns()) || (matvec_nnz_ != A->getNnz());
+    if (matrix_changed || values_changed_)
+    {
+      workspace_->resetMatvecSetup();
+    }
+    cusparseSpMatDescr_t matA        = workspace_->getSpmvMatrixDescriptor();
+    void*                buffer_spmv = workspace_->getSpmvBuffer();
+    if (!workspace_->matvecSetup())
     {
+      // setup first, allocate, etc.
       status = cusparseCreateCsr(&matA,
                                  A->getNumRows(),
                                  A->getNumColumns(),
@@ -99,11 +105,6 @@ namespace ReSolve
                                  CUSPARSE_INDEX_BASE_ZERO,
                                  CUDA_R_64F);
       error_sum += status;
-      values_changed_ = false;
-    }
-    if (!workspace_->matvecSetup())
-    {
-      // setup first, allocate, etc.
       size_t bufferSize = 0;
 
       status = cusparseSpMV_bufferSize(handle_cusparse,
@@ -123,6 +124,13 @@ namespace ReSolve
       workspace_->setSpmvBuffer(buffer_spmv);
 
       workspace_->matvecSetupDone();
+
+      matrix_for_matvec_ = A;
+      matvec_num_rows_   = A->getNumRows();
+      matvec_num_cols_   = A->getNumColumns();
+      matvec_nnz_        = A->getNnz();
+
+      values_changed_ = false;
     }
 
     status = cusparseSpMV(handle_cusparse,

resolve/matrix/MatrixHandlerCuda.hpp

@@ -56,7 +56,12 @@ namespace ReSolve
 
   private:
     LinAlgWorkspaceCUDA* workspace_{nullptr};
-    bool                 values_changed_{true}; ///< needed for matvec
+    bool                 values_changed_{true}; ///< Flag to indicate if matrix values changed since the cached SpMV setup.
+
+    matrix::Sparse* matrix_for_matvec_{nullptr}; ///< Matrix used for cached SpMV setup.
+    index_type      matvec_num_rows_{0};         ///< Number of rows in cached SpMV matrix.
+    index_type      matvec_num_cols_{0};         ///< Number of columns in cached SpMV matrix.
+    index_type      matvec_nnz_{0};              ///< Number of nonzeros in cached SpMV matrix.
 
     MemoryHandler mem_; ///< Device memory manager object
   };

resolve/matrix/MatrixHandlerHip.cpp

@@ -78,6 +78,13 @@ namespace ReSolve
 
     rocsparse_handle handle_rocsparse = workspace_->getRocsparseHandle();
 
+    // Rebuild cached SpMV setup if the matrix object or dimensions changed.
+    bool matrix_changed = (matrix_for_matvec_ != A) || (matvec_num_rows_ != A->getNumRows()) || (matvec_num_cols_ != A->getNumColumns()) || (matvec_nnz_ != A->getNnz());
+    if (matrix_changed || values_changed_)
+    {
+      workspace_->resetMatvecSetup();
+    }
+
     rocsparse_mat_info  infoA  = workspace_->getSpmvMatrixInfo();
     rocsparse_mat_descr descrA = workspace_->getSpmvMatrixDescriptor();
 
@@ -106,6 +113,12 @@ namespace ReSolve
       workspace_->setSpmvMatrixDescriptor(descrA);
       workspace_->setSpmvMatrixInfo(infoA);
       workspace_->matvecSetupDone();
+
+      matrix_for_matvec_ = A;
+      matvec_num_rows_   = A->getNumRows();
+      matvec_num_cols_   = A->getNumColumns();
+      matvec_nnz_        = A->getNnz();
+      values_changed_    = false;
     }
 
     status = rocsparse_dcsrmv(handle_rocsparse,

resolve/matrix/MatrixHandlerHip.hpp

@@ -56,7 +56,12 @@ namespace ReSolve
 
   private:
     LinAlgWorkspaceHIP* workspace_{nullptr};
-    bool                values_changed_{true}; ///< needed for matvec
+    bool                values_changed_{true}; ///< Flag to indicate if matrix values changed since the cached SpMV setup.
+
+    matrix::Sparse* matrix_for_matvec_{nullptr}; ///< Matrix used for cached SpMV setup.
+    index_type      matvec_num_rows_{0};         ///< Number of rows in cached SpMV matrix.
+    index_type      matvec_num_cols_{0};         ///< Number of columns in cached SpMV matrix.
+    index_type      matvec_nnz_{0};              ///< Number of nonzeros in cached SpMV matrix.
 
     MemoryHandler mem_; ///< Device memory manager object
   };

resolve/workspace/LinAlgWorkspaceCUDA.cpp

@@ -204,6 +204,23 @@ namespace ReSolve
     matvec_setup_done_ = true;
   }
 
+  /**
+   * @brief Reset cached SpMV resources.
+   */
+  void LinAlgWorkspaceCUDA::resetMatvecSetup()
+  {
+    if (matvec_setup_done_)
+    {
+      cusparseDestroySpMat(mat_A_);
+      matvec_setup_done_ = false;
+    }
+    if (buffer_spmv_ != nullptr)
+    {
+      mem_.deleteOnDevice(buffer_spmv_);
+      buffer_spmv_ = nullptr;
+    }
+  }
+
   void LinAlgWorkspaceCUDA::initializeHandles()
   {
     cusparseCreate(&handle_cusparse_);

resolve/workspace/LinAlgWorkspaceCUDA.hpp

@@ -47,6 +47,13 @@ namespace ReSolve
 
     bool matvecSetup();
     void matvecSetupDone();
+    /**
+     * @brief Reset the cached CUDA SpMV setup.
+     *
+     * Destroys the cached sparse matrix descriptor and frees the SpMV buffer so
+     * the next matvec call can rebuild the SpMV setup if the matrix or its dimensions have changed.
+     */
+    void resetMatvecSetup();
 
   private:
     // handles