Matrix handler now checks in the background if the matrix is CSR. (#190)

* Matrix handler checks if matrix is CSR. No other formats are supported by the handler. * Sparse matrix base class has enum that indicates sparse format. * Remove string input parameter specifying matrix format in matvec function.
ORNL · Sep 26, 2024 · bd7b1c1 · bd7b1c1
1 parent a1864d6
commit bd7b1c1
Show file tree

Hide file tree

Showing 39 changed files with 370 additions and 271 deletions.
diff --git a/examples/r_KLU_GLU.cpp b/examples/r_KLU_GLU.cpp
@@ -137,7 +137,7 @@ int main(int argc, char *argv[])
     vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
     real_type bnorm = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE));
     matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
-    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
 
 
     matrix_handler->matrixInfNorm(A, &norm_A, ReSolve::memory::DEVICE); 

diff --git a/examples/r_KLU_GLU_matrix_values_update.cpp b/examples/r_KLU_GLU_matrix_values_update.cpp
@@ -148,7 +148,7 @@ int main(int argc, char *argv[])
 
 
       matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
-      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
 
       std::cout << "\t 2-Norm of the residual: " 
                 << std::scientific << std::setprecision(16) 

diff --git a/examples/r_KLU_KLU.cpp b/examples/r_KLU_KLU.cpp
@@ -129,7 +129,7 @@ int main(int argc, char *argv[])
 
     matrix_handler->setValuesChanged(true, ReSolve::memory::HOST);
 
-    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", ReSolve::memory::HOST); 
+    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::HOST); 
     norm_r = vector_handler->infNorm(vec_r, ReSolve::memory::HOST);
 
     std::cout << "\t2-Norm of the residual: " 

diff --git a/examples/r_KLU_KLU_standalone.cpp b/examples/r_KLU_KLU_standalone.cpp
@@ -92,7 +92,7 @@ int main(int argc, char *argv[])
 
   matrix_handler->setValuesChanged(true, ReSolve::memory::HOST);
 
-  matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", ReSolve::memory::HOST); 
+  matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::HOST); 
 
   std::cout << "\t 2-Norm of the residual: " 
             << std::scientific << std::setprecision(16) 

diff --git a/examples/r_KLU_cusolverrf_redo_factorization.cpp b/examples/r_KLU_cusolverrf_redo_factorization.cpp
@@ -155,7 +155,7 @@ int main(int argc, char *argv[] )
 
     matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
 
-    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
     res_nrm = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE));
     b_nrm   = sqrt(vector_handler->dot(vec_rhs, vec_rhs, ReSolve::memory::DEVICE));
     std::cout << "\t2-Norm of the residual: " 
@@ -180,7 +180,7 @@ int main(int argc, char *argv[] )
 
       matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
 
-      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
       res_nrm = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE));
 
       std::cout <<"\t New residual norm: "

diff --git a/examples/r_KLU_rf.cpp b/examples/r_KLU_rf.cpp
@@ -148,7 +148,7 @@ int main(int argc, char *argv[] )
 
     matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
 
-    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
 
     std::cout << "\t 2-Norm of the residual: " 
               << std::scientific << std::setprecision(16) 

diff --git a/examples/r_KLU_rf_FGMRES.cpp b/examples/r_KLU_rf_FGMRES.cpp
@@ -128,7 +128,7 @@ int main(int argc, char *argv[])
       norm_b = vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE);
       norm_b = sqrt(norm_b);
       matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
-      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
 
       matrix_handler->matrixInfNorm(A, &norm_A, ReSolve::memory::DEVICE); 
       norm_x = vector_handler->infNorm(vec_x, ReSolve::memory::DEVICE);
@@ -174,7 +174,7 @@ int main(int argc, char *argv[])
       FGMRES->resetMatrix(A);
       FGMRES->setupPreconditioner("LU", Rf);
 
-      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
 
       matrix_handler->matrixInfNorm(A, &norm_A, ReSolve::memory::DEVICE); 
       norm_x = vector_handler->infNorm(vec_x, ReSolve::memory::DEVICE);

diff --git a/examples/r_KLU_rf_FGMRES_reuse_factorization.cpp b/examples/r_KLU_rf_FGMRES_reuse_factorization.cpp
@@ -130,7 +130,7 @@ int main(int argc, char *argv[])
       norm_b = vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE);
       norm_b = sqrt(norm_b);
       matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
-      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
       std::cout << "\t 2-Norm of the residual : " 
                 << std::scientific << std::setprecision(16) 
                 << sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE))/norm_b << "\n";
@@ -181,7 +181,7 @@ int main(int argc, char *argv[])
       matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
       FGMRES->resetMatrix(A);
 
-      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
 
       std::cout << "\t 2-Norm of the residual (before IR): " 
                 << std::scientific << std::setprecision(16) 

diff --git a/examples/r_KLU_rocSolverRf_FGMRES.cpp b/examples/r_KLU_rocSolverRf_FGMRES.cpp
@@ -140,7 +140,7 @@ int main(int argc, char *argv[])
       norm_b = vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE);
       norm_b = sqrt(norm_b);
       matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
-      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
       std::cout << "\t2-Norm of the residual: "
         << std::scientific << std::setprecision(16) 
         << sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE))/norm_b << "\n";
@@ -176,7 +176,7 @@ int main(int argc, char *argv[])
       FGMRES->resetMatrix(A);
       FGMRES->setupPreconditioner("LU", Rf);
 
-      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+      matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
       real_type rnrm = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE));
       std::cout << "\t 2-Norm of the residual (before IR): " 
         << std::scientific << std::setprecision(16) 

diff --git a/examples/r_KLU_rocsolverrf.cpp b/examples/r_KLU_rocsolverrf.cpp
@@ -138,7 +138,7 @@ int main(int argc, char *argv[] )
     vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
     real_type bnorm = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE));
     matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
-    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
 
     std::cout << "\t 2-Norm of the residual: " 
               << std::scientific << std::setprecision(16) 

diff --git a/examples/r_KLU_rocsolverrf_redo_factorization.cpp b/examples/r_KLU_rocsolverrf_redo_factorization.cpp
@@ -140,7 +140,7 @@ int main(int argc, char *argv[] )
 
     matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
 
-    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
     res_nrm = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE));
     b_nrm = sqrt(vector_handler->dot(vec_rhs, vec_rhs, ReSolve::memory::DEVICE));
     std::cout << "\t 2-Norm of the residual: " 
@@ -163,7 +163,7 @@ int main(int argc, char *argv[] )
 
          matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
 
-         matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", ReSolve::memory::DEVICE); 
+         matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); 
          res_nrm = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE));
 
          std::cout<<"\t New residual norm: "

diff --git a/examples/r_SysSolver.cpp b/examples/r_SysSolver.cpp
@@ -129,7 +129,7 @@ int main(int argc, char *argv[])
 
     matrix_handler->setValuesChanged(true, ReSolve::memory::HOST);
 
-    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", ReSolve::memory::HOST); 
+    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::HOST); 
 
     std::cout << "\t 2-Norm of the residual: " 
               << std::scientific << std::setprecision(16) 

diff --git a/resolve/LinSolverIterativeFGMRES.cpp b/resolve/LinSolverIterativeFGMRES.cpp
@@ -124,7 +124,7 @@ namespace ReSolve
     vec_V_->setToZero(memspace_);
 
     rhs->deepCopyVectorData(vec_V_->getData(memspace_), 0, memspace_);  
-    matrix_handler_->matvec(A_, x, vec_V_, &MINUSONE, &ONE, "csr", memspace_); 
+    matrix_handler_->matvec(A_, x, vec_V_, &MINUSONE, &ONE, memspace_); 
     rnorm = 0.0;
     bnorm = vector_handler_->dot(rhs, rhs, memspace_);
     rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_);
@@ -189,7 +189,7 @@ namespace ReSolve
 
         vec_v->setData( vec_V_->getVectorData(i + 1, memspace_), memspace_);
 
-        matrix_handler_->matvec(A_, vec_z, vec_v, &ONE, &ZERO,"csr", memspace_); 
+        matrix_handler_->matvec(A_, vec_z, vec_v, &ONE, &ZERO, memspace_); 
 
         // orthogonalize V[i+1], form a column of h_H_
 
@@ -274,7 +274,7 @@ namespace ReSolve
       }
 
       rhs->deepCopyVectorData(vec_V_->getData(memspace_), 0, memspace_);  
-      matrix_handler_->matvec(A_, x, vec_V_, &MINUSONE, &ONE,"csr", memspace_); 
+      matrix_handler_->matvec(A_, x, vec_V_, &MINUSONE, &ONE, memspace_); 
       rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_);
       // rnorm = ||V_1||
       rnorm = std::sqrt(rnorm);

diff --git a/resolve/LinSolverIterativeRandFGMRES.cpp b/resolve/LinSolverIterativeRandFGMRES.cpp
@@ -154,7 +154,7 @@ namespace ReSolve
     vec_V_->setToZero(memspace_);
 
     rhs->deepCopyVectorData(vec_V_->getData(memspace_), 0, memspace_);  
-    matrix_handler_->matvec(A_, x, vec_V_, &MINUSONE, &ONE, "csr", memspace_); 
+    matrix_handler_->matvec(A_, x, vec_V_, &MINUSONE, &ONE, memspace_); 
 
     vec_v->setData(vec_V_->getVectorData(0, memspace_), memspace_);
     vec_s->setData(vec_S_->getVectorData(0, memspace_), memspace_);
@@ -228,7 +228,7 @@ namespace ReSolve
         // V_{i+1}=A*Z_i
         vec_v->setData(vec_V_->getVectorData(i + 1, memspace_), memspace_);
 
-        matrix_handler_->matvec(A_, vec_z, vec_v, &ONE, &ZERO, "csr", memspace_); 
+        matrix_handler_->matvec(A_, vec_z, vec_v, &ONE, &ZERO, memspace_); 
 
         // orthogonalize V[i+1], form a column of h_H_
         // this is where it differs from normal solver GS
@@ -337,7 +337,7 @@ namespace ReSolve
       }
 
       rhs->deepCopyVectorData(vec_V_->getData(memspace_), 0, memspace_);  
-      matrix_handler_->matvec(A_, x, vec_V_, &MINUSONE, &ONE,"csr", memspace_); 
+      matrix_handler_->matvec(A_, x, vec_V_, &MINUSONE, &ONE, memspace_); 
       if (outer_flag) {
 
         sketching_handler_->reset();

diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp
@@ -682,7 +682,7 @@ namespace ReSolve
       return -1.0;
     }
     matrixHandler_->setValuesChanged(true, ms);
-    matrixHandler_->matvec(A_, x, resVector_, &ONE, &MINUSONE, "csr", ms);
+    matrixHandler_->matvec(A_, x, resVector_, &ONE, &MINUSONE, ms);
     resnorm = std::sqrt(vectorHandler_->dot(resVector_, resVector_, ms));
     return resnorm/norm_b;
   }
@@ -711,7 +711,7 @@ namespace ReSolve
       return -1.0;
     }
     matrixHandler_->setValuesChanged(true, ms);
-    matrixHandler_->matvec(A_, x, resVector_, &ONE, &MINUSONE, "csr", ms);
+    matrixHandler_->matvec(A_, x, resVector_, &ONE, &MINUSONE, ms);
     resnorm = vectorHandler_->infNorm(resVector_, ms);
     norm_x  = vectorHandler_->infNorm(x, ms);
     matrixHandler_->matrixInfNorm(A_, &norm_A, ms);

diff --git a/resolve/matrix/Coo.cpp b/resolve/matrix/Coo.cpp
@@ -12,10 +12,12 @@ namespace ReSolve
 
   matrix::Coo::Coo()
   {
+    sparse_format_ = TRIPLET;
   }
 
   matrix::Coo::Coo(index_type n, index_type m, index_type nnz) : Sparse(n, m, nnz)
   {
+    sparse_format_ = TRIPLET;
   }
 
   matrix::Coo::Coo(index_type n, 
@@ -24,6 +26,7 @@ namespace ReSolve
                    bool symmetric,
                    bool expanded) : Sparse(n, m, nnz, symmetric, expanded)
   {
+    sparse_format_ = TRIPLET;
   }
 
   /**
@@ -41,7 +44,8 @@ namespace ReSolve
                    memory::MemorySpace memspaceDst)
     : Sparse(n, m, nnz, symmetric, expanded)
   {
-    using namespace memory;
+    sparse_format_ = TRIPLET;
+
     int control = -1;
     if ((memspaceSrc == memory::HOST)   && (memspaceDst == memory::HOST))  { control = 0;}
     if ((memspaceSrc == memory::HOST)   && (memspaceDst == memory::DEVICE)){ control = 1;}

diff --git a/resolve/matrix/Csc.cpp b/resolve/matrix/Csc.cpp
@@ -10,6 +10,7 @@ namespace ReSolve
 
   matrix::Csc::Csc()
   {
+    sparse_format_ = COMPRESSED_SPARSE_COLUMN;
   }
 
   matrix::Csc::Csc(index_type n, index_type m, index_type nnz) : Sparse(n, m, nnz)
@@ -22,6 +23,7 @@ namespace ReSolve
                        bool symmetric,
                        bool expanded) : Sparse(n, m, nnz, symmetric, expanded)
   {
+    sparse_format_ = COMPRESSED_SPARSE_COLUMN;
   }
 
   matrix::Csc::~Csc()

diff --git a/resolve/matrix/Csr.cpp b/resolve/matrix/Csr.cpp
@@ -13,10 +13,12 @@ namespace ReSolve
 
   matrix::Csr::Csr()
   {
+    sparse_format_ = COMPRESSED_SPARSE_ROW;
   }
 
   matrix::Csr::Csr(index_type n, index_type m, index_type nnz) : Sparse(n, m, nnz)
   {
+    sparse_format_ = COMPRESSED_SPARSE_ROW;
   }
 
   matrix::Csr::Csr(index_type n, 
@@ -25,6 +27,7 @@ namespace ReSolve
                    bool symmetric,
                    bool expanded) : Sparse(n, m, nnz, symmetric, expanded)
   {
+    sparse_format_ = COMPRESSED_SPARSE_ROW;
   }
 
   /**
@@ -53,7 +56,8 @@ namespace ReSolve
                    memory::MemorySpace memspaceDst)
     : Sparse(n, m, nnz, symmetric, expanded)
   {
-    using namespace memory;
+    sparse_format_ = COMPRESSED_SPARSE_ROW;
+
     int control = -1;
     if ((memspaceSrc == memory::HOST)   && (memspaceDst == memory::HOST))  { control = 0;}
     if ((memspaceSrc == memory::HOST)   && (memspaceDst == memory::DEVICE)){ control = 1;}

diff --git a/resolve/matrix/MatrixHandler.cpp b/resolve/matrix/MatrixHandler.cpp
@@ -116,7 +116,6 @@ namespace ReSolve {
    * @param[out] vec_result - Vector where the result is stored
    * @param[in]  alpha - scalar parameter
    * @param[in]  beta  - scalar parameter
-   * @param[in]  matrixFormat - Only CSR format is supported at this time
    * @param[in]  memspace     - Device where the product is computed
    * @return result := alpha * A * x + beta * result
    */
@@ -125,16 +124,15 @@ namespace ReSolve {
                             vector_type* vec_result, 
                             const real_type* alpha, 
                             const real_type* beta,
-                            std::string matrixFormat, 
                             memory::MemorySpace memspace)
   {
     using namespace ReSolve::memory;
     switch (memspace) {
       case HOST:
-        return cpuImpl_->matvec(A, vec_x, vec_result, alpha, beta, matrixFormat);
+        return cpuImpl_->matvec(A, vec_x, vec_result, alpha, beta);
         break;
       case DEVICE:
-        return devImpl_->matvec(A, vec_x, vec_result, alpha, beta, matrixFormat);
+        return devImpl_->matvec(A, vec_x, vec_result, alpha, beta);
         break;
     }
     return 1;

diff --git a/resolve/matrix/MatrixHandler.hpp b/resolve/matrix/MatrixHandler.hpp
@@ -60,7 +60,6 @@ namespace ReSolve {
                  vector_type* vec_result,
                  const real_type* alpha,
                  const real_type* beta,
-                 std::string matrix_type,
                  memory::MemorySpace memspace);
       int matrixInfNorm(matrix::Sparse *A, real_type* norm, memory::MemorySpace memspace);
       void setValuesChanged(bool toWhat, memory::MemorySpace memspace);