Add 3 new element-wise vector operations

CHANGELOG.md

@@ -18,6 +18,8 @@
 
 - Made Vector::copyDataTo able to copy from device to host and vice versa
 
+- Added `diagSolve`, `max`, and `abs` vector operations.
+
 ## Changes to Re::Solve in release 0.99.2
 
 ### Major Features

resolve/cuda/cudaVectorKernels.cu

@@ -12,6 +12,8 @@
 #include <resolve/cuda/cudaKernels.h>
 #include <resolve/cuda/cudaVectorKernels.h>
 
+#include "resolve/Common.hpp"
+
 namespace ReSolve
 {
   namespace cuda
@@ -81,8 +83,83 @@ namespace ReSolve
         }
       }
 
+      /**
+       * @brief Multiplies vector by an inverse of a diagonal matrix.
+       *
+       * @param[in]  n       - size of the vectors
+       * @param[in]  d_val   - diagonal matrix values
+       * @param[in, out] vec - vector to be divided. Changes in place.
+       *
+       * @todo Decide how to allow user to configure grid and block sizes.
+       */
+      __global__ void diagSolve(index_type       n,
+                                const real_type* d_val,
+                                real_type*       vec)
+      {
+        // Get the index of the element to be processed
+        index_type idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+        // Check if the index is within bounds
+        if (idx < n)
+        {
+          // Divide the vector element by the corresponding diag value
+          vec[idx] /= d_val[idx];
+        }
+      }
+
+      /**
+       * @brief Computes the element-wise max of two vectors.
+       *
+       * @param[in]  n   - size of the vectors
+       * @param[in]  x   - First vector values
+       * @param[in]  y   - Second vector values
+       * @param[out] out - Output vector values
+       *
+       * @todo Decide how to allow user to configure grid and block sizes.
+       */
+      __global__ void max(index_type       n,
+                          const real_type* x,
+                          const real_type* y,
+                          real_type*       out)
+      {
+        // Get the index of the element to be processed
+        index_type idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+        // Check if the index is within bounds
+        if (idx < n)
+        {
+          // Compute maximum of elements
+          out[idx] = fmax(x[idx], y[idx]);
+        }
+      }
+
+      /**
+       * @brief Computes the element-wise absolute value of a vector.
+       *
+       * @param[in]  n   - size of the vector
+       * @param[in]  in  - Vector input
+       * @param[out] out - Vector output
+       *
+       * @todo Decide how to allow user to configure grid and block sizes.
+       */
+      __global__ void abs(index_type       n,
+                          const real_type* in,
+                          real_type*       out)
+      {
+        // Get the index of the element to be processed
+        index_type idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+        // Check if the index is within bounds
+        if (idx < n)
+        {
+          // Compute absolute value of element
+          out[idx] = fabs(in[idx]);
+        }
+      }
     } // namespace kernels
 
+    constexpr index_type block_size = 256;
+
     void setArrayConst(index_type n, real_type val, real_type* arr)
     {
       index_type num_blocks;
@@ -118,5 +195,61 @@ namespace ReSolve
       // Launch the kernel
       kernels::scale<<<num_blocks, block_size>>>(n, diag, vec);
     }
+
+    /**
+     * @brief Wrapper that divides a vector's elements by another vector's elements
+     *
+     * @param[in]  n       - size of the vectors
+     * @param[in]  diag - diag values
+     * @param[in, out] vec - vector to be divided. Changes in place.
+     *
+     * @todo Decide how to allow user to configure grid and block sizes.
+     */
+    void diagSolve(index_type       n,
+                   const real_type* diag,
+                   real_type*       vec)
+    {
+      int num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::diagSolve<<<num_blocks, block_size>>>(n, diag, vec);
+    }
+
+    /**
+     * @brief Wrapper that computes the element-wise max of two vectors.
+     *
+     * @param[in]  n   - size of the vectors
+     * @param[in]  x   - First vector values
+     * @param[in]  y   - Second vector values
+     * @param[out] out - Output vector values
+     *
+     * @todo Decide how to allow user to configure grid and block sizes.
+     */
+    void max(index_type       n,
+             const real_type* x,
+             const real_type* y,
+             real_type*       out)
+    {
+      int num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::max<<<num_blocks, block_size>>>(n, x, y, out);
+    }
+
+    /**
+     * @brief Wrapper that computes the element-wise absolute value of a vector.
+     *
+     * @param[in]  n   - size of the vector
+     * @param[in]  in  - Vector input
+     * @param[out] out - Vector output
+     *
+     * @todo Decide how to allow user to configure grid and block sizes.
+     */
+    void abs(index_type       n,
+             const real_type* in,
+             real_type*       out)
+    {
+      int num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::abs<<<num_blocks, block_size>>>(n, in, out);
+    }
   } // namespace cuda
 } // namespace ReSolve

resolve/cuda/cudaVectorKernels.h

@@ -18,5 +18,8 @@ namespace ReSolve
     void setArrayConst(index_type n, real_type val, real_type* arr);
     void addConst(index_type n, real_type val, real_type* arr);
     void scale(index_type n, const real_type* diag, real_type* vec);
+    void diagSolve(index_type n, const real_type* diag, real_type* vec);
+    void max(index_type n, const real_type* x, const real_type* y, real_type* out);
+    void abs(index_type n, const real_type* in, real_type* out);
   } // namespace cuda
 } // namespace ReSolve

resolve/hip/hipVectorKernels.h

@@ -18,5 +18,8 @@ namespace ReSolve
     void setArrayConst(index_type n, real_type val, real_type* arr);
     void addConst(index_type n, real_type val, real_type* arr);
     void scale(index_type n, const real_type* diag, real_type* vec);
+    void diagSolve(index_type n, const real_type* diag, real_type* vec);
+    void max(index_type n, const real_type* x, const real_type* y, real_type* out);
+    void abs(index_type n, const real_type* in, real_type* out);
   } // namespace hip
 } // namespace ReSolve

resolve/hip/hipVectorKernels.hip

@@ -69,8 +69,84 @@ namespace ReSolve {
           vec[idx] *= diag[idx];
         }
       }
+
+      /**
+       * @brief Multiplies vector by an inverse of a diagonal matrix.
+       *
+       * @param[in]  n       - size of the vectors
+       * @param[in]  d_val   - diagonal matrix values
+       * @param[in, out] vec - vector to be divided. Changes in place.
+       *
+       * @todo Decide how to allow user to configure grid and block sizes.
+       */
+      __global__ void diagSolve(index_type       n,
+                                const real_type* d_val,
+                                real_type*       vec)
+      {
+        // Get the index of the element to be processed
+        index_type idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+        // Check if the index is within bounds
+        if (idx < n)
+        {
+          // Divide the vector element by the corresponding diag value
+          vec[idx] /= d_val[idx];
+        }
+      }
+
+      /**
+       * @brief Computes the element-wise max of two vectors.
+       *
+       * @param[in]  n   - size of the vectors
+       * @param[in]  x   - First vector values
+       * @param[in]  y   - Second vector values
+       * @param[out] out - Output vector values
+       *
+       * @todo Decide how to allow user to configure grid and block sizes.
+       */
+      __global__ void max(index_type       n,
+                          const real_type* x,
+                          const real_type* y,
+                          real_type*       out)
+      {
+        // Get the index of the element to be processed
+        index_type idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+        // Check if the index is within bounds
+        if (idx < n)
+        {
+          // Compute maximum of elements
+          out[idx] = fmax(x[idx], y[idx]);
+        }
+      }
+
+      /**
+       * @brief Computes the element-wise absolute value of a vector.
+       *
+       * @param[in]  n   - size of the vector
+       * @param[in]  in  - Vector input
+       * @param[out] out - Vector output
+       *
+       * @todo Decide how to allow user to configure grid and block sizes.
+       */
+      __global__ void abs(index_type       n,
+                          const real_type* in,
+                          real_type*       out)
+      {
+        // Get the index of the element to be processed
+        index_type idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+        // Check if the index is within bounds
+        if (idx < n)
+        {
+          // Compute absolute value of element
+          out[idx] = fabs(in[idx]);
+        }
+      }
     } // namespace kernels
 
+    constexpr index_type block_size = 256;
+
     void setArrayConst(index_type n, real_type c, real_type* v)
     {
       index_type num_blocks;
@@ -106,5 +182,61 @@ namespace ReSolve {
       // Launch the kernel
       kernels::scale<<<num_blocks, block_size>>>(n, diag, vec);
     }
+
+    /**
+     * @brief Wrapper that divides a vector's elements by another vector's elements
+     *
+     * @param[in]  n       - size of the vectors
+     * @param[in]  diag - diag values
+     * @param[in, out] vec - vector to be divided. Changes in place.
+     *
+     * @todo Decide how to allow user to configure grid and block sizes.
+     */
+    void diagSolve(index_type       n,
+                           const real_type* diag,
+                           real_type*       vec)
+    {
+      int num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::diagSolve<<<num_blocks, block_size>>>(n, diag, vec);
+    }
+
+    /**
+     * @brief Wrapper that computes the element-wise max of two vectors.
+     *
+     * @param[in]  n   - size of the vectors
+     * @param[in]  x   - First vector values
+     * @param[in]  y   - Second vector values
+     * @param[out] out - Output vector values
+     *
+     * @todo Decide how to allow user to configure grid and block sizes.
+     */
+    void max(index_type       n,
+             const real_type* x,
+             const real_type* y,
+             real_type*       out)
+    {
+      int num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::max<<<num_blocks, block_size>>>(n, x, y, out);
+    }
+
+    /**
+     * @brief Wrapper that computes the element-wise absolute value of a vector.
+     *
+     * @param[in]  n   - size of the vector
+     * @param[in]  in  - Vector input
+     * @param[out] out - Vector output
+     *
+     * @todo Decide how to allow user to configure grid and block sizes.
+     */
+    void abs(index_type       n,
+             const real_type* in,
+             real_type*       out)
+    {
+      int num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::abs<<<num_blocks, block_size>>>(n, in, out);
+    }
   } // namespace hip
 } // namespace ReSolve

resolve/hykkt/permutation/Permutation.cpp

@@ -247,8 +247,8 @@ namespace ReSolve
                                real_type*      old_val,
                                real_type*      new_val)
     {
-      index_type  length;
-      index_type* apply_perm_;
+      index_type  length{0};
+      index_type* apply_perm_{nullptr};
       switch (permutation)
       {
       case PERM_V:

resolve/hykkt/ruiz/CMakeLists.txt

@@ -18,7 +18,9 @@ set(RUIZ_HEADER_INSTALL
 # Build shared library ReSolve::hykkt
 add_library(resolve_hykkt_ruiz SHARED ${RUIZ_SRC})
 
-target_link_libraries(resolve_hykkt_ruiz PRIVATE resolve_logger resolve_vector)
+target_link_libraries(
+  resolve_hykkt_ruiz PRIVATE resolve_logger resolve_vector resolve_matrix
+)
 
 # Link to CUDA ReSolve backend if CUDA is support enabled
 if(RESOLVE_USE_CUDA)

resolve/hykkt/spgemm/CMakeLists.txt

@@ -21,6 +21,7 @@ add_library(resolve_hykkt_spgemm SHARED ${SPGEMM_SRC})
 
 target_link_libraries(
   resolve_hykkt_spgemm PRIVATE resolve_logger ${suitesparse_cholmod}
+                               resolve_matrix resolve_vector
 )
 target_include_directories(
   resolve_hykkt_spgemm PUBLIC ${SUITESPARSE_INCLUDE_DIR}