Add 3 new element-wise vector operations

CHANGELOG.md

@@ -12,6 +12,8 @@
 
 - Removed unnecessary full facotorization in the examples and made the input 1 based.
 
+- Added `elementwiseDivide`, `elementwiseMax`, and `abs` vector operations.
+
 ## Changes to Re::Solve in release 0.99.2
 
 ### Major Features

resolve/cuda/cudaVectorKernels.cu

@@ -81,6 +81,75 @@ namespace ReSolve
         }
       }
 
+      /**
+       * @brief Divides a vector's elements by another vector's elements
+       *
+       * @param[in]  n       - size of the vectors
+       * @param[in]  d_val   - divisor 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 elementwiseDivide(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 divisor value
+          vec[idx] /= d_val[idx];
+        }
+      }
+
+      /**
+       * @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, out] y - Second vector values. Changes in place.
+       *
+       * @todo Decide how to allow user to configure grid and block sizes.
+       */
+      __global__ void elementwiseMax(index_type       n,
+                                     const real_type* x,
+                                     real_type*       y)
+      {
+        // 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
+          y[idx] = max(x[idx], y[idx]);
+        }
+      }
+
+      /**
+       * @brief Computes the element-wise absolute value of a vector.
+       *
+       * @param[in]      n - size of the vector
+       * @param[in, out] x - Vector values. Changes in place.
+       *
+       * @todo Decide how to allow user to configure grid and block sizes.
+       */
+      __global__ void abs(index_type n,
+                          real_type* x)
+      {
+        // 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
+          x[idx] = fabs(x[idx]);
+        }
+      }
     } // namespace kernels
 
     void setArrayConst(index_type n, real_type val, real_type* arr)
@@ -118,5 +187,63 @@ 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]  divisor - divisor 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 elementwiseDivide(index_type       n,
+                           const real_type* divisor,
+                           real_type*       vec)
+    {
+      // Define block size and number of blocks
+      const int block_size = 256;
+      int       num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::elementwiseDivide<<<num_blocks, block_size>>>(n, divisor, 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, out] y - Second vector values. Changes in place.
+     *
+     * @todo Decide how to allow user to configure grid and block sizes.
+     */
+    void elementwiseMax(index_type       n,
+                        const real_type* x,
+                        real_type*       y)
+    {
+      // Define block size and number of blocks
+      const int block_size = 256;
+      int       num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::elementwiseMax<<<num_blocks, block_size>>>(n, x, y);
+    }
+
+    /**
+     * @brief Wrapper that computes the element-wise absolute value of a vector.
+     *
+     * @param[in]      n - size of the vector
+     * @param[in, out] x - Vector values. Changes in place.
+     *
+     * @todo Decide how to allow user to configure grid and block sizes.
+     */
+    void abs(index_type n,
+             real_type* x)
+    {
+      // Define block size and number of blocks
+      const int block_size = 256;
+      int       num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::abs<<<num_blocks, block_size>>>(n, x);
+    }
   } // 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 elementwiseDivide(index_type n, const real_type* divisor, real_type* vec);
+    void elementwiseMax(index_type n, const real_type* x, real_type* y);
+    void abs(index_type n, real_type* x);
   } // 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 elementwiseDivide(index_type n, const real_type* divisor, real_type* vec);
+    void elementwiseMax(index_type n, const real_type* x, real_type* y);
+    void abs(index_type n, real_type* x);
   } // namespace hip
 } // namespace ReSolve

resolve/hip/hipVectorKernels.hip

@@ -69,6 +69,76 @@ namespace ReSolve {
           vec[idx] *= diag[idx];
         }
       }
+
+      /**
+       * @brief Divides a vector's elements by another vector's elements
+       *
+       * @param[in]  n       - size of the vectors
+       * @param[in]  d_val   - divisor 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 elementwiseDivide(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 divisor value
+          vec[idx] /= d_val[idx];
+        }
+      }
+
+      /**
+       * @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, out] y - Second vector values. Changes in place.
+       *
+       * @todo Decide how to allow user to configure grid and block sizes.
+       */
+      __global__ void elementwiseMax(index_type       n,
+                                     const real_type* x,
+                                     real_type*       y)
+      {
+        // 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
+          y[idx] = max(x[idx], y[idx]);
+        }
+      }
+
+      /**
+       * @brief Computes the element-wise absolute value of a vector.
+       *
+       * @param[in]      n - size of the vector
+       * @param[in, out] x - Vector values. Changes in place.
+       *
+       * @todo Decide how to allow user to configure grid and block sizes.
+       */
+      __global__ void abs(index_type n,
+                          real_type* x)
+      {
+        // 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
+          x[idx] = fabs(x[idx]);
+        }
+      }
     } // namespace kernels
 
     void setArrayConst(index_type n, real_type c, real_type* v)
@@ -106,5 +176,63 @@ 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]  divisor - divisor 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 elementwiseDivide(index_type       n,
+                           const real_type* divisor,
+                           real_type*       vec)
+    {
+      // Define block size and number of blocks
+      const int block_size = 256;
+      int       num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::elementwiseDivide<<<num_blocks, block_size>>>(n, divisor, 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, out] y - Second vector values. Changes in place.
+     *
+     * @todo Decide how to allow user to configure grid and block sizes.
+     */
+    void elementwiseMax(index_type       n,
+                        const real_type* x,
+                        real_type*       y)
+    {
+      // Define block size and number of blocks
+      const int block_size = 256;
+      int       num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::elementwiseMax<<<num_blocks, block_size>>>(n, x, y);
+    }
+
+    /**
+     * @brief Wrapper that computes the element-wise absolute value of a vector.
+     *
+     * @param[in]      n - size of the vector
+     * @param[in, out] x - Vector values. Changes in place.
+     *
+     * @todo Decide how to allow user to configure grid and block sizes.
+     */
+    void abs(index_type n,
+             real_type* x)
+    {
+      // Define block size and number of blocks
+      const int block_size = 256;
+      int       num_blocks = (n + block_size - 1) / block_size;
+      // Launch the kernel
+      kernels::abs<<<num_blocks, block_size>>>(n, x);
+    }
   } // namespace hip
 } // namespace ReSolve

resolve/vector/VectorHandler.cpp

@@ -307,6 +307,88 @@ namespace ReSolve
     return 1;
   }
 
+  /**
+   * @brief Divide the elements of a vector by the elements of another vector
+   *
+   * @param[in] divisor - vector divisor
+   * @param[in,out] vec - vector to be divided
+   * @param[in] memspace - Device where the operation is computed
+   *
+   * @pre The two vectors must be the same size
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandler::elementwiseDivide(vector::Vector* divisor, vector::Vector* vec, memory::MemorySpace memspace)
+  {
+    assert(divisor->getSize() == vec->getSize() && "Diagonal vector must be of the same size as the vector.");
+    assert(divisor->getData(memspace) != nullptr && "Diagonal vector data is null!\n");
+    assert(vec->getData(memspace) != nullptr && "Vector data is null!\n");
+    using namespace ReSolve::memory;
+    switch (memspace)
+    {
+    case HOST:
+      return cpuImpl_->elementwiseDivide(divisor, vec);
+      break;
+    case DEVICE:
+      return devImpl_->elementwiseDivide(divisor, vec);
+      break;
+    }
+    return 1;
+  }
+
+  /**
+   * @brief Takes the element-wise max between two vectors.
+   *
+   * @param[in] x        - The first vector
+   * @param[in,out] y    - The second vector (result is returned in y)
+   * @param[in] memspace - Device where the operation is computed
+   *
+   * @pre The two vectors must be the same size
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandler::elementwiseMax(vector::Vector* x, vector::Vector* y, memory::MemorySpace memspace)
+  {
+    assert(x->getSize() == y->getSize() && "Vectors must be the same size.");
+    assert(x->getData(memspace) != nullptr && "Vector x data is null!\n");
+    assert(y->getData(memspace) != nullptr && "Vector y data is null!\n");
+    using namespace ReSolve::memory;
+    switch (memspace)
+    {
+    case HOST:
+      return cpuImpl_->elementwiseMax(x, y);
+      break;
+    case DEVICE:
+      return devImpl_->elementwiseMax(x, y);
+      break;
+    }
+    return 1;
+  }
+
+  /**
+   * @brief Computes the element-wise absolute value of a vector.
+   *
+   * @param[in,out] x    - Input and output vector
+   * @param[in] memspace - Device where the operation is computed
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandler::abs(vector::Vector* x, memory::MemorySpace memspace)
+  {
+    assert(x->getData(memspace) != nullptr && "Vector data is null!\n");
+    using namespace ReSolve::memory;
+    switch (memspace)
+    {
+    case HOST:
+      return cpuImpl_->abs(x);
+      break;
+    case DEVICE:
+      return devImpl_->abs(x);
+      break;
+    }
+    return 1;
+  }
+
   /**
    * @brief If CUDA support is enabled in the handler.
    *