diff --git a/CHANGELOG.md b/CHANGELOG.md
index cabfe5013..75c3a6701 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/resolve/cuda/cudaVectorKernels.cu b/resolve/cuda/cudaVectorKernels.cu
index 0a879c61d..c604d7e75 100644
--- a/resolve/cuda/cudaVectorKernels.cu
+++ b/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
diff --git a/resolve/cuda/cudaVectorKernels.h b/resolve/cuda/cudaVectorKernels.h
index f8bef0d5b..fccef3447 100644
--- a/resolve/cuda/cudaVectorKernels.h
+++ b/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
diff --git a/resolve/hip/hipVectorKernels.h b/resolve/hip/hipVectorKernels.h
index f5790017d..c7057a9f3 100644
--- a/resolve/hip/hipVectorKernels.h
+++ b/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
diff --git a/resolve/hip/hipVectorKernels.hip b/resolve/hip/hipVectorKernels.hip
index d904c69ca..51ad35c65 100644
--- a/resolve/hip/hipVectorKernels.hip
+++ b/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
diff --git a/resolve/vector/VectorHandler.cpp b/resolve/vector/VectorHandler.cpp
index 86e30ba4b..6b6285d15 100644
--- a/resolve/vector/VectorHandler.cpp
+++ b/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.
    *
diff --git a/resolve/vector/VectorHandler.hpp b/resolve/vector/VectorHandler.hpp
index cbea0c28e..2708a99d1 100644
--- a/resolve/vector/VectorHandler.hpp
+++ b/resolve/vector/VectorHandler.hpp
@@ -60,6 +60,10 @@ namespace ReSolve
               memory::MemorySpace memspace);
 
     int scale(vector::Vector* diag, vector::Vector* vec, memory::MemorySpace memspace);
+    int elementwiseDivide(vector::Vector* divisor, vector::Vector* vec, memory::MemorySpace memspace);
+    int elementwiseMax(vector::Vector* x, vector::Vector* y, memory::MemorySpace memspace);
+
+    int abs(vector::Vector* x, memory::MemorySpace memspace);
 
     /** infNorm:
      * Returns infinity norm of a vector (i.e., entry with max abs value)
diff --git a/resolve/vector/VectorHandlerCpu.cpp b/resolve/vector/VectorHandlerCpu.cpp
index de7ff3869..894e280ef 100644
--- a/resolve/vector/VectorHandlerCpu.cpp
+++ b/resolve/vector/VectorHandlerCpu.cpp
@@ -298,4 +298,73 @@ namespace ReSolve
     return 0;
   }
 
+  /**
+   * @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
+   *
+   * @pre The two vectors must be the same size
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandlerCpu::elementwiseDivide(vector::Vector* divisor, vector::Vector* vec)
+  {
+    real_type* divisor_data = divisor->getData(memory::HOST);
+    real_type* vec_data     = vec->getData(memory::HOST);
+    index_type n            = vec->getSize();
+
+    for (index_type i = 0; i < n; ++i)
+    {
+      vec_data[i] /= divisor_data[i];
+    }
+    vec->setDataUpdated(memory::HOST);
+    return 0;
+  }
+
+  /**
+   * @brief Take the element-wise max of two vectors.
+   * Each element of the output will be the maximum value of the corresponding elements in the input vectors.
+   *
+   * @param[in]     x - First vector
+   * @param[in,out] y - Second vector (output)
+   *
+   * @pre The two vectors must be the same size
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandlerCpu::elementwiseMax(vector::Vector* x, vector::Vector* y)
+  {
+    real_type* x_data = x->getData(memory::HOST);
+    real_type* y_data = y->getData(memory::HOST);
+    index_type n      = y->getSize();
+
+    for (index_type i = 0; i < n; ++i)
+    {
+      y_data[i] = std::max(x_data[i], y_data[i]);
+    }
+    y->setDataUpdated(memory::HOST);
+    return 0;
+  }
+
+  /**
+   * @brief Take the element-wise absolute value of a vector.
+   *
+   * @param[in,out] x - Input and output vector
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandlerCpu::abs(vector::Vector* x)
+  {
+    real_type* x_data = x->getData(memory::HOST);
+    index_type n      = x->getSize();
+
+    for (index_type i = 0; i < n; ++i)
+    {
+      x_data[i] = std::abs(x_data[i]);
+    }
+    x->setDataUpdated(memory::HOST);
+    return 0;
+  }
+
 } // namespace ReSolve
diff --git a/resolve/vector/VectorHandlerCpu.hpp b/resolve/vector/VectorHandlerCpu.hpp
index bf3c2ba83..60103ec3d 100644
--- a/resolve/vector/VectorHandlerCpu.hpp
+++ b/resolve/vector/VectorHandlerCpu.hpp
@@ -58,6 +58,12 @@ namespace ReSolve
 
     virtual int scale(vector::Vector* diag, vector::Vector* vec);
 
+    virtual int elementwiseDivide(vector::Vector* divisor, vector::Vector* vec);
+
+    virtual int elementwiseMax(vector::Vector* x, vector::Vector* y);
+
+    virtual int abs(vector::Vector* x);
+
   private:
     LinAlgWorkspaceCpu* workspace_;
   };
diff --git a/resolve/vector/VectorHandlerCuda.cpp b/resolve/vector/VectorHandlerCuda.cpp
index e60cad347..5ec030615 100644
--- a/resolve/vector/VectorHandlerCuda.cpp
+++ b/resolve/vector/VectorHandlerCuda.cpp
@@ -304,4 +304,62 @@ namespace ReSolve
     return 0;
   }
 
+  /**
+   * @brief Divide a vector's elements by another vector's elements in CUDA
+   *
+   * @param[in]  divisor  - vector of divisors
+   * @param[in, out]  vec - vector to be divided
+   *
+   * @pre The diagonal vector must be of the same size as the vector.
+   * @pre vec is undivided
+   * @post vec is divided
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandlerCuda::elementwiseDivide(vector::Vector* divisor, vector::Vector* vec)
+  {
+    real_type* divisor_data = divisor->getData(memory::DEVICE);
+    real_type* vec_data     = vec->getData(memory::DEVICE);
+    index_type n            = vec->getSize();
+    cuda::elementwiseDivide(n, divisor_data, vec_data);
+    vec->setDataUpdated(memory::DEVICE);
+    return 0;
+  }
+
+  /**
+   * @brief Calculate element-wise maximum between two vectors in CUDA
+   *
+   * @param[in]      x - The first vector
+   * @param[in, out] y - The second vector (result is returned in y)
+   *
+   * @pre The two vectors must be the same size
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandlerCuda::elementwiseMax(vector::Vector* x, vector::Vector* y)
+  {
+    real_type* x_data = x->getData(memory::DEVICE);
+    real_type* y_data = y->getData(memory::DEVICE);
+    index_type n      = y->getSize();
+    cuda::elementwiseMax(n, x_data, y_data);
+    y->setDataUpdated(memory::DEVICE);
+    return 0;
+  }
+
+  /**
+   * @brief Calculate element-wise absolute value of a vector in CUDA
+   *
+   * @param[in, out] x - The vector (result is returned in x)
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandlerCuda::abs(vector::Vector* x)
+  {
+    real_type* x_data = x->getData(memory::DEVICE);
+    index_type n      = x->getSize();
+    cuda::abs(n, x_data);
+    x->setDataUpdated(memory::DEVICE);
+    return 0;
+  }
+
 } // namespace ReSolve
diff --git a/resolve/vector/VectorHandlerCuda.hpp b/resolve/vector/VectorHandlerCuda.hpp
index 5a5b500c7..58dc67278 100644
--- a/resolve/vector/VectorHandlerCuda.hpp
+++ b/resolve/vector/VectorHandlerCuda.hpp
@@ -67,6 +67,35 @@ namespace ReSolve
      */
     virtual int scale(vector::Vector* diag, vector::Vector* vec);
 
+    /**
+     * @brief elementwiseDivide: divides a vector's elements by another's
+     *
+     * @param[in] divisor vector of size n x 1
+     * @param[in,out] vec vector of size n x 1 (this is where the result is stored)
+     *
+     * @return 0 if successful, 1 otherwise
+     */
+    virtual int elementwiseDivide(vector::Vector* divisor, vector::Vector* vec);
+
+    /**
+     * @brief elementwiseMax: calculate the element-wise maximum of two vectors
+     *
+     * @param[in]     x vector of size n x 1
+     * @param[in,out] y vector of size n x 1 (this is where the result is stored)
+     *
+     * @return 0 if successful, 1 otherwise
+     */
+    virtual int elementwiseMax(vector::Vector* x, vector::Vector* y);
+
+    /**
+     * @brief abs: calculate the element-wise absolute value of a vector
+     *
+     * @param[in,out] x vector of size n x 1 (this is where the result is stored)
+     *
+     * @return 0 if successful, 1 otherwise
+     */
+    virtual int abs(vector::Vector* x);
+
   private:
     MemoryHandler        mem_; ///< Device memory manager object
     LinAlgWorkspaceCUDA* workspace_;
diff --git a/resolve/vector/VectorHandlerHip.cpp b/resolve/vector/VectorHandlerHip.cpp
index f3930d8ac..60a09505e 100644
--- a/resolve/vector/VectorHandlerHip.cpp
+++ b/resolve/vector/VectorHandlerHip.cpp
@@ -290,4 +290,62 @@ namespace ReSolve
     return 0;
   }
 
+  /**
+   * @brief Divide a vector's elements by another vector's elements in HIP
+   *
+   * @param[in]  divisor  - vector of divisors
+   * @param[in, out]  vec - vector to be divided
+   *
+   * @pre The diagonal vector must be of the same size as the vector.
+   * @pre vec is undivided
+   * @post vec is divided
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandlerHip::elementwiseDivide(vector::Vector* divisor, vector::Vector* vec)
+  {
+    real_type* divisor_data = divisor->getData(memory::DEVICE);
+    real_type* vec_data     = vec->getData(memory::DEVICE);
+    index_type n            = vec->getSize();
+    hip::elementwiseDivide(n, divisor_data, vec_data);
+    vec->setDataUpdated(memory::DEVICE);
+    return 0;
+  }
+
+  /**
+   * @brief Calculate element-wise maximum between two vectors in HIP
+   *
+   * @param[in]      x - The first vector
+   * @param[in, out] y - The second vector (result is returned in y)
+   *
+   * @pre The two vectors must be the same size
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandlerHip::elementwiseMax(vector::Vector* x, vector::Vector* y)
+  {
+    real_type* x_data = x->getData(memory::DEVICE);
+    real_type* y_data = y->getData(memory::DEVICE);
+    index_type n      = y->getSize();
+    hip::elementwiseMax(n, x_data, y_data);
+    y->setDataUpdated(memory::DEVICE);
+    return 0;
+  }
+
+  /**
+   * @brief Calculate element-wise absolute value of a vector in HIP
+   *
+   * @param[in, out] x - The vector (result is returned in x)
+   *
+   * @return 0 if successful, 1 otherwise
+   */
+  int VectorHandlerHip::abs(vector::Vector* x)
+  {
+    real_type* x_data = x->getData(memory::DEVICE);
+    index_type n      = x->getSize();
+    hip::abs(n, x_data);
+    x->setDataUpdated(memory::DEVICE);
+    return 0;
+  }
+
 } // namespace ReSolve
diff --git a/resolve/vector/VectorHandlerHip.hpp b/resolve/vector/VectorHandlerHip.hpp
index a1601dddb..e083518a6 100644
--- a/resolve/vector/VectorHandlerHip.hpp
+++ b/resolve/vector/VectorHandlerHip.hpp
@@ -67,6 +67,35 @@ namespace ReSolve
      */
     virtual int scale(vector::Vector* diag, vector::Vector* vec);
 
+    /**
+     * @brief elementwiseDivide: divides a vector's elements by another's
+     *
+     * @param[in] divisor vector of size n x 1
+     * @param[in,out] vec vector of size n x 1 (this is where the result is stored)
+     *
+     * @return 0 if successful, 1 otherwise
+     */
+    virtual int elementwiseDivide(vector::Vector* divisor, vector::Vector* vec);
+
+    /**
+     * @brief elementwiseMax: calculate the element-wise maximum of two vectors
+     *
+     * @param[in]     x vector of size n x 1
+     * @param[in,out] y vector of size n x 1 (this is where the result is stored)
+     *
+     * @return 0 if successful, 1 otherwise
+     */
+    virtual int elementwiseMax(vector::Vector* x, vector::Vector* y);
+
+    /**
+     * @brief abs: calculate the element-wise absolute value of a vector
+     *
+     * @param[in,out] x vector of size n x 1 (this is where the result is stored)
+     *
+     * @return 0 if successful, 1 otherwise
+     */
+    virtual int abs(vector::Vector* x);
+
   private:
     LinAlgWorkspaceHIP* workspace_;
     MemoryHandler       mem_; ///< Device memory manager object
diff --git a/resolve/vector/VectorHandlerImpl.hpp b/resolve/vector/VectorHandlerImpl.hpp
index 9b6c0a073..7290b08f6 100644
--- a/resolve/vector/VectorHandlerImpl.hpp
+++ b/resolve/vector/VectorHandlerImpl.hpp
@@ -45,6 +45,15 @@ namespace ReSolve
     // Scale a vector by a diagonal matrix
     virtual int scale(vector::Vector* diag, vector::Vector* vec) = 0;
 
+    // Divide the elements of a vector by the elements of another vector
+    virtual int elementwiseDivide(vector::Vector* divisor, vector::Vector* vec) = 0;
+
+    // Compute element-wise max of two vectors
+    virtual int elementwiseMax(vector::Vector* x, vector::Vector* y) = 0;
+
+    // Compute element-wise absolute value of a vector
+    virtual int abs(vector::Vector* x) = 0;
+
     /** gemv:
      * if `transpose = N` (no), `x = beta*x +  alpha*V*y`,
      * where `x` is `[n x 1]`, `V` is `[n x k]` and `y` is `[k x 1]`.
diff --git a/tests/unit/vector/VectorHandlerTests.hpp b/tests/unit/vector/VectorHandlerTests.hpp
index 7372ec976..52cc8bd05 100644
--- a/tests/unit/vector/VectorHandlerTests.hpp
+++ b/tests/unit/vector/VectorHandlerTests.hpp
@@ -280,6 +280,140 @@ namespace ReSolve
         return status.report(__func__);
       }
 
+      TestOutcome elementwiseDivide(index_type N)
+      {
+        TestStatus status;
+
+        vector::Vector divisor(N);
+        vector::Vector vec(N);
+
+        // divisor[i] = i + 1, vec[i] = 3.0
+        // expected result vec[i] = 3.0 / (i + 1)
+        divisor.allocate(memspace_);
+        vec.allocate(memspace_);
+
+        vec.setToConst(3.0, memspace_);
+
+        auto divisor_data = std::unique_ptr<real_type[]>(new real_type[N]);
+        for (size_t i = 0; i < static_cast<size_t>(N); ++i)
+        {
+          divisor_data[i] = (real_type) (i + 1);
+        }
+        divisor.copyDataFrom(divisor_data.get(), memory::HOST, memspace_);
+
+        handler_.elementwiseDivide(&divisor, &vec, memspace_);
+
+        if (memspace_ == memory::DEVICE)
+        {
+          vec.syncData(memory::HOST);
+        }
+
+        for (index_type i = 0; i < N; ++i)
+        {
+          if (!isEqual(vec.getData(memory::HOST)[i], (real_type) 3.0 / (i + 1)))
+          {
+            std::cout << "Solution vector element vec[" << i << "] = " << vec.getData(memory::HOST)[i]
+                      << ", expected: " << (real_type) 3.0 / (i + 1) << "\n";
+            status *= false;
+            break;
+          }
+        }
+
+        return status.report(__func__);
+      }
+
+      TestOutcome elementwiseMax(index_type N)
+      {
+        TestStatus status;
+
+        vector::Vector x(N);
+        vector::Vector y(N);
+
+        x.allocate(memspace_);
+        y.allocate(memspace_);
+
+        auto x_data = std::unique_ptr<real_type[]>(new real_type[N]);
+        auto y_data = std::unique_ptr<real_type[]>(new real_type[N]);
+        for (size_t i = 0; i < static_cast<size_t>(N); ++i)
+        {
+          if (i % 3 == 0)
+          {
+            x_data[i] = (real_type) (i + 1);
+            y_data[i] = (real_type) i * .5;
+          }
+          else
+          {
+            x_data[i] = -(real_type) (i + 1);
+            y_data[i] = (real_type) (i + 1);
+          }
+        }
+        x.copyDataFrom(x_data.get(), memory::HOST, memspace_);
+        y.copyDataFrom(y_data.get(), memory::HOST, memspace_);
+
+        handler_.elementwiseMax(&x, &y, memspace_);
+
+        if (memspace_ == memory::DEVICE)
+        {
+          y.syncData(memory::HOST);
+        }
+
+        for (index_type i = 0; i < N; ++i)
+        {
+          if (!isEqual(y.getData(memory::HOST)[i], (real_type) (i + 1)))
+          {
+            std::cout << "Solution vector element y[" << i << "] = " << y.getData(memory::HOST)[i]
+                      << ", expected: " << (real_type) (i + 1) << "\n";
+            status *= false;
+            break;
+          }
+        }
+
+        return status.report(__func__);
+      }
+
+      TestOutcome abs(index_type N)
+      {
+        TestStatus status;
+
+        vector::Vector x(N);
+
+        x.allocate(memspace_);
+
+        auto x_data = std::unique_ptr<real_type[]>(new real_type[N]);
+        for (size_t i = 0; i < static_cast<size_t>(N); ++i)
+        {
+          if (i % 3 == 0)
+          {
+            x_data[i] = -(real_type) i;
+          }
+          else
+          {
+            x_data[i] = (real_type) i;
+          }
+        }
+        x.copyDataFrom(x_data.get(), memory::HOST, memspace_);
+
+        handler_.abs(&x, memspace_);
+
+        if (memspace_ == memory::DEVICE)
+        {
+          x.syncData(memory::HOST);
+        }
+
+        for (index_type i = 0; i < N; ++i)
+        {
+          if (!isEqual(x.getData(memory::HOST)[i], (real_type) i))
+          {
+            std::cout << "Solution vector element y[" << i << "] = " << x.getData(memory::HOST)[i]
+                      << ", expected: " << (real_type) i << "\n";
+            status *= false;
+            break;
+          }
+        }
+
+        return status.report(__func__);
+      }
+
     private:
       ReSolve::VectorHandler&      handler_;
       ReSolve::memory::MemorySpace memspace_{memory::HOST};
diff --git a/tests/unit/vector/runVectorHandlerTests.cpp b/tests/unit/vector/runVectorHandlerTests.cpp
index 74fae7e41..c9f575eda 100644
--- a/tests/unit/vector/runVectorHandlerTests.cpp
+++ b/tests/unit/vector/runVectorHandlerTests.cpp
@@ -25,6 +25,9 @@ int main(int, char**)
     result += test.massAxpy(100, 10);
     result += test.massDot(100, 10);
     result += test.scale(100);
+    result += test.elementwiseDivide(100);
+    result += test.elementwiseMax(100);
+    result += test.abs(100);
 
     std::cout << "\n";
   }
@@ -48,6 +51,9 @@ int main(int, char**)
     result += test.massDot(1000, 30);
     result += test.infNorm(1000);
     result += test.scale(1000);
+    result += test.elementwiseDivide(1000);
+    result += test.elementwiseMax(1000);
+    result += test.abs(1000);
 
     std::cout << "\n";
   }
@@ -72,6 +78,9 @@ int main(int, char**)
     result += test.massDot(1000, 30);
     result += test.infNorm(1000);
     result += test.scale(1000);
+    result += test.elementwiseDivide(1000);
+    result += test.elementwiseMax(1000);
+    result += test.abs(1000);
 
     std::cout << "\n";
   }