Declutter functionality test

examples/CMakeLists.txt

@@ -133,19 +133,22 @@ set(CONSUMER_PATH ${CMAKE_INSTALL_PREFIX}/share/examples)
 install(PROGRAMS test.sh DESTINATION ${CONSUMER_PATH})
 
 # Select consumer app
-if(RESOLVE_USE_CUDA)
-  set(RESOLVE_CONSUMER_APP "testKLU_Rf_FGMRES.cpp")
-elseif(RESOLVE_USE_HIP)
-  set(RESOLVE_CONSUMER_APP "testKLU_RocSolver.cpp")
-elseif(RESOLVE_USE_KLU)
-  set(RESOLVE_CONSUMER_APP "testKLU.cpp")
-else()
-  set(RESOLVE_CONSUMER_APP "testVersion.cpp")
-endif()
+if(RESOLVE_USE_KLU)
+  if(RESOLVE_USE_CUDA)
+    set(RESOLVE_CONSUMER_APP "testRefactor.cpp")
+  elseif(RESOLVE_USE_HIP)
+    set(RESOLVE_CONSUMER_APP "testRefactor.cpp")
+  else()
+    set(RESOLVE_CONSUMER_APP "testKlu.cpp")
+  endif()
+else(RESOLVE_USE_KLU)
+  set(RESOLVE_CONSUMER_APP "testSysGmres.cpp")
+endif(RESOLVE_USE_KLU)
 
 # Install directory with example on how to consume ReSolve
 install(DIRECTORY resolve_consumer DESTINATION share/examples)
 install(FILES ${PROJECT_SOURCE_DIR}/tests/functionality/${RESOLVE_CONSUMER_APP} DESTINATION share/examples/resolve_consumer RENAME consumer.cpp)
+install(FILES ${PROJECT_SOURCE_DIR}/tests/functionality/TestHelper.hpp DESTINATION share/examples/resolve_consumer)
 
 # Shell script argumets: 
 #    1. Path to where resolve is installed.

examples/r_SysSolverHipRefine.cpp

@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
   vector_type* vec_x   = nullptr;
 
   ReSolve::SystemSolver* solver = new ReSolve::SystemSolver(workspace_HIP);
-  solver->setRefinementMethod("fgmres", "CGS2");
+  solver->setRefinementMethod("fgmres", "cgs2");
 
   for (int i = 0; i < numSystems; ++i)
   {

examples/resolve_consumer/CMakeLists.txt

@@ -18,4 +18,4 @@ target_link_libraries(consume.exe PRIVATE ReSolve::ReSolve)
 enable_testing()
 
 # RESOLVE_DATA is set in test.sh and is the file path the matrix data files used in the testKLU_Rf_FGMRES
-add_test(NAME resolve_consumer COMMAND $<TARGET_FILE:consume.exe> "${RESOLVE_DATA}")
+add_test(NAME resolve_consumer COMMAND $<TARGET_FILE:consume.exe> "-d" "${RESOLVE_DATA}" "-i")

resolve/CMakeLists.txt

@@ -10,6 +10,7 @@ add_subdirectory(utilities)
 
 # C++ files
 set(ReSolve_SRC
+    Common.hpp
     LinSolver.cpp
     LinSolverDirect.cpp
     LinSolverIterative.cpp
@@ -175,4 +176,8 @@ install(TARGETS ReSolve
         LIBRARY DESTINATION lib)
 
 # install include headers
-install(FILES ${ReSolve_HEADER_INSTALL} DESTINATION include/resolve)
+install(DIRECTORY ${CMAKE_SOURCE_DIR}/resolve/ # source directory
+        DESTINATION include/resolve # target directory
+        FILES_MATCHING # install only matched files
+        PATTERN "*.hpp" # select header files
+)

resolve/LinSolverDirectCuSolverGLU.cpp

@@ -151,7 +151,6 @@ namespace ReSolve
         Mshifts[row]++;
       }
     }
-    //Mshifts.~vector(); 
   }
 
   int LinSolverDirectCuSolverGLU::refactorize()
@@ -190,6 +189,7 @@ namespace ReSolve
                                             &r_nrm_inf_,
                                             info_M_,
                                             glu_buffer_);
+    x->setDataUpdated(memory::DEVICE);
     return status_cusolver_; 
   }
 

resolve/LinSolverDirectCuSolverRf.cpp

@@ -1,5 +1,8 @@
+#include <cassert>
+
 #include <resolve/vector/Vector.hpp>
 #include <resolve/matrix/Csr.hpp>
+#include <resolve/matrix/Csc.hpp>
 #include "LinSolverDirectCuSolverRf.hpp"
 
 namespace ReSolve 
@@ -28,16 +31,52 @@ namespace ReSolve
                                        index_type* Q,
                                        vector_type* /* rhs */)
   {
+    assert(A->getSparseFormat() == matrix::Sparse::COMPRESSED_SPARSE_ROW &&
+           "Matrix A has to be in CSR format for cusolverRf input.\n");
+    assert(L->getSparseFormat() == U->getSparseFormat() &&
+           "Matrices L and U have to be in the same format for cusolverRf input.\n");
+
+    int error_sum = 0;
+    this->A_ = A;
+    index_type n = A_->getNumRows();
+
     //remember - P and Q are generally CPU variables
-    // factorization data is stored in the handle. If function is called again, destroy the old handle to get rid of old data. 
+    // factorization data is stored in the handle. 
+    // If function is called again, destroy the old handle to get rid of old data. 
     if (setup_completed_) {
       cusolverRfDestroy(handle_cusolverrf_);
       cusolverRfCreate(&handle_cusolverrf_);
     }
 
-    int error_sum = 0;
-    this->A_ = (matrix::Csr*) A;
-    index_type n = A_->getNumRows();
+    matrix::Csc* L_csc = nullptr;
+    matrix::Csc* U_csc = nullptr;
+    matrix::Csr* L_csr = nullptr;
+    matrix::Csr* U_csr = nullptr;
+
+    matrix::Sparse::SparseFormat matrix_format = L->getSparseFormat();
+    // std::cout << "Matrix format is " << matrix_format << "\n";
+
+    switch (L->getSparseFormat()) {
+      case matrix::Sparse::COMPRESSED_SPARSE_COLUMN:
+        // std::cout << "converting L and U factors from CSC to CSR format ...\n";
+        L_csc = static_cast<matrix::Csc*>(L);
+        U_csc = static_cast<matrix::Csc*>(U);
+        L_csr = new matrix::Csr(L_csc->getNumRows(), L_csc->getNumColumns(), L_csc->getNnz());
+        U_csr = new matrix::Csr(U_csc->getNumRows(), U_csc->getNumColumns(), U_csc->getNnz());
+        csc2csr(L_csc, L_csr);
+        csc2csr(U_csc, U_csr);
+        L_csr->syncData(memory::DEVICE);
+        U_csr->syncData(memory::DEVICE);
+        break;
+      case matrix::Sparse::COMPRESSED_SPARSE_ROW:
+        L_csr = dynamic_cast<matrix::Csr*>(L);
+        U_csr = dynamic_cast<matrix::Csr*>(U);
+        break;
+      default:
+        out::error() << "Matrix type for L and U factors not recognized!\n";
+        out::error() << "Refactorization not completed.\n";
+        return 1;
+    }
 
     if (d_P_ == nullptr){
       mem_.allocateArrayOnDevice(&d_P_, n);
@@ -64,14 +103,14 @@ namespace ReSolve
                                                A_->getRowData(memory::DEVICE),
                                                A_->getColData(memory::DEVICE),
                                                A_->getValues( memory::DEVICE),
-                                               L->getNnz(),
-                                               L->getRowData(memory::DEVICE),
-                                               L->getColData(memory::DEVICE),
-                                               L->getValues( memory::DEVICE),
-                                               U->getNnz(),
-                                               U->getRowData(memory::DEVICE),
-                                               U->getColData(memory::DEVICE),
-                                               U->getValues( memory::DEVICE),
+                                               L_csr->getNnz(),
+                                               L_csr->getRowData(memory::DEVICE),
+                                               L_csr->getColData(memory::DEVICE),
+                                               L_csr->getValues( memory::DEVICE),
+                                               U_csr->getNnz(),
+                                               U_csr->getRowData(memory::DEVICE),
+                                               U_csr->getColData(memory::DEVICE),
+                                               U_csr->getValues( memory::DEVICE),
                                                d_P_,
                                                d_Q_,
                                                handle_cusolverrf_);
@@ -81,16 +120,35 @@ namespace ReSolve
     status_cusolverrf_ = cusolverRfAnalyze(handle_cusolverrf_);
     error_sum += status_cusolverrf_;
 
-    this->A_ = A;
-    //default
-
     const cusolverRfFactorization_t fact_alg =
       CUSOLVERRF_FACTORIZATION_ALG0;  // 0 - default, 1 or 2
     const cusolverRfTriangularSolve_t solve_alg =
       CUSOLVERRF_TRIANGULAR_SOLVE_ALG1;  //  1- default, 2 or 3 // 1 causes error
     this->setAlgorithms(fact_alg, solve_alg);
 
     setup_completed_ = true;
+
+    // Remove temporary objects upon setup completion
+    switch (L->getSparseFormat()) {
+      case matrix::Sparse::COMPRESSED_SPARSE_COLUMN:
+        delete L_csr;
+        delete U_csr;
+        L_csr = nullptr;
+        U_csr = nullptr;
+        L_csc = nullptr;
+        U_csc = nullptr;
+        break;
+      case matrix::Sparse::COMPRESSED_SPARSE_ROW:
+        L_csr = nullptr;
+        U_csr = nullptr;
+        L_csc = nullptr;
+        U_csc = nullptr;
+        break;
+      default:
+        break;
+    }
+    // delete L_csr;
+    // delete U_csr;
 
     return error_sum;
   }
@@ -280,5 +338,94 @@ namespace ReSolve
     params_list_["zero_pivot"]  = ZERO_PIVOT;
     params_list_["pivot_boost"] = PIVOT_BOOST;
   }
+
+  /**
+   * @brief Convert CSC to CSR matrix on the host
+   * 
+   * @authors Slaven Peles <[email protected]>, Daniel Reynolds (SMU), and
+   * David Gardner and Carol Woodward (LLNL)
+   */
+  int LinSolverDirectCuSolverRf::csc2csr(matrix::Csc* A_csc, matrix::Csr* A_csr)
+  {
+    // int error_sum = 0; TODO: Collect error output!
+    assert(A_csc->getNnz() == A_csr->getNnz());
+    assert(A_csc->getNumRows() == A_csr->getNumRows());
+    assert(A_csr->getNumColumns() == A_csc->getNumColumns());
+
+    A_csr->allocateMatrixData(memory::HOST);
+
+    index_type nnz = A_csc->getNnz();
+    index_type n   = A_csc->getNumColumns();
+
+    index_type* rowIdxCsc = A_csc->getRowData(memory::HOST);
+    index_type* colPtrCsc = A_csc->getColData(memory::HOST);
+    real_type*  valuesCsc = A_csc->getValues( memory::HOST);
+
+    index_type* rowPtrCsr = A_csr->getRowData(memory::HOST);
+    index_type* colIdxCsr = A_csr->getColData(memory::HOST);
+    real_type*  valuesCsr = A_csr->getValues( memory::HOST);
+
+    // Set all CSR row pointers to zero
+    for (index_type i = 0; i <= n; ++i) {
+      rowPtrCsr[i] = 0;
+    }
+
+    // Set all CSR values and column indices to zero
+    for (index_type i = 0; i < nnz; ++i) {
+      colIdxCsr[i] = 0;
+      valuesCsr[i] = 0.0;
+    }
+
+    // Compute number of entries per row
+    for (index_type i = 0; i < nnz; ++i) {
+      rowPtrCsr[rowIdxCsc[i]]++;
+    }
+
+    // Compute cumualtive sum of nnz per row
+    for (index_type row = 0, rowsum = 0; row < n; ++row)
+    {
+      // Store value in row pointer to temp
+      index_type temp  = rowPtrCsr[row];
+
+      // Copy cumulative sum to the row pointer
+      rowPtrCsr[row] = rowsum;
+
+      // Update row sum
+      rowsum += temp;
+    }
+    rowPtrCsr[n] = nnz;
+
+    for (index_type col = 0; col < n; ++col)
+    {
+      // Compute positions of column indices and values in CSR matrix and store them there
+      // Overwrites CSR row pointers in the process
+      // adding to them the number of elements in that row
+      for (index_type jj = colPtrCsc[col]; jj < colPtrCsc[col+1]; jj++)
+      {
+          index_type row  = rowIdxCsc[jj];
+          index_type dest = rowPtrCsr[row];
+
+          colIdxCsr[dest] = col;
+          valuesCsr[dest] = valuesCsc[jj];
+
+          rowPtrCsr[row]++;
+      }
+    }
+
+    // Restore CSR row pointer values
+    // All values in rowPtrCsr have shifted by the number of elements in that row
+    // for i>=1: new rowPtrCsr[i] = old rowPtrCsr[i-1] and new rowPtrCsr[0]=0 
+    for (index_type row = 0, last = 0; row <= n; row++)
+    {
+        index_type temp  = rowPtrCsr[row];
+        rowPtrCsr[row] = last;
+        last    = temp;
+    }
+
+    // Mark data on the host as updated
+    A_csr->setUpdated(memory::HOST);
+
+    return 0;
+  }
 
 } // namespace resolve

resolve/LinSolverDirectCuSolverRf.hpp

@@ -18,6 +18,8 @@ namespace ReSolve
   namespace matrix
   {
     class Sparse;
+    class Csr;
+    class Csc;
   }
 
   // Forward declaration of ReSolve handlers workspace
@@ -53,14 +55,15 @@ namespace ReSolve
       bool getCliParamBool(const std::string id) const override;
       int printCliParam(const std::string id) const override;
 
+    private:
+      void initParamList();
+      int csc2csr(matrix::Csc* A_csc, matrix::Csr* A_csr);
+
     private:
       enum ParamaterIDs {ZERO_PIVOT=0, PIVOT_BOOST};
       real_type zero_pivot_{0.0};  ///< The value below which zero pivot is flagged. 
       real_type pivot_boost_{0.0}; ///< The value which is substituted for zero pivot.
 
-    private:
-      void initParamList();
-
       cusolverRfHandle_t handle_cusolverrf_;
       cusolverStatus_t status_cusolverrf_;
 

resolve/LinSolverDirectRocSparseILU0.cpp

@@ -300,10 +300,11 @@ namespace ReSolve
    * are currently redundant code (like an if (true)).
    * In the future, they will be expanded to include more options.
    * 
-   * @param id - string ID for parameter to get.
+   * @param[in] id - string ID for parameter to get.
+   * @param[in] value unused/ignored
    * @return int Value of the int parameter to return.
    */
-  int LinSolverDirectRocSparseILU0::setCliParam(const std::string id, const std::string value)
+  int LinSolverDirectRocSparseILU0::setCliParam(const std::string id, const std::string /* value */)
   {
     switch (getParamId(id))
     {

resolve/SystemSolver.cpp

@@ -388,17 +388,7 @@ namespace ReSolve
       status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_);
     }
     if (refactorizationMethod_ == "cusolverrf") {
-      matrix::Csc* L_csc = dynamic_cast<matrix::Csc*>(L_);
-      matrix::Csc* U_csc = dynamic_cast<matrix::Csc*>(U_);
-      L_csc->syncData(memory::DEVICE);
-      U_csc->syncData(memory::DEVICE);
-      matrix::Csr* L_csr = new matrix::Csr(L_csc->getNumRows(), L_csc->getNumColumns(), L_csc->getNnz());
-      matrix::Csr* U_csr = new matrix::Csr(U_csc->getNumRows(), U_csc->getNumColumns(), U_csc->getNnz());
-      matrixHandler_->csc2csr(L_csc, L_csr, memory::DEVICE);
-      matrixHandler_->csc2csr(U_csc, U_csr, memory::DEVICE);
-      status += refactorizationSolver_->setup(A_, L_csr, U_csr, P_, Q_);
-      delete L_csr;
-      delete U_csr;
+      status += refactorizationSolver_->setup(A_, L_, U_, P_, Q_);
 
       LinSolverDirectCuSolverRf* Rf = dynamic_cast<LinSolverDirectCuSolverRf*>(refactorizationSolver_);
       Rf->setNumericalProperties(1e-14, 1e-1);
@@ -771,15 +761,15 @@ namespace ReSolve
   {
     // Map string input to the Gram-Schmidt variant enum
     GramSchmidt::GSVariant gs_variant;
-    if (variant == "CGS2") {
+    if (variant == "cgs2") {
       gs_variant = GramSchmidt::CGS2;
-    } else if (variant == "MGS") {
+    } else if (variant == "mgs") {
       gs_variant = GramSchmidt::MGS;
-    } else if (variant == "MGS_TWO_SYNC") {
+    } else if (variant == "mgs_two_sync") {
       gs_variant = GramSchmidt::MGS_TWO_SYNC;
-    } else if (variant == "MGS_PM") {
+    } else if (variant == "mgs_pm") {
       gs_variant = GramSchmidt::MGS_PM;
-    } else if (variant == "CGS1") {
+    } else if (variant == "cgs1") {
       gs_variant = GramSchmidt::CGS1;
     } else {
       out::warning() << "Gram-Schmidt variant " << variant << " not recognized.\n";