Implement A := c*A+I and A := c*A+B methods in MatrixHandler

CHANGELOG.md

@@ -4,6 +4,10 @@
 
 - Added cmake-format
 
+- The MatrixHandler::scaleAddI function was added to calculate A := cA + I
+
+- The MatrixHandler::scaleAddB function was added to calculate A := cA + B
+
 ## Changes to Re::Solve in release 0.99.2
 
 ### Major Features

resolve/LinSolverDirectRocSolverRf.cpp

@@ -1,5 +1,7 @@
 #include "LinSolverDirectRocSolverRf.hpp"
 
+#include <cassert>
+
 #include <resolve/hip/hipKernels.h>
 
 #include <resolve/Profiling.hpp>

resolve/MemoryUtils.hpp

@@ -38,7 +38,7 @@ namespace ReSolve
    *
    * @brief Provides basic memory allocation, free and copy functions.
    *
-   * This class provedes abstractions for memory management functiosn for
+   * This class provedes abstractions for memory management functions for
    * different GPU programming models.
    *
    * @tparam Policy - Memory management policy (vendor specific)

resolve/matrix/MatrixHandler.cpp

@@ -351,6 +351,57 @@ namespace ReSolve
     }
   }
 
+  /**
+   * @brief Multiply csr matrix by a constant and add I.
+   * @param[in,out] A - Sparse CSR matrix
+   * @param[in] alpha - scalar parameter
+   * @param[in] memspace - Device where the operation is computed
+   * @return 0 if successful, 1 otherwise
+   */
+  int MatrixHandler::scaleAddI(matrix::Csr* A, real_type alpha, memory::MemorySpace memspace)
+  {
+    assert(A->getSparseFormat() == matrix::Sparse::COMPRESSED_SPARSE_ROW && "Matrix is assumed to be in CSR format.\n");
+    assert((A->getValues(memspace) != nullptr || A->getNnz() == 0) && "Matrix values are null!\n");
+    assert(A->getNumRows() == A->getNumColumns() && "Matrix must be square.\n");
+    using namespace ReSolve::memory;
+    switch (memspace)
+    {
+    case HOST:
+      return cpuImpl_->scaleAddI(A, alpha);
+      break;
+    case DEVICE:
+      return devImpl_->scaleAddI(A, alpha);
+      break;
+    }
+    return 1;
+  }
+
+  /**
+   * @brief Multiply csr matrix by a constant and add B.
+   * @param[in,out] A - Sparse CSR matrix
+   * @param[in] alpha - scalar parameter
+   * @param[in] B - Sparse CSR matrix
+   * @param[in] memspace - Device where the operation is computed
+   * @return 0 if successful, 1 otherwise
+   */
+  int MatrixHandler::scaleAddB(matrix::Csr* A, real_type alpha, matrix::Csr* B, memory::MemorySpace memspace)
+  {
+    assert(A->getSparseFormat() == matrix::Sparse::COMPRESSED_SPARSE_ROW && "Matrix is assumed to be in CSR format.\n");
+    assert(A->getValues(memspace) != nullptr && "Matrix values are null!\n");
+    assert(A->getNumRows() == A->getNumColumns() && "Matrix must be square.\n");
+    using namespace ReSolve::memory;
+    switch (memspace)
+    {
+    case HOST:
+      return cpuImpl_->scaleAddB(A, alpha, B);
+      break;
+    case DEVICE:
+      return devImpl_->scaleAddB(A, alpha, B);
+      break;
+    }
+    return 1;
+  }
+
   /**
    * @brief If CUDA support is enabled in the handler.
    *

resolve/matrix/MatrixHandler.hpp

@@ -62,6 +62,10 @@ namespace ReSolve
 
     void addConst(matrix::Sparse* A, real_type alpha, memory::MemorySpace memspace);
 
+    int scaleAddI(matrix::Csr* A, real_type alpha, memory::MemorySpace memspace);
+
+    int scaleAddB(matrix::Csr* A, real_type alpha, matrix::Csr* B, memory::MemorySpace memspace);
+
     /// Should compute vec_result := alpha*A*vec_x + beta*vec_result, but at least on cpu alpha and beta are flipped
     int  matvec(matrix::Sparse*     A,
                 vector_type*        vec_x,