From 4b88a60794d1e0714e85811a2a3fb4bb30da4c0a Mon Sep 17 00:00:00 2001 From: pelesh Date: Wed, 8 Jul 2026 13:07:13 -0400 Subject: [PATCH 01/10] Add VectorHandler class to GK --- GridKit/LinearAlgebra/Vector/CMakeLists.txt | 4 +- .../LinearAlgebra/Vector/VectorHandler.cpp | 353 ++++++++++++ .../LinearAlgebra/Vector/VectorHandler.hpp | 88 +++ .../LinearAlgebra/Vector/VectorHandlerCpu.cpp | 388 +++++++++++++ .../LinearAlgebra/Vector/VectorHandlerCpu.hpp | 74 +++ .../LinearAlgebra/Vector/CMakeLists.txt | 8 + .../Vector/VectorHandlerTests.hpp | 509 ++++++++++++++++++ .../Vector/runVectorHandlerTests.cpp | 44 ++ 8 files changed, 1466 insertions(+), 2 deletions(-) create mode 100644 GridKit/LinearAlgebra/Vector/VectorHandler.cpp create mode 100644 GridKit/LinearAlgebra/Vector/VectorHandler.hpp create mode 100644 GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp create mode 100644 GridKit/LinearAlgebra/Vector/VectorHandlerCpu.hpp create mode 100644 tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp create mode 100644 tests/UnitTests/LinearAlgebra/Vector/runVectorHandlerTests.cpp diff --git a/GridKit/LinearAlgebra/Vector/CMakeLists.txt b/GridKit/LinearAlgebra/Vector/CMakeLists.txt index 314eed26e..4bc319e9c 100644 --- a/GridKit/LinearAlgebra/Vector/CMakeLists.txt +++ b/GridKit/LinearAlgebra/Vector/CMakeLists.txt @@ -1,7 +1,7 @@ gridkit_add_library( dense_vector - SOURCES Vector.cpp - HEADERS Vector.hpp + SOURCES Vector.cpp VectorHandler.cpp VectorHandlerCpu.cpp + HEADERS Vector.hpp VectorHandler.hpp VectorHandlerCpu.hpp LINK_LIBRARIES GridKit::utilities_logger GridKit::cpu_backend) # # If no GPU is enabled link to dummy device backend diff --git a/GridKit/LinearAlgebra/Vector/VectorHandler.cpp b/GridKit/LinearAlgebra/Vector/VectorHandler.cpp new file mode 100644 index 000000000..817ee4df8 --- /dev/null +++ b/GridKit/LinearAlgebra/Vector/VectorHandler.cpp @@ -0,0 +1,353 @@ +#include "VectorHandler.hpp" + +#include +#include + +#include +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + using out = GridKit::Utilities::Logger; + + /** + * @brief dot product of two vectors i.e, a = x^Ty + * + * @param[in] x The first vector + * @param[in] y The second vector + * @param[in] memspace Memory space the operation is computed in (HOST or DEVICE) + * + * @return dot product of _x_ and _y_ + */ + template + ScalarT VectorHandler::dot(Vector* x, Vector* y, memory::MemorySpace memspace) + { + switch (memspace) + { + case memory::HOST: + return cpuImpl_.dot(x, y); + case memory::DEVICE: + out::error() << "VectorHandler::dot - DEVICE memory space not yet supported\n"; + return static_cast(NAN); + } + return static_cast(NAN); + } + + /** + * @brief scale a vector by a constant i.e, x = alpha*x where alpha is a constant + * + * @param[in] alpha The constant + * @param[in,out] x The vector + * @param[in] memspace Memory space the operation is computed in (HOST or DEVICE) + */ + template + void VectorHandler::scal(const ScalarT alpha, Vector* x, memory::MemorySpace memspace) + { + switch (memspace) + { + case memory::HOST: + cpuImpl_.scal(alpha, x); + break; + case memory::DEVICE: + out::error() << "VectorHandler::scal - DEVICE memory space not yet supported\n"; + break; + } + } + + /** + * @brief compute infinity norm of a vector (i.e., find an entry with largest absolute value) + * + * @param[in] x The vector + * @param[in] memspace Memory space the operation is computed in (HOST or DEVICE) + * + * @return infinity norm of _x_ + */ + template + ScalarT VectorHandler::amax(Vector* x, memory::MemorySpace memspace) + { + switch (memspace) + { + case memory::HOST: + return cpuImpl_.amax(x); + case memory::DEVICE: + out::error() << "VectorHandler::amax - DEVICE memory space not yet supported\n"; + return static_cast(NAN); + } + return static_cast(NAN); + } + + /** + * @brief axpy i.e, y = alpha*x+y where alpha is a constant + * + * @param[in] alpha The constant + * @param[in] x The first vector + * @param[in,out] y The second vector (result is returned in y) + * @param[in] memspace Memory space the operation is computed in (HOST or DEVICE) + */ + template + void VectorHandler::axpy(const ScalarT alpha, + Vector* x, + Vector* y, + memory::MemorySpace memspace) + { + switch (memspace) + { + case memory::HOST: + cpuImpl_.axpy(alpha, x, y); + break; + case memory::DEVICE: + out::error() << "VectorHandler::axpy - DEVICE memory space not yet supported\n"; + break; + } + } + + /** + * @brief gemv computes dense matrix-vector product. + * + * 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]`. + * If `transpose = T` (yes), `x := beta*x + alpha*V^T*y`, + * where `x` is `[k x 1]`, `V` is `[n x k]` and `y` is `[n x 1]`. + * + * @param[in] transpose - yes (T) or no (N) + * @param[in] k - Number of columns in (non-transposed) matrix to use + * @param[in] alpha - Constant scalar + * @param[in] beta - Constant scalar + * @param[in] V - Multivector containing the matrix, organized columnwise + * @param[in] y - Vector, k x 1 if N and n x 1 if T + * @param[in,out] x - Vector, n x 1 if N and k x 1 if T + * @param[in] memspace - Memory space the operation is computed in (HOST or DEVICE) + * + * @note Parameter k is not the total number of columns in V but the number + * of columns to use in matrix-vector product. + */ + template + void VectorHandler::gemv(char transpose, + IdxT k, + const ScalarT alpha, + const ScalarT beta, + Vector* V, + Vector* y, + Vector* x, + memory::MemorySpace memspace) + { + switch (memspace) + { + case memory::HOST: + cpuImpl_.gemv(transpose, k, alpha, beta, V, y, x); + break; + case memory::DEVICE: + out::error() << "VectorHandler::gemv - DEVICE memory space not yet supported\n"; + break; + } + } + + /** + * @brief Multivector axpy: y := y - \sum_i alpha_i x_i + * + * @param[in] size number of elements in y + * @param[in] alpha vector size k x 1 + * @param[in] k number of vectors in the multivector x + * @param[in] x (multi)vector [size x k] + * @param[in,out] y vector size size x 1 (this is where the result is stored) + * @param[in] memspace Memory space the operation is computed in (HOST or DEVICE) + * + * @pre _k_ > 0, _size_ > 0, _size_ = x->getSize() + */ + template + void VectorHandler::axpyMulti(IdxT size, + Vector* alpha, + IdxT k, + Vector* x, + Vector* y, + memory::MemorySpace memspace) + { + assert(y->getSize() == x->getSize() && "Sizes of x and y must match!\n"); + assert(alpha->getSize() == k && "Size of alpha must match k!\n"); + + switch (memspace) + { + case memory::HOST: + cpuImpl_.axpyMulti(size, alpha, k, x, y); + break; + case memory::DEVICE: + out::error() << "VectorHandler::axpyMulti - DEVICE memory space not yet supported\n"; + break; + } + } + + /** + * @brief Multivector dot product, i.e V^T x + * + * Computes V^T x with k vectors from multivector V. Result is stored in `res`. + * + * @param[in] size - Number of elements in a single vector in V + * @param[in] V - Multivector; k vectors of size n x 1 each + * @param[in] k - Number of vectors in V to use + * @param[in] x - Multivector; 2 vectors of size n x 1 each + * @param[out] res - Multivector; 2 vectors size k x 1 each + * @param[in] memspace - Memory space the operation is computed in (HOST or DEVICE) + * + * @pre _size_ > 0, _k_ > 0, size = x->getSize(). + * @pre _res_ needs to be allocated to k x 2 size. + */ + template + void VectorHandler::dot2Multi(IdxT size, + Vector* V, + IdxT k, + Vector* x, + Vector* res, + memory::MemorySpace memspace) + { + assert(x->getSize() == V->getSize() && "Sizes of V and x do not match!\n"); + assert(res->getSize() == k && "Size of `res` must match k!\n"); + + switch (memspace) + { + case memory::HOST: + cpuImpl_.dot2Multi(size, V, k, x, res); + break; + case memory::DEVICE: + out::error() << "VectorHandler::dot2Multi - DEVICE memory space not yet supported\n"; + break; + } + } + + /** + * @brief Scale a vector by a diagonal matrix + * + * @param[in] diag - vector representing the diagonal matrix + * @param[in,out] vec - vector to be scaled + * @param[in] memspace - Memory space the operation is computed in (HOST or DEVICE) + * + * @pre The diagonal vector must be of the same size as the vector. + */ + template + void VectorHandler::scal(Vector* diag, Vector* vec, memory::MemorySpace memspace) + { + assert(diag->getSize() == vec->getSize() && "Diagonal vector must be of the same size as the vector."); + + switch (memspace) + { + case memory::HOST: + cpuImpl_.scal(diag, vec); + break; + case memory::DEVICE: + out::error() << "VectorHandler::scal - DEVICE memory space not yet supported\n"; + break; + } + } + + /** + * @brief Scale a vector by a diagonal matrix represented by a contiguous subvector of an input vector + * + * @param[in] diag - vector representing the diagonal matrix + * @param[in,out] vec - vector to be scaled + * @param[in] diag_offset - the index of diag where the diagonal matrix begins (inclusive) + * @param[in] memspace - Memory space the operation is computed in (HOST or DEVICE) + */ + template + void VectorHandler::scal(Vector* diag, + Vector* vec, + IdxT diag_offset, + memory::MemorySpace memspace) + { + switch (memspace) + { + case memory::HOST: + cpuImpl_.scal(diag, vec, diag_offset); + break; + case memory::DEVICE: + out::error() << "VectorHandler::scal - DEVICE memory space not yet supported\n"; + break; + } + } + + /** + * @brief Multiplies vector by an inverse of a diagonal matrix. + * + * @param[in] diag - diagonal matrix stored in a vector object + * @param[in,out] vec - vector to be divided + * @param[in] memspace - Memory space the operation is computed in (HOST or DEVICE) + * + * @pre The two vectors must be the same size + * + * @return 0 if successful, 1 otherwise + */ + template + int VectorHandler::diagSolve(Vector* diag, Vector* vec, memory::MemorySpace memspace) + { + assert(diag->getSize() == vec->getSize() && "Diagonal vector must be of the same size as the vector."); + + switch (memspace) + { + case memory::HOST: + return cpuImpl_.diagSolve(diag, vec); + case memory::DEVICE: + out::error() << "VectorHandler::diagSolve - DEVICE memory space not yet supported\n"; + return 1; + } + return 1; + } + + /** + * @brief Takes the element-wise max between two vectors. + * + * @param[in] x - The first vector + * @param[in] y - The second vector + * @param[out] out - The output vector + * @param[in] memspace - Memory space the operation is computed in (HOST or DEVICE) + * + * @pre The two vectors must be the same size + * + * @return 0 if successful, 1 otherwise + */ + template + int VectorHandler::max(Vector* x, Vector* y, Vector* out, memory::MemorySpace memspace) + { + assert(x->getSize() == y->getSize() && "Vectors must be the same size."); + assert(x->getSize() == out->getSize() && "Vectors must be the same size."); + + switch (memspace) + { + case memory::HOST: + return cpuImpl_.max(x, y, out); + case memory::DEVICE: + GridKit::LinearAlgebra::out::error() << "VectorHandler::max - DEVICE memory space not yet supported\n"; + return 1; + } + return 1; + } + + /** + * @brief Computes the element-wise absolute value of a vector. + * + * @param[in] in - Input vector + * @param[out] out - Output vector + * @param[in] memspace - Memory space the operation is computed in (HOST or DEVICE) + * + * @return 0 if successful, 1 otherwise + */ + template + int VectorHandler::abs(Vector* in, Vector* out, memory::MemorySpace memspace) + { + assert(in->getSize() == out->getSize() && "Vector sizes do not match!\n"); + + switch (memspace) + { + case memory::HOST: + return cpuImpl_.abs(in, out); + case memory::DEVICE: + GridKit::LinearAlgebra::out::error() << "VectorHandler::abs - DEVICE memory space not yet supported\n"; + return 1; + } + return 1; + } + + // Available template instantiations + template class VectorHandler; + + } // namespace LinearAlgebra +} // namespace GridKit diff --git a/GridKit/LinearAlgebra/Vector/VectorHandler.hpp b/GridKit/LinearAlgebra/Vector/VectorHandler.hpp new file mode 100644 index 000000000..b718458c2 --- /dev/null +++ b/GridKit/LinearAlgebra/Vector/VectorHandler.hpp @@ -0,0 +1,88 @@ +#pragma once + +#include +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + template + class Vector; + + /** + * @brief Dispatches dense (multi)vector operations to the implementation + * appropriate for the requested memory space. + * + * Vector kernels operating on HOST memory are always available. Support + * for DEVICE memory space is provided only when GridKit is built with GPU + * support enabled and a device implementation for the requested kernel is + * available. + * + * @author Slaven Peles + */ + template + class VectorHandler + { + public: + VectorHandler() = default; + ~VectorHandler() = default; + + // y := alpha * x + y + void axpy(const ScalarT alpha, + Vector* x, + Vector* y, + memory::MemorySpace memspace); + + // Dot product of two vectors + ScalarT dot(Vector* x, Vector* y, memory::MemorySpace memspace); + + // Scale vector by scalar + void scal(const ScalarT alpha, Vector* x, memory::MemorySpace memspace); + + // Scale vector by diagonal matrix represented as a vector (i.e., vec = diag*vec) + void scal(Vector* diag, Vector* vec, memory::MemorySpace memspace); + void scal(Vector* diag, + Vector* vec, + IdxT diag_offset, + memory::MemorySpace memspace); + + // axpy for multivectors + void axpyMulti(IdxT size, + Vector* alpha, + IdxT k, + Vector* x, + Vector* y, + memory::MemorySpace memspace); + + // Dot product of two vectors with a multivector V + void dot2Multi(IdxT size, + Vector* V, + IdxT k, + Vector* x, + Vector* res, + memory::MemorySpace memspace); + + // Dense matrix-vector product. + void gemv(char transpose, + IdxT k, // number of vectors from multivector V to use + const ScalarT alpha, + const ScalarT beta, + Vector* V, + Vector* y, + Vector* x, + memory::MemorySpace memspace); + + int diagSolve(Vector* diag, Vector* vec, memory::MemorySpace memspace); + int max(Vector* x, Vector* y, Vector* out, memory::MemorySpace memspace); + int abs(Vector* in, Vector* out, memory::MemorySpace memspace); + + // Vector infinity norm + ScalarT amax(Vector* x, memory::MemorySpace memspace); + + private: + VectorHandlerCpu cpuImpl_; + }; // class VectorHandler + + } // namespace LinearAlgebra +} // namespace GridKit diff --git a/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp new file mode 100644 index 000000000..b46c977c4 --- /dev/null +++ b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp @@ -0,0 +1,388 @@ +#include "VectorHandlerCpu.hpp" + +#include +#include + +#include +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + using out = GridKit::Utilities::Logger; + + /** + * @brief dot product of two vectors i.e, a = x^Ty + * + * @param[in] x The first vector + * @param[in] y The second vector + * + * @return dot product of _x_ and _y_ + */ + template + ScalarT VectorHandlerCpu::dot(Vector* x, Vector* y) + { + const ScalarT* x_data = x->getData(memory::HOST); + const ScalarT* y_data = y->getData(memory::HOST); + ScalarT sum = 0.0; + ScalarT c = 0.0; + // Kahan summation to reduce round-off error accumulation + for (IdxT i = 0; i < x->getSize(); ++i) + { + ScalarT y = (x_data[i] * y_data[i]) - c; + ScalarT t = sum + y; + c = (t - sum) - y; + sum = t; + } + return sum; + } + + /** + * @brief scale a vector by a constant i.e, x = alpha*x where alpha is a constant + * + * @param[in] alpha The constant + * @param[in,out] x The vector + */ + template + void VectorHandlerCpu::scal(const ScalarT alpha, Vector* x) + { + ScalarT* x_data = x->getData(memory::HOST); + + for (IdxT i = 0; i < x->getSize(); ++i) + { + x_data[i] *= alpha; + } + x->setDataUpdated(memory::HOST); + } + + /** + * @brief compute infinity norm of a vector (i.e., find an entry with largest absolute value) + * + * @param[in] x vector + * + * @return infinity norm of _x_ + */ + template + ScalarT VectorHandlerCpu::amax(Vector* x) + { + const ScalarT* x_data = x->getData(memory::HOST); + + ScalarT vecmax = std::abs(x_data[0]); + ScalarT v; + for (IdxT i = 1; i < x->getSize(); ++i) + { + v = std::abs(x_data[i]); + if (v > vecmax) + { + vecmax = v; + } + } + return vecmax; + } + + /** + * @brief axpy i.e, y = alpha*x+y where alpha is a constant + * + * @param[in] alpha The constant + * @param[in] x The first vector + * @param[in,out] y The second vector (result is returned in y) + */ + template + void VectorHandlerCpu::axpy(const ScalarT alpha, Vector* x, Vector* y) + { + ScalarT* x_data = x->getData(memory::HOST); + ScalarT* y_data = y->getData(memory::HOST); + for (IdxT i = 0; i < x->getSize(); ++i) + { + y_data[i] = alpha * x_data[i] + y_data[i]; + } + y->setDataUpdated(memory::HOST); + } + + /** + * @brief gemv computes matrix-vector product where both matrix and vectors are dense. + * i.e., x = beta*x + alpha*V*y + * + * @param[in] transpose - transposed = 'T' or not 'N' + * @param[in] k Number of columns in (non-transposed) matrix to use + * @param[in] alpha Constant scalar + * @param[in] beta Constant scalar + * @param[in] V Multivector containing the matrix, organized columnwise + * @param[in] y Vector, k x 1 if N and n x 1 if T + * @param[in,out] x Vector, n x 1 if N and k x 1 if T + * + * @note Parameter k is not the total number of columns in V but the number + * of columns to use in matrix-vector product. + * + * @pre Number of columns in V >= k + * @pre If transpose = N, size of y must equal k. If transpose = T, size of + * x must equal k. + */ + template + void VectorHandlerCpu::gemv(char transpose, + IdxT k, + const ScalarT alpha, + const ScalarT beta, + Vector* V, + Vector* y, + Vector* x) + { + // x = beta*x + alpha*V*y OR x = beta*x + alpha*V^Ty + const ScalarT* V_data = V->getData(memory::HOST); + const ScalarT* y_data = y->getData(memory::HOST); + ScalarT* x_data = x->getData(memory::HOST); + const IdxT n = V->getSize(); + + IdxT i, j; + ScalarT sum; + switch (transpose) + { + case 'T': + assert((V->getSize() == y->getSize()) + && "gemv: Size mismatch! Size of V does not match size of y."); + for (i = 0; i < k; ++i) + { + sum = beta * x_data[i]; + ScalarT c = 0.0; + for (j = 0; j < n; ++j) + { + ScalarT y = (alpha * V_data[i * n + j] * y_data[j]) - c; + ScalarT t = sum + y; + c = (t - sum) - y; + sum = t; + } + x_data[i] = sum; + } + break; + case 'N': + assert((V->getSize() == x->getSize()) + && "gemv: Size mismatch! Size of V does not match size of x."); + for (i = 0; i < n; ++i) + { + sum = beta * x_data[i]; + ScalarT c = 0.0; + for (j = 0; j < k; ++j) + { + ScalarT y = (alpha * V_data[n * j + i] * y_data[j]) - c; + ScalarT t = sum + y; + c = (t - sum) - y; + sum = t; + } + x_data[i] = sum; + } + break; + default: + out::error() << "Unrecognized transpose option " << transpose + << " in gemv. Valid options are 'N' (not transposed) and 'T' (transposed).\n"; + } // switch + x->setDataUpdated(memory::HOST); + } + + /** + * @brief mass (bulk) axpy i.e, y = y - x*alpha where alpha is a vector + * + * @param[in] size number of elements in y + * @param[in] alpha vector size k x 1 + * @param[in] x (multi)vector size size x k + * @param[in,out] y vector size size x 1 (this is where the result is stored) + * + * @pre _k_ > 0, _size_ > 0, _size_ = x->getSize() + */ + template + void VectorHandlerCpu::axpyMulti(IdxT size, + Vector* alpha, + IdxT k, + Vector* x, + Vector* y) + { + ScalarT* alpha_data = alpha->getData(memory::HOST); + ScalarT* y_data = y->getData(memory::HOST); + ScalarT* x_data = x->getData(memory::HOST); + IdxT i, j; + ScalarT sum; + + for (i = 0; i < size; ++i) + { + sum = 0.0; + for (j = 0; j < k; ++j) + { + sum += x_data[j * size + i] * alpha_data[j]; + } + y_data[i] = y_data[i] - sum; + } + y->setDataUpdated(memory::HOST); + } + + /** + * @brief mass (bulk) dot product i.e, V^T x, where V is n x k dense multivector (a dense + * multivector consisting of k vectors size n) and x is k x 2 dense multivector (a multivector + * consisting of two vectors size n each) + * + * @param[in] size Number of elements in a single vector in V + * @param[in] V Multivector; k vectors size n x 1 each + * @param[in] k Number of vectors in V + * @param[in] x Multivector; 2 vectors size n x 1 each + * @param[out] res Multivector; 2 vectors size k x 1 each (result is returned in res) + * + * @pre _size_ > 0, _k_ > 0, size = x->getSize(), _res_ needs to be allocated + */ + template + void VectorHandlerCpu::dot2Multi(IdxT size, + Vector* V, + IdxT k, + Vector* x, + Vector* res) + { + ScalarT* res_data = res->getData(memory::HOST); + const ScalarT* x_data = x->getData(memory::HOST); + const ScalarT* V_data = V->getData(memory::HOST); + + ScalarT c0 = 0.0; + ScalarT cq = 0.0; + + for (IdxT i = 0; i < k; ++i) + { + res_data[i] = 0.0; + res_data[i + k] = 0.0; + + // Make sure we don't accumulate round-off errors + for (IdxT j = 0; j < size; ++j) + { + ScalarT y0 = (V_data[i * size + j] * x_data[j]) - c0; + ScalarT yq = (V_data[i * size + j] * x_data[j + size]) - cq; + ScalarT t0 = res_data[i] + y0; + ScalarT tq = res_data[i + k] + yq; + c0 = (t0 - res_data[i]) - y0; + cq = (tq - res_data[i + k]) - yq; + + res_data[i] = t0; + res_data[i + k] = tq; + } + } + res->setDataUpdated(memory::HOST); + } + + /** + * @brief Scale a vector by a diagonal matrix + * + * @param[in] diag Diagonal vector + * @param[in,out] vec Vector to be scaled + */ + template + void VectorHandlerCpu::scal(Vector* diag, Vector* vec) + { + const ScalarT* diag_data = diag->getData(memory::HOST); + ScalarT* vec_data = vec->getData(memory::HOST); + IdxT n = vec->getSize(); + + for (IdxT i = 0; i < n; ++i) + { + vec_data[i] *= diag_data[i]; + } + vec->setDataUpdated(memory::HOST); + } + + /** + * @brief Scale a vector by a diagonal matrix represented by a contiguous subvector of diag + * + * @param[in] diag Diagonal vector + * @param[in,out] vec Vector to be scaled + * @param[in] diag_offset - the index of diag where the diagonal matrix begins + */ + template + void VectorHandlerCpu::scal(Vector* diag, Vector* vec, IdxT diag_offset) + { + const ScalarT* diag_data = &diag->getData(memory::HOST)[diag_offset]; + ScalarT* vec_data = vec->getData(memory::HOST); + IdxT n = vec->getSize(); + + for (IdxT i = 0; i < n; ++i) + { + vec_data[i] *= diag_data[i]; + } + vec->setDataUpdated(memory::HOST); + } + + /** + * @brief Multiplies vector by an inverse of a diagonal matrix. + * + * @param[in] diag - diagonal matrix stored in a vector object + * @param[in,out] vec - vector to be divided + * + * @pre The two vectors must be the same size + * + * @return 0 if successful, 1 otherwise + */ + template + int VectorHandlerCpu::diagSolve(Vector* diag, Vector* vec) + { + ScalarT* diag_data = diag->getData(memory::HOST); + ScalarT* vec_data = vec->getData(memory::HOST); + IdxT n = vec->getSize(); + + for (IdxT i = 0; i < n; ++i) + { + vec_data[i] /= diag_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 input vector + * @param[in] y - Second input vector + * @param[out] out - Output vector + * + * @pre The three vectors must be the same size + * + * @return 0 if successful, 1 otherwise + */ + template + int VectorHandlerCpu::max(Vector* x, Vector* y, Vector* out) + { + const ScalarT* x_data = x->getData(memory::HOST); + const ScalarT* y_data = y->getData(memory::HOST); + ScalarT* out_data = out->getData(memory::HOST); + IdxT n = y->getSize(); + + for (IdxT i = 0; i < n; ++i) + { + out_data[i] = std::max(x_data[i], y_data[i]); + } + out->setDataUpdated(memory::HOST); + return 0; + } + + /** + * @brief Take the element-wise absolute value of a vector. + * + * @param[in] in - Input vector + * @param[out] out - Output vector + * + * @return 0 if successful, 1 otherwise + */ + template + int VectorHandlerCpu::abs(Vector* in, Vector* out) + { + const ScalarT* in_data = in->getData(memory::HOST); + ScalarT* out_data = out->getData(memory::HOST); + IdxT n = in->getSize(); + + for (IdxT i = 0; i < n; ++i) + { + out_data[i] = std::abs(in_data[i]); + } + out->setDataUpdated(memory::HOST); + return 0; + } + + // Available template instantiations + template class VectorHandlerCpu; + + } // namespace LinearAlgebra +} // namespace GridKit diff --git a/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.hpp b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.hpp new file mode 100644 index 000000000..6fd527dbd --- /dev/null +++ b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.hpp @@ -0,0 +1,74 @@ +#pragma once + +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + template + class Vector; + + /** + * @brief Implements dense (multi)vector kernels that operate on HOST data. + * + * All methods in this class assume the vectors passed in have up-to-date + * data on HOST and operate directly on the raw HOST arrays. + * + * @author Slaven Peles + */ + template + class VectorHandlerCpu + { + public: + VectorHandlerCpu() = default; + ~VectorHandlerCpu() = default; + + // y = alpha * x + y + void axpy(const ScalarT alpha, Vector* x, Vector* y); + + // dot: x \cdot y + ScalarT dot(Vector* x, Vector* y); + + // scal: x = alpha * x + void scal(const ScalarT alpha, Vector* x); + + // vector infinity norm + ScalarT amax(Vector* x); + + // mass axpy: y = y - x*alpha, where x is [n x k] and alpha is [k x 1]; x is stored columnwise + void axpyMulti(IdxT size, Vector* alpha, IdxT k, Vector* x, Vector* y); + + // mass dot: V^T x, where V is [n x k] and x is [k x 2], everything is stored and returned columnwise + void dot2Multi(IdxT size, Vector* V, IdxT k, Vector* x, Vector* res); + + /** 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]`. + * if `transpose = T` (yes), `x = beta*x + alpha*V^T*y`, + * where `x` is `[k x 1]`, `V` is `[n x k]` and `y` is `[n x 1]`. + */ + void gemv(char transpose, + IdxT k, + const ScalarT alpha, + const ScalarT beta, + Vector* V, + Vector* y, + Vector* x); + + // Scale a vector by a diagonal matrix + void scal(Vector* diag, Vector* vec); + void scal(Vector* diag, Vector* vec, IdxT diag_offset); + + // Divide the elements of a vector by the elements of another vector + int diagSolve(Vector* diag, Vector* vec); + + // Compute element-wise max of two vectors + int max(Vector* x, Vector* y, Vector* out); + + // Compute element-wise absolute value of a vector + int abs(Vector* in, Vector* out); + }; + + } // namespace LinearAlgebra +} // namespace GridKit diff --git a/tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt b/tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt index 11105d0f4..889206aea 100644 --- a/tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt +++ b/tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt @@ -5,3 +5,11 @@ target_link_libraries( add_test(NAME VectorTest COMMAND $) install(TARGETS test_vector RUNTIME DESTINATION bin) + +add_executable(test_vector_handler runVectorHandlerTests.cpp) +target_link_libraries( + test_vector_handler GridKit::dense_vector GridKit::testing) + +add_test(NAME VectorHandlerTest COMMAND $) + +install(TARGETS test_vector_handler RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp b/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp new file mode 100644 index 000000000..b8d1a3ea0 --- /dev/null +++ b/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp @@ -0,0 +1,509 @@ +#pragma once +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + using namespace LinearAlgebra; + + /** + * @class Tests for the vector handler. + */ + template + class VectorHandlerTests + { + public: + VectorHandlerTests(VectorHandler& handler, memory::MemorySpace memspace = memory::HOST) + : handler_(handler), + memspace_(memspace) + { + } + + virtual ~VectorHandlerTests() + { + } + + TestOutcome vectorHandlerConstructor() + { + TestStatus status; + status.skipTest(); + + return status.report(__func__); + } + + /** + * @brief Test vector infinity norm. + */ + TestOutcome amax(IdxT N) + { + TestStatus status; + status = true; + + Vector x(N); + + ScalarT* data = new ScalarT[N]; + for (IdxT i = 0; i < N; ++i) + { + data[i] = static_cast(0.1) * static_cast(i); + } + x.allocate(memspace_); + x.copyFromExternal(data, memory::HOST, memspace_); + + ScalarT result = handler_.amax(&x, memspace_); + ScalarT answer = static_cast(N - 1) * static_cast(0.1); + + if (!isEqual(result, answer)) + { + std::cout << "The result " << result << " is incorrect. " + << "Expected answer is " << answer << "\n"; + status *= false; + } + + delete[] data; + return status.report(__func__); + } + + /** + * @brief Test axpy: y = alpha*x + y. + */ + TestOutcome axpy(IdxT N) + { + TestStatus status; + status = true; + + Vector x(N); + Vector y(N); + + x.allocate(memspace_); + y.allocate(memspace_); + + x.setToConst(3.0, memspace_); + y.setToConst(1.0, memspace_); + + ScalarT alpha = static_cast(0.5); + + // the result is a vector with y[i] = 2.5 for all i; + handler_.axpy(alpha, &x, &y, memspace_); + status *= verifyAnswer(y, 2.5); + + return status.report(__func__); + } + + /** + * @brief Test dot product. + */ + TestOutcome dot(IdxT N) + { + TestStatus status; + status = true; + + Vector x(N); + Vector y(N); + + x.allocate(memspace_); + y.allocate(memspace_); + + x.setToConst(0.25, memspace_); + y.setToConst(4.0, memspace_); + + // the answer is N + ScalarT answer = static_cast(N); + ScalarT result = handler_.dot(&x, &y, memspace_); + + if (!isEqual(result, answer)) + { + std::cout << "The result " << result << " is incorrect. " + << "Expected answer is " << answer << "\n"; + status *= false; + } + + return status.report(__func__); + } + + /** + * @brief Test scal: x = alpha*x. + */ + TestOutcome scal(IdxT N) + { + TestStatus status; + status = true; + + Vector x(N); + + x.allocate(memspace_); + + x.setToConst(1.25, memspace_); + + ScalarT alpha = 3.5; + + // the answer is x[i] = 4.375; + ScalarT answer = 4.375; + handler_.scal(alpha, &x, memspace_); + status *= verifyAnswer(x, answer); + + return status.report(__func__); + } + + /** + * @brief Test mass (bulk) axpy: y = y - x*alpha, where x is [N x K] and alpha is [K x 1]. + */ + TestOutcome axpyMulti(IdxT N, IdxT K) + { + TestStatus status; + status = true; + + Vector x(N, K); + Vector y(N); + Vector alpha(K); + + x.allocate(memspace_); + y.allocate(memspace_); + alpha.allocate(memspace_); + + alpha.setToConst(-1.0, memspace_); + y.setToConst(2.0, memspace_); + + for (IdxT ii = 0; ii < K; ++ii) + { + ScalarT c; + if (ii % 2 == 0) + { + c = -1.0; + } + else + { + c = static_cast(0.5); + } + x.setToConst(ii, c, memspace_); + } + + IdxT r = K % 2; + ScalarT res = static_cast((std::floor(static_cast(K) / 2.0) + static_cast(r)) * 1.0 + + std::floor(static_cast(K) / 2.0) * (-0.5)); + + handler_.axpyMulti(N, &alpha, K, &x, &y, memspace_); + status *= verifyAnswer(y, 2.0 - res); + + return status.report(__func__); + } + + /** + * @brief Test mass (bulk) dot product: V^T x. + */ + TestOutcome massDot(IdxT N, IdxT K) + { + TestStatus status; + status = true; + + Vector x(N, K); + Vector y(N, 2); + Vector res(K, 2); + x.allocate(memspace_); + y.allocate(memspace_); + res.allocate(memspace_); + // res is a write-only output; mark it up to date before it is used. + res.setToZero(memspace_); + + x.setToConst(1.0, memspace_); + y.setToConst(-1.0, memspace_); + handler_.dot2Multi(N, &x, K, &y, &res, memspace_); + + status *= verifyAnswer(res, (-1.0) * static_cast(N)); + + return status.report(__func__); + } + + /** + * @brief Test dense matrix-vector product (gemv), both transposed and not. + */ + TestOutcome gemv(IdxT N, IdxT K) + { + TestStatus status; + status = true; + + Vector V(N, K); + Vector yN(K); ///< For the test with NO TRANSPOSE + Vector xN(N); + Vector yT(N); ///< for the test with TRANSPOSE + Vector xT(K); + + V.allocate(memspace_); + yN.allocate(memspace_); + xN.allocate(memspace_); + yT.allocate(memspace_); + xT.allocate(memspace_); + + V.setToConst(1.0, memspace_); + yN.setToConst(-1.0, memspace_); + xN.setToConst(0.5, memspace_); + yT.setToConst(-1.0, memspace_); + xT.setToConst(0.5, memspace_); + + ScalarT alpha = -1.0; + ScalarT beta = 1.0; + handler_.gemv('N', K, alpha, beta, &V, &yN, &xN, memspace_); + status *= verifyAnswer(xN, static_cast(K) + static_cast(0.5)); + handler_.gemv('T', K, alpha, beta, &V, &yT, &xT, memspace_); + status *= verifyAnswer(xT, static_cast(N) + static_cast(0.5)); + + return status.report(__func__); + } + + /** + * @brief Test scaling a vector by a diagonal matrix. + */ + TestOutcome scale(IdxT N) + { + TestStatus status; + status = true; + + Vector diag(N); + Vector vec(N); + + // diag[i] = i + 1, vec[i] = 3.0 + // expected result vec[i] = (i + 1) * 3.0 + diag.allocate(memspace_); + vec.allocate(memspace_); + + vec.setToConst(3.0, memspace_); + + auto diag_data = std::unique_ptr(new ScalarT[N]); + for (IdxT i = 0; i < N; ++i) + { + diag_data[i] = static_cast(i + 1); + } + diag.copyFromExternal(diag_data.get(), memory::HOST, memspace_); + + handler_.scal(&diag, &vec, memspace_); + + if (memspace_ == memory::DEVICE) + { + vec.syncData(memory::HOST); + } + + for (IdxT i = 0; i < N; ++i) + { + if (!isEqual(vec.getData(memory::HOST)[i], static_cast(i + 1) * 3.0)) + { + std::cout << "Solution vector element vec[" << i << "] = " << vec.getData(memory::HOST)[i] + << ", expected: " << static_cast(i + 1) * 3.0 << "\n"; + status *= false; + break; + } + } + + return status.report(__func__); + } + + /** + * @brief Test dividing a vector by a diagonal matrix. + */ + TestOutcome diagSolve(IdxT N) + { + TestStatus status; + status = true; + + Vector diag(N); + Vector vec(N); + + // diag[i] = i + 1, vec[i] = 3.0 + // expected result vec[i] = 3.0 / (i + 1) + diag.allocate(memspace_); + vec.allocate(memspace_); + + vec.setToConst(3.0, memspace_); + + auto diag_data = std::unique_ptr(new ScalarT[N]); + for (IdxT i = 0; i < N; ++i) + { + diag_data[i] = static_cast(i + 1); + } + diag.copyFromExternal(diag_data.get(), memory::HOST, memspace_); + + handler_.diagSolve(&diag, &vec, memspace_); + + if (memspace_ == memory::DEVICE) + { + vec.syncData(memory::HOST); + } + + for (IdxT i = 0; i < N; ++i) + { + if (!isEqual(vec.getData(memory::HOST)[i], static_cast(3.0) / static_cast(i + 1))) + { + std::cout << "Solution vector element vec[" << i << "] = " << vec.getData(memory::HOST)[i] + << ", expected: " << static_cast(3.0) / static_cast(i + 1) << "\n"; + status *= false; + break; + } + } + + return status.report(__func__); + } + + /** + * @brief Test element-wise max of two vectors. + */ + TestOutcome max(IdxT N) + { + TestStatus status; + status = true; + + Vector x(N); + Vector y(N); + Vector z(N); + + x.allocate(memspace_); + y.allocate(memspace_); + z.allocate(memspace_); + // z is a write-only output; mark it up to date before it is used. + z.setToZero(memspace_); + + auto x_data = std::unique_ptr(new ScalarT[N]); + auto y_data = std::unique_ptr(new ScalarT[N]); + for (IdxT i = 0; i < N; ++i) + { + if (i % 3 == 0) + { + x_data[i] = static_cast(i + 1); + y_data[i] = static_cast(i) * 0.5; + } + else + { + x_data[i] = -static_cast(i + 1); + y_data[i] = static_cast(i + 1); + } + } + x.copyFromExternal(x_data.get(), memory::HOST, memspace_); + y.copyFromExternal(y_data.get(), memory::HOST, memspace_); + + handler_.max(&x, &y, &z, memspace_); + handler_.max(&x, &y, &y, memspace_); + + if (memspace_ == memory::DEVICE) + { + y.syncData(memory::HOST); + z.syncData(memory::HOST); + } + + for (IdxT i = 0; i < N; ++i) + { + if (!isEqual(y.getData(memory::HOST)[i], static_cast(i + 1))) + { + std::cout << "Solution vector element y[" << i << "] = " << y.getData(memory::HOST)[i] + << ", expected: " << static_cast(i + 1) << "\n"; + status *= false; + break; + } + + if (!isEqual(z.getData(memory::HOST)[i], static_cast(i + 1))) + { + std::cout << "Solution vector element z[" << i << "] = " << z.getData(memory::HOST)[i] + << ", expected: " << static_cast(i + 1) << "\n"; + status *= false; + break; + } + } + + return status.report(__func__); + } + + /** + * @brief Test element-wise absolute value of a vector. + */ + TestOutcome abs(IdxT N) + { + TestStatus status; + status = true; + + Vector x(N); + Vector y(N); + + x.allocate(memspace_); + y.allocate(memspace_); + // y is a write-only output; mark it up to date before it is used. + y.setToZero(memspace_); + + auto x_data = std::unique_ptr(new ScalarT[N]); + for (IdxT i = 0; i < N; ++i) + { + if (i % 3 == 0) + { + x_data[i] = -static_cast(i); + } + else + { + x_data[i] = static_cast(i); + } + } + x.copyFromExternal(x_data.get(), memory::HOST, memspace_); + + handler_.abs(&x, &y, memspace_); + handler_.abs(&x, &x, memspace_); + + if (memspace_ == memory::DEVICE) + { + x.syncData(memory::HOST); + y.syncData(memory::HOST); + } + + for (IdxT i = 0; i < N; ++i) + { + if (!isEqual(x.getData(memory::HOST)[i], static_cast(i))) + { + std::cout << "Solution vector element x[" << i << "] = " << x.getData(memory::HOST)[i] + << ", expected: " << static_cast(i) << "\n"; + status *= false; + break; + } + + if (!isEqual(y.getData(memory::HOST)[i], static_cast(i))) + { + std::cout << "Solution vector element y[" << i << "] = " << y.getData(memory::HOST)[i] + << ", expected: " << static_cast(i) << "\n"; + status *= false; + break; + } + } + + return status.report(__func__); + } + + private: + VectorHandler& handler_; + memory::MemorySpace memspace_{memory::HOST}; + + // we can verify through norm but that would defeat the purpose of testing vector handler ... + bool verifyAnswer(Vector& x, ScalarT answer) + { + bool success = true; + + if (memspace_ == memory::DEVICE) + { + x.syncData(memory::HOST); + } + + for (IdxT i = 0; i < x.getSize(); ++i) + { + if (!isEqual(x.getData(memory::HOST)[i], answer)) + { + std::cout << std::setprecision(16); + success = false; + std::cout << "Solution vector element x[" << i << "] = " << x.getData(memory::HOST)[i] + << ", expected: " << answer << "\n"; + break; + } + } + return success; + } + }; // class + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/LinearAlgebra/Vector/runVectorHandlerTests.cpp b/tests/UnitTests/LinearAlgebra/Vector/runVectorHandlerTests.cpp new file mode 100644 index 000000000..acb10edf4 --- /dev/null +++ b/tests/UnitTests/LinearAlgebra/Vector/runVectorHandlerTests.cpp @@ -0,0 +1,44 @@ +#include +#include +#include + +#include "VectorHandlerTests.hpp" + +int main(int, char**) +{ + GridKit::Testing::TestingResults result; + + { + std::cout << "Running vector handler tests on CPU:\n"; + + GridKit::LinearAlgebra::VectorHandler handler; + + GridKit::Testing::VectorHandlerTests test(handler); + result += test.vectorHandlerConstructor(); + result += test.dot(50); + result += test.axpy(50); + result += test.scal(50); + result += test.amax(50); + result += test.gemv(500, 10); + result += test.axpyMulti(100, 10); + result += test.massDot(100, 10); + result += test.scale(100); + result += test.diagSolve(100); + result += test.max(100); + result += test.abs(100); + + std::cout << "\n"; + } + +#ifdef GRIDKIT_ENABLE_CUDA + // TODO: Once a CUDA VectorHandler implementation is available, run the + // same tests here with a device-backed VectorHandler and memory::DEVICE. +#endif + +#ifdef GRIDKIT_ENABLE_HIP + // TODO: Once a HIP VectorHandler implementation is available, run the + // same tests here with a device-backed VectorHandler and memory::DEVICE. +#endif + + return result.summary(); +} From 818a3404a67ffa2946f892eff2d76cffa35813b0 Mon Sep 17 00:00:00 2001 From: pelesh Date: Wed, 8 Jul 2026 13:36:54 -0400 Subject: [PATCH 02/10] [skip ci] Update GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com> --- GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp index b46c977c4..f4c009acd 100644 --- a/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp +++ b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp @@ -101,7 +101,7 @@ namespace GridKit } /** - * @brief gemv computes matrix-vector product where both matrix and vectors are dense. + * @brief gemv computes matrix-vector product where the matrix and vectors are dense. * i.e., x = beta*x + alpha*V*y * * @param[in] transpose - transposed = 'T' or not 'N' From 18df88a8b24ca8b6df4f5603fb8d71f5fb28265f Mon Sep 17 00:00:00 2001 From: pelesh Date: Wed, 8 Jul 2026 13:37:23 -0400 Subject: [PATCH 03/10] [skip ci] Update GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com> --- GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp index f4c009acd..a7ce43e29 100644 --- a/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp +++ b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp @@ -306,6 +306,9 @@ namespace GridKit /** * @brief Multiplies vector by an inverse of a diagonal matrix. + * This is equivalent to solving a system with the original matrix, + * with the vector as the right hand side (typically used in this context). + * * * @param[in] diag - diagonal matrix stored in a vector object * @param[in,out] vec - vector to be divided From d71c49ef3854512d9d52262b7db09c443375c1a9 Mon Sep 17 00:00:00 2001 From: pelesh Date: Wed, 8 Jul 2026 14:24:08 -0400 Subject: [PATCH 04/10] [skip ci] Update tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com> --- tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp b/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp index b8d1a3ea0..154ae3373 100644 --- a/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp +++ b/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp @@ -222,7 +222,7 @@ namespace GridKit } /** - * @brief Test dense matrix-vector product (gemv), both transposed and not. + * @brief Test dense matrix-vector product (gemv), transposed and not. */ TestOutcome gemv(IdxT N, IdxT K) { From 65ef47cf8370992cd8f5563127392e764a3abd00 Mon Sep 17 00:00:00 2001 From: pelesh Date: Wed, 8 Jul 2026 14:25:15 -0400 Subject: [PATCH 05/10] [skip ci] Update tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com> --- tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp b/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp index 154ae3373..281c3c722 100644 --- a/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp +++ b/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp @@ -481,7 +481,7 @@ namespace GridKit VectorHandler& handler_; memory::MemorySpace memspace_{memory::HOST}; - // we can verify through norm but that would defeat the purpose of testing vector handler ... + // We could verify through norm but that would defeat the purpose of testing vector handler. bool verifyAnswer(Vector& x, ScalarT answer) { bool success = true; From f3022040c2b697512e255ffb77702726fe9d9380 Mon Sep 17 00:00:00 2001 From: pelesh Date: Wed, 8 Jul 2026 14:29:47 -0400 Subject: [PATCH 06/10] [skip ci] Minor comment clarification. --- tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp b/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp index 281c3c722..2d9cf772a 100644 --- a/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp +++ b/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp @@ -304,7 +304,10 @@ namespace GridKit } /** - * @brief Test dividing a vector by a diagonal matrix. + * @brief Test elementwise division of two vectors. + * + * Mathematically, this is equivalent to multiplying inverse of a + * diagonal matrix to a vector. */ TestOutcome diagSolve(IdxT N) { From 6f4073e7fb78ac227c6b9ef9467b6e4485b6df55 Mon Sep 17 00:00:00 2001 From: pelesh Date: Wed, 8 Jul 2026 17:54:07 -0400 Subject: [PATCH 07/10] MemoryHandler -> MemoryManager --- GridKit/LinearAlgebra/SparseMatrix/CooMatrix.hpp | 2 +- GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp | 2 +- GridKit/LinearAlgebra/Vector/Vector.hpp | 2 +- GridKit/MemoryUtilities/MemoryUtils.hpp | 6 +++--- .../UnitTests/LinearAlgebra/SparseMatrix/SparseCooTests.hpp | 2 +- .../UnitTests/LinearAlgebra/SparseMatrix/SparseCsrTests.hpp | 2 +- tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp | 4 ++-- 7 files changed, 10 insertions(+), 10 deletions(-) diff --git a/GridKit/LinearAlgebra/SparseMatrix/CooMatrix.hpp b/GridKit/LinearAlgebra/SparseMatrix/CooMatrix.hpp index 54313e863..77cc7a541 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/CooMatrix.hpp +++ b/GridKit/LinearAlgebra/SparseMatrix/CooMatrix.hpp @@ -84,7 +84,7 @@ namespace GridKit IdxT* map_to_sorted_ = {nullptr}; ///< map from orginal to sorted IdxT* map_to_dedup_ = {nullptr}; ///< map from sorted to deduplicated - MemoryHandler mem_; ///< Device memory manager object + MemoryManager mem_; ///< Device memory manager object }; } // namespace LinearAlgebra } // namespace GridKit diff --git a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp index d661a07cd..74029a0da 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp +++ b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp @@ -102,7 +102,7 @@ namespace GridKit bool owns_gpu_sparsity_pattern_{false}; ///< for row/col data bool owns_gpu_values_{false}; ///< for nonzero values - MemoryHandler mem_; ///< Device memory manager object + MemoryManager mem_; ///< Device memory manager object }; } // namespace LinearAlgebra } // namespace GridKit diff --git a/GridKit/LinearAlgebra/Vector/Vector.hpp b/GridKit/LinearAlgebra/Vector/Vector.hpp index f337d71e2..05214ffaf 100644 --- a/GridKit/LinearAlgebra/Vector/Vector.hpp +++ b/GridKit/LinearAlgebra/Vector/Vector.hpp @@ -95,7 +95,7 @@ namespace GridKit bool owns_gpu_data_{true}; ///< data owneship flag for DEVICE data bool owns_cpu_data_{true}; ///< data ownership flag for HOST data - MemoryHandler mem_; ///< Device memory manager object + MemoryManager mem_; ///< Device memory manager object }; } // namespace LinearAlgebra } // namespace GridKit diff --git a/GridKit/MemoryUtilities/MemoryUtils.hpp b/GridKit/MemoryUtilities/MemoryUtils.hpp index ca9419932..926fdc2bc 100644 --- a/GridKit/MemoryUtilities/MemoryUtils.hpp +++ b/GridKit/MemoryUtilities/MemoryUtils.hpp @@ -155,10 +155,10 @@ namespace GridKit // Check if GPU support is enabled in GridKit and set appropriate device memory manager. #if defined GridKit_ENABLE_CUDA #include -using MemoryHandler = GridKit::MemoryUtils; +using MemoryManager = GridKit::MemoryUtils; #elif defined GridKit_ENABLE_HIP #include -using MemoryHandler = GridKit::MemoryUtils; +using MemoryManager = GridKit::MemoryUtils; #else #error Unrecognized device, probably bug in CMake configuration #endif @@ -167,6 +167,6 @@ using MemoryHandler = GridKit::MemoryUtils; // If no GPU support is present, set device memory manager to a dummy object. #include -using MemoryHandler = GridKit::MemoryUtils; +using MemoryManager = GridKit::MemoryUtils; #endif diff --git a/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCooTests.hpp b/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCooTests.hpp index baed39c19..9ee684253 100644 --- a/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCooTests.hpp +++ b/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCooTests.hpp @@ -172,7 +172,7 @@ namespace GridKit private: memory::MemorySpace memspace_; - MemoryHandler mem_; + MemoryManager mem_; }; // class SparseCooTests } // namespace Testing } // namespace GridKit diff --git a/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCsrTests.hpp b/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCsrTests.hpp index d3231c542..cc9cac5b6 100644 --- a/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCsrTests.hpp +++ b/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCsrTests.hpp @@ -392,7 +392,7 @@ namespace GridKit private: memory::MemorySpace memspace_; - MemoryHandler mem_; + MemoryManager mem_; }; // class SparseTests } // namespace Testing } // namespace GridKit diff --git a/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp b/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp index afc4a7f30..c6b48a189 100644 --- a/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp +++ b/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp @@ -221,7 +221,7 @@ namespace GridKit * @brief Test copying data from vector to an array. * * This creates a vector, copies data to it, and then copies the data to an array - * in the current memory space. Finally, it uses the MemoryHandler to copy the data + * in the current memory space. Finally, it uses the MemoryManager to copy the data * to HOST for verification. * * @param[in] N Number of elements in the vector. @@ -366,7 +366,7 @@ namespace GridKit private: memory::MemorySpace memspace_{memory::HOST}; - MemoryHandler mh_; + MemoryManager mh_; /// Check if vector elements are set to the same number bool verifyAnswer(Vector& x, ScalarT answer) From e223b93265b409654eb1f71fde0bad4bfff91aed Mon Sep 17 00:00:00 2001 From: pelesh Date: Wed, 8 Jul 2026 21:54:38 +0000 Subject: [PATCH 08/10] Apply pre-commit fixes --- GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp | 2 +- tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp index a7ce43e29..979ea1205 100644 --- a/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp +++ b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp @@ -306,7 +306,7 @@ namespace GridKit /** * @brief Multiplies vector by an inverse of a diagonal matrix. - * This is equivalent to solving a system with the original matrix, + * This is equivalent to solving a system with the original matrix, * with the vector as the right hand side (typically used in this context). * * diff --git a/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp b/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp index 2d9cf772a..748a065f0 100644 --- a/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp +++ b/tests/UnitTests/LinearAlgebra/Vector/VectorHandlerTests.hpp @@ -305,7 +305,7 @@ namespace GridKit /** * @brief Test elementwise division of two vectors. - * + * * Mathematically, this is equivalent to multiplying inverse of a * diagonal matrix to a vector. */ From c3929342ff4daa26d970cb734e648c241884c786 Mon Sep 17 00:00:00 2001 From: pelesh Date: Thu, 9 Jul 2026 07:19:31 -0400 Subject: [PATCH 09/10] Protect MemoryManager with a namespace. --- GridKit/MemoryUtilities/MemoryUtils.hpp | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/GridKit/MemoryUtilities/MemoryUtils.hpp b/GridKit/MemoryUtilities/MemoryUtils.hpp index 926fdc2bc..4e6c7bf30 100644 --- a/GridKit/MemoryUtilities/MemoryUtils.hpp +++ b/GridKit/MemoryUtilities/MemoryUtils.hpp @@ -155,10 +155,8 @@ namespace GridKit // Check if GPU support is enabled in GridKit and set appropriate device memory manager. #if defined GridKit_ENABLE_CUDA #include -using MemoryManager = GridKit::MemoryUtils; #elif defined GridKit_ENABLE_HIP #include -using MemoryManager = GridKit::MemoryUtils; #else #error Unrecognized device, probably bug in CMake configuration #endif @@ -167,6 +165,18 @@ using MemoryManager = GridKit::MemoryUtils; // If no GPU support is present, set device memory manager to a dummy object. #include -using MemoryManager = GridKit::MemoryUtils; #endif + +namespace GridKit +{ +#ifdef GridKit_ENABLE_GPU +#if defined GridKit_ENABLE_CUDA + using MemoryManager = MemoryUtils; +#elif defined GridKit_ENABLE_HIP + using MemoryManager = MemoryUtils; +#endif +#else + using MemoryManager = MemoryUtils; +#endif +} // namespace GridKit From 23722226a8c4f9bdc04d2f649ed0263065343914 Mon Sep 17 00:00:00 2001 From: pelesh Date: Thu, 9 Jul 2026 07:29:10 -0400 Subject: [PATCH 10/10] Allow for int indices in VectorHandler --- GridKit/LinearAlgebra/Vector/Vector.cpp | 6 +++--- GridKit/LinearAlgebra/Vector/VectorHandler.cpp | 1 + GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp | 1 + 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/GridKit/LinearAlgebra/Vector/Vector.cpp b/GridKit/LinearAlgebra/Vector/Vector.cpp index bbbf82c98..819f2de48 100644 --- a/GridKit/LinearAlgebra/Vector/Vector.cpp +++ b/GridKit/LinearAlgebra/Vector/Vector.cpp @@ -39,8 +39,8 @@ namespace GridKit : n_capacity_(n), k_(k), n_size_(n), - gpu_updated_(new bool[k]), - cpu_updated_(new bool[k]) + gpu_updated_(new bool[static_cast(k)]), + cpu_updated_(new bool[static_cast(k)]) { setHostUpdated(false); setDeviceUpdated(false); @@ -1091,7 +1091,7 @@ namespace GridKit // template class Vector; template class Vector; - // template class Vector; + template class Vector; } // namespace LinearAlgebra } // namespace GridKit diff --git a/GridKit/LinearAlgebra/Vector/VectorHandler.cpp b/GridKit/LinearAlgebra/Vector/VectorHandler.cpp index 817ee4df8..edd941d99 100644 --- a/GridKit/LinearAlgebra/Vector/VectorHandler.cpp +++ b/GridKit/LinearAlgebra/Vector/VectorHandler.cpp @@ -348,6 +348,7 @@ namespace GridKit // Available template instantiations template class VectorHandler; + template class VectorHandler; } // namespace LinearAlgebra } // namespace GridKit diff --git a/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp index 979ea1205..331e61685 100644 --- a/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp +++ b/GridKit/LinearAlgebra/Vector/VectorHandlerCpu.cpp @@ -386,6 +386,7 @@ namespace GridKit // Available template instantiations template class VectorHandlerCpu; + template class VectorHandlerCpu; } // namespace LinearAlgebra } // namespace GridKit