diff --git a/CHANGELOG.md b/CHANGELOG.md index 0c31a6e41..066786e09 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -61,6 +61,7 @@ - Added `BusToSignalAdapter` component for communicating bus voltages and injection currents. - Added cmake-format hooks, including in pre-commit. - Added off-nominal tap ratio and phase shift support to the PhasorDynamics `Branch` model. +- Added portable Vector class to GridKit ## v0.1 diff --git a/GridKit/LinearAlgebra/CMakeLists.txt b/GridKit/LinearAlgebra/CMakeLists.txt index 4586c24e9..98d64e05a 100644 --- a/GridKit/LinearAlgebra/CMakeLists.txt +++ b/GridKit/LinearAlgebra/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(SparseMatrix) add_subdirectory(DenseMatrix) +add_subdirectory(Vector) install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/MemoryUtils.hpp DESTINATION include/GridKit/LinearAlgebra) diff --git a/GridKit/LinearAlgebra/Vector/CMakeLists.txt b/GridKit/LinearAlgebra/Vector/CMakeLists.txt new file mode 100644 index 000000000..d4384dc15 --- /dev/null +++ b/GridKit/LinearAlgebra/Vector/CMakeLists.txt @@ -0,0 +1,7 @@ +gridkit_add_library(dense_vector + SOURCES + Vector.cpp + HEADERS + Vector.hpp + LINK_LIBRARIES + GridKit::utilities_logger) diff --git a/GridKit/LinearAlgebra/Vector/Vector.cpp b/GridKit/LinearAlgebra/Vector/Vector.cpp new file mode 100644 index 000000000..abb77228e --- /dev/null +++ b/GridKit/LinearAlgebra/Vector/Vector.cpp @@ -0,0 +1,1007 @@ +#include +#include + +#include +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + + using out = GridKit::Utilities::Logger; + + /** + * @brief Single vector constructor. + * + * @param[in] n - Number of elements in the vector + */ + template + Vector::Vector(IdxT n) + : n_capacity_(n), + k_(1), + n_size_(n), + gpu_updated_(new bool[1]), + cpu_updated_(new bool[1]) + { + gpu_updated_[0] = false; + cpu_updated_[0] = false; + } + + /** + * @brief Multivector constructor. + * + * @param[in] n - Number of elements in the vector + * @param[in] k - Number of vectors in multivector + */ + template + Vector::Vector(IdxT n, IdxT k) + : n_capacity_(n), + k_(k), + n_size_(n), + gpu_updated_(new bool[k]), + cpu_updated_(new bool[k]) + { + setHostUpdated(false); + setDeviceUpdated(false); + } + + /** + * @brief destructor. + * + */ + template + Vector::~Vector() + { + if (owns_cpu_data_ && h_data_) + mem_.deleteOnHost(h_data_); + if (owns_gpu_data_ && d_data_) + mem_.deleteOnDevice(d_data_); + delete[] gpu_updated_; + delete[] cpu_updated_; + } + + /** + * @brief Get capacity of a single vector. + * + * Vector memory is allocated to `n_capacity_*k_`. This is the maximum + * number of elements that the (multi)vector can hold. + * + * @return `n_capacity_` the maximum number of elements in the vector. + */ + template + IdxT Vector::getCapacity() const + { + return n_capacity_; + } + + /** + * @brief Get the number of elements in a single vector. + * + * For vectors with changing sizes, set the vector capacity to + * the maximum expected size. + * + * @return `n_size_` number of elements currently in the vector. + */ + template + IdxT Vector::getSize() const + { + return n_size_; + } + + /** + * @brief Get the number of vectors in multivector. + * + * @return _k_, number of vectors in the multivector, + * or 1 if the vector is not a multivector. + */ + template + IdxT Vector::getNumVectors() const + { + return k_; + } + + /** + * @brief Set the vector data pointer (HOST or DEVICE) to an external data. + * + * @param[in] data - Pointer to data + * @param[in] memspace - Memory space (HOST or DEVICE) + * + * @pre Vector data pointers must be null. If the vector data already + * exists this function returns error message. + * + * @warning This function DOES NOT ALLOCATE any data, it only assigns the + * pointer. + * + * @warning This is an expert level method. Use only if you know what + * you are doing. + */ + template + int Vector::setData(ScalarT* data, memory::MemorySpace memspace) + { + using namespace memory; + if (d_data_ || h_data_) + { + out::error() << "Vector::setData - data already exists, ignoring call\n"; + return 1; + } + + switch (memspace) + { + case HOST: + h_data_ = data; + setHostUpdated(true); + setDeviceUpdated(false); + owns_cpu_data_ = false; + break; + case DEVICE: + d_data_ = data; + setHostUpdated(false); + setDeviceUpdated(true); + owns_gpu_data_ = false; + break; + } + return 0; + } + + /** + * @brief Set the flag to indicate that the data (HOST or DEVICE) has been + * updated. + * + * Use this function if you update vector elements by accessing the raw data + * pointer. + * + * @param[in] memspace - Memory space (HOST or DEVICE) + * + * @warning This is an expert level method. Use only if you know what + * you are doing. + */ + template + int Vector::setDataUpdated(memory::MemorySpace memspace) + { + assert(cpu_updated_ && gpu_updated_ && "Update flags not allocated"); + + using namespace memory; + switch (memspace) + { + case HOST: + setHostUpdated(true); + setDeviceUpdated(false); + break; + case DEVICE: + setHostUpdated(false); + setDeviceUpdated(true); + break; + } + return 0; + } + + /** + * @brief Set the flag to indicate that the data (HOST or DEVICE) for + * vector `j` in the multivector has been updated. + * + * Use this function if you update vector elements by accessing the raw data + * pointer. + * + * @param[in] memspace - Memory space (HOST or DEVICE) + * + * @warning This is an expert level method. Use only if you know what + * you are doing. + */ + template + int Vector::setDataUpdated(IdxT j, memory::MemorySpace memspace) + { + assert(cpu_updated_ && gpu_updated_ && "Update flags not allocated"); + + using namespace memory; + switch (memspace) + { + case HOST: + cpu_updated_[j] = true; + gpu_updated_[j] = false; + break; + case DEVICE: + gpu_updated_[j] = true; + cpu_updated_[j] = false; + break; + } + return 0; + } + + /** + * @brief Copy data from another vector. + * + * @param[in] source - Vector whose data will be copied + * @param[in] memspaceSrc - Memory space of the data source (HOST or DEVICE) + * @param[in] memspaceDst - Memory space to copy data to (HOST or DEVICE) + * + * @pre Size of _source_ is greater than or equal to the current vector size. + */ + template + int Vector::copyFromExternal(const Vector& source, memory::MemorySpace memspaceSrc, memory::MemorySpace memspaceDst) + { + const ScalarT* source_data = source.getData(memspaceSrc); + return copyFromExternal(source_data, memspaceSrc, memspaceDst); + } + + /** + * @brief Copy vector data from an input array. + * + * Destination memory must be pre-allocated via allocate() before calling + * this function. + * + * @param[in] source - Array to copy from + * @param[in] memspaceSrc - Memory space of the source array (HOST or DEVICE) + * @param[in] memspaceDst - Memory space to copy data to (HOST or DEVICE) + * + * @return 0 if successful, 1 otherwise. + */ + template + int Vector::copyFromExternal(const ScalarT* source, + memory::MemorySpace memspaceSrc, + memory::MemorySpace memspaceDst) + { + switch (memspaceDst) + { + case memory::HOST: + if (h_data_ == nullptr) + { + out::error() << "Vector::copyFromExternal - host destination not allocated\n"; + return 1; + } + break; + case memory::DEVICE: + if (d_data_ == nullptr) + { + out::error() << "Vector::copyFromExternal - device destination not allocated\n"; + return 1; + } + break; + } + + switch (memspaceSrc) + { + case memory::HOST: + switch (memspaceDst) + { + case memory::HOST: + mem_.copyArrayHostToHost(h_data_, source, n_size_ * k_); + setHostUpdated(true); + setDeviceUpdated(false); + break; + case memory::DEVICE: + mem_.copyArrayHostToDevice(d_data_, source, n_size_ * k_); + setHostUpdated(false); + setDeviceUpdated(true); + break; + default: + return 1; + } + break; + case memory::DEVICE: + switch (memspaceDst) + { + case memory::HOST: + mem_.copyArrayDeviceToHost(h_data_, source, n_size_ * k_); + setHostUpdated(true); + setDeviceUpdated(false); + break; + case memory::DEVICE: + mem_.copyArrayDeviceToDevice(d_data_, source, n_size_ * k_); + setHostUpdated(false); + setDeviceUpdated(true); + break; + default: + return 1; + } + break; + default: + return 1; + } + return 0; + } + + /** + * @brief get a pointer to HOST or DEVICE vector data. + * + * @param[in] memspace - Memory space of the pointer (HOST or DEVICE) + * + * @return pointer to the vector data (HOST or DEVICE). In case of multivectors, + * vectors are stored column-wise. + * + * @note This function gives you access to the pointer, not to a copy. + * If you change the values using the pointer, the vector values will + * change too. Make sure to use setDataUpdated function to set the update + * flags correctly after changing the values. + */ + template + ScalarT* Vector::getData(memory::MemorySpace memspace) + { + using memory::DEVICE; + using memory::HOST; + + switch (memspace) + { + case HOST: + if (cpu_updated_[0] == false) + { + out::error() << "Vector::getData - host data is stale. Perhaps you need to call syncData?\n"; + return nullptr; + } + return h_data_; + case DEVICE: + if (gpu_updated_[0] == false) + { + out::error() << "Vector::getData - device data is stale. Perhaps you need to call syncData?\n"; + return nullptr; + } + return d_data_; + default: + return nullptr; + } + } + + /** + * @brief get a pointer to HOST or DEVICE vector data. + * + * @param[in] memspace - Memory space of the pointer (HOST or DEVICE) + * + * @return pointer to the vector data (HOST or DEVICE). In case of multivectors, + * vectors are stored column-wise. + */ + template + const ScalarT* Vector::getData(memory::MemorySpace memspace) const + { + using memory::DEVICE; + using memory::HOST; + + switch (memspace) + { + case HOST: + if (cpu_updated_[0] == false) + { + out::error() << "Vector::getData - host data is stale. Perhaps you need to call syncData?\n"; + return nullptr; + } + return h_data_; + case DEVICE: + if (gpu_updated_[0] == false) + { + out::error() << "Vector::getData - device data is stale. Perhaps you need to call syncData?\n"; + return nullptr; + } + return d_data_; + default: + return nullptr; + } + } + + /** + * @brief Get a pointer to HOST or DEVICE data of a vector in a multivector. + * + * @param[in] j - Index of a vector in multivector + * @param[in] memspace - Memory space of the pointer (HOST or DEVICE) + * + * @return Pointer to the _j_th vector data (HOST or DEVICE). + * + * @pre `j` < `k_`, i.e., `j` is smaller than the number of vectors. + * + * @note This function gives you access to the pointer, not to a copy. + * If you change the values using the pointer, the vector values will + * change too. Call setDataUpdated() to update the staleness flags. + */ + template + ScalarT* Vector::getData(IdxT j, memory::MemorySpace memspace) + { + using memory::DEVICE; + using memory::HOST; + + if (k_ <= j) + { + out::error() << "Vector::getData - vector index " << j << " out of range, multivector has only " << k_ << " vectors\n"; + return nullptr; + } + + switch (memspace) + { + case HOST: + if (cpu_updated_[j] == false) + { + out::error() << "Vector::getData - host data for vector " << j << " is stale. Perhaps you need to call syncData?\n"; + return nullptr; + } + return &h_data_[j * n_size_]; + case DEVICE: + if (gpu_updated_[j] == false) + { + out::error() << "Vector::getData - device data for vector " << j << " is stale. Perhaps you need to call syncData?\n"; + return nullptr; + } + return &d_data_[j * n_size_]; + default: + return nullptr; + } + } + + /** + * @brief Get a const pointer to HOST or DEVICE data of a vector in a multivector. + * + * @param[in] j - Index of a vector in multivector + * @param[in] memspace - Memory space of the pointer (HOST or DEVICE) + * + * @return Const pointer to the _j_th vector data (HOST or DEVICE). + * + * @pre `j` < `k_`, i.e., `j` is smaller than the number of vectors. + */ + template + const ScalarT* Vector::getData(IdxT j, memory::MemorySpace memspace) const + { + using memory::DEVICE; + using memory::HOST; + + if (k_ <= j) + { + out::error() << "Vector::getData - vector index " << j << " out of range, multivector has only " << k_ << " vectors\n"; + return nullptr; + } + + switch (memspace) + { + case HOST: + if (cpu_updated_[j] == false) + { + out::error() << "Vector::getData - host data for vector " << j << " is stale. Perhaps you need to call syncData?\n"; + return nullptr; + } + return &h_data_[j * n_size_]; + case DEVICE: + if (gpu_updated_[j] == false) + { + out::error() << "Vector::getData - device data for vector " << j << " is stale. Perhaps you need to call syncData?\n"; + return nullptr; + } + return &d_data_[j * n_size_]; + default: + return nullptr; + } + } + + /** + * @brief Sync out of date memory space with the updated one. + * + * syncData is the only function that can set data on both HOST and DEVICE + * to the same values. + * + * @param[in] memspaceDst - Memory space to sync + * + * @return 0 if successful, 1 otherwise. + * + * @warning This function can be called only when all vectors in a + * multivector have the same update status. Otherwise, you need to sync + * vectors in a multivector individually. + * + */ + template + int Vector::syncData(memory::MemorySpace memspaceDst) + { + using namespace memory; + + bool all_cpu_updated = cpu_updated_[0]; + bool all_gpu_updated = gpu_updated_[0]; + + // Verify that all vectors in multivector have the same update status. + for (IdxT i = 1; i < k_; ++i) + { + if (gpu_updated_[i] != all_gpu_updated) + { + out::error() << "Vector::syncData - inconsistent update state across device columns.\n" + << "Use syncData(j, memspace) for individual vectors\n"; + return 1; + } + if (cpu_updated_[i] != all_cpu_updated) + { + out::error() << "Vector::syncData - inconsistent update state across host columns.\n" + << "Use syncData(j, memspace) for individual vectors\n"; + return 1; + } + } + + switch (memspaceDst) + { + case DEVICE: // cpu -> gpu + if (gpu_updated_[0]) + { + out::error() << "Vector::syncData - device already up to date\n"; + return 1; + } + if (!cpu_updated_[0]) + { + out::error() << "Vector::syncData - host data is stale, cannot sync to device\n"; + return 1; + } + if (d_data_ == nullptr) + { + out::error() << "Vector::syncData - device data not allocated\n"; + return 1; + } + mem_.copyArrayHostToDevice(d_data_, h_data_, n_size_ * k_); + setDeviceUpdated(true); + break; + case HOST: // gpu -> cpu + if (cpu_updated_[0]) + { + out::error() << "Vector::syncData - host already up to date\n"; + return 1; + } + if (!gpu_updated_[0]) + { + out::error() << "Vector::syncData - device data is stale, cannot sync to host\n"; + return 1; + } + if (h_data_ == nullptr) + { + out::error() << "Vector::syncData - host data not allocated\n"; + return 1; + } + mem_.copyArrayDeviceToHost(h_data_, d_data_, n_size_ * k_); + setHostUpdated(true); + break; + default: + return 1; + } + return 0; + } + + /** + * @brief Sync out of date memory space with the updated one. + * + * syncData is the only function that can set data on both HOST and DEVICE + * to the same values. + * + * @param[in] memspaceDst - Memory space to sync + * + * @return 0 if successful, 1 otherwise. + * + * @warning This function can be called only when all vectors in a + * multivector have the same update status. Otherwise, you need to sync + * vectors in a multivector individually. + * + */ + template + int Vector::syncData(IdxT j, memory::MemorySpace memspaceDst) + { + using namespace memory; + + switch (memspaceDst) + { + case DEVICE: // cpu->gpu + if (gpu_updated_[j]) + { + out::error() << "Vector::syncData - device already up to date\n"; + return 1; + } + if (!cpu_updated_[j]) + { + out::error() << "Vector::syncData - host data is stale, cannot sync to device\n"; + return 1; + } + if (d_data_ == nullptr) + { + out::error() << "Vector::syncData - device data not allocated\n"; + return 1; + } + mem_.copyArrayHostToDevice(&d_data_[j * n_size_], &h_data_[j * n_size_], n_size_); + gpu_updated_[j] = true; + break; + case HOST: // gpu -> cpu + if (cpu_updated_[j]) + { + out::error() << "Vector::syncData - host already up to date\n"; + return 1; + } + if (!gpu_updated_[j]) + { + out::error() << "Vector::syncData - device data is stale, cannot sync to host\n"; + return 1; + } + if (h_data_ == nullptr) + { + out::error() << "Vector::syncData - host data not allocated\n"; + return 1; + } + mem_.copyArrayDeviceToHost(&h_data_[j * n_size_], &d_data_[j * n_size_], n_size_); + cpu_updated_[j] = true; + break; + default: + return 1; + } + return 0; + } + + /** + * @brief Allocate vector data for HOST or DEVICE + * + * @param[in] memspace - Memory space of the data to be allocated + * + */ + template + int Vector::allocate(memory::MemorySpace memspace) + { + using namespace memory; + switch (memspace) + { + case HOST: + if (!owns_cpu_data_) + { + out::error() << "Vector::allocate - cannot reallocate host data," + << " vector does not own it\n"; + return 1; + } + mem_.deleteOnHost(h_data_); + mem_.allocateArrayOnHost(&h_data_, n_capacity_ * k_); + owns_cpu_data_ = true; + setHostUpdated(false); + break; + case DEVICE: + if (!owns_gpu_data_) + { + out::error() << "Vector::allocate - cannot reallocate device data," + << " vector does not own it\n"; + return 1; + } + mem_.deleteOnDevice(d_data_); + mem_.allocateArrayOnDevice(&d_data_, n_capacity_ * k_); + owns_gpu_data_ = true; + setDeviceUpdated(false); + break; + } + return 0; + } + + /** + * @brief Set vector data to zero. + * + * In case of multivectors, the entire multivector is set to zero. + * + * @param[in] memspace - Memory space of the data to be zeroed (HOST or DEVICE) + * + */ + template + int Vector::setToZero(memory::MemorySpace memspace) + { + using namespace memory; + switch (memspace) + { + case HOST: + if (h_data_ == nullptr) + { + out::error() << "Vector::setToZero - host data not allocated\n"; + return 1; + } + mem_.setZeroArrayOnHost(h_data_, n_capacity_ * k_); + setHostUpdated(true); + setDeviceUpdated(false); + break; + case DEVICE: + if (d_data_ == nullptr) + { + out::error() << "Vector::setToZero - device data not allocated\n"; + return 1; + } + mem_.setZeroArrayOnDevice(d_data_, n_capacity_ * k_); + setHostUpdated(false); + setDeviceUpdated(true); + break; + } + return 0; + } + + /** + * @brief set the data of a single vector in a multivector to zero. + * + * @param[in] j - Index of a vector in the multivector + * @param[in] memspace - Memory space of the data to be zeroed (HOST or DEVICE) + * + * @pre `j` < `k_`, i.e., `j` is smaller than the number of vectors. + */ + template + int Vector::setToZero(IdxT j, memory::MemorySpace memspace) + { + using namespace memory; + switch (memspace) + { + case HOST: + if (h_data_ == nullptr) + { + out::error() << "Vector::setToZero - host data not allocated\n"; + return 1; + } + mem_.setZeroArrayOnHost(&h_data_[j * n_size_], n_size_); + cpu_updated_[j] = true; + gpu_updated_[j] = false; + break; + case DEVICE: + if (d_data_ == nullptr) + { + out::error() << "Vector::setToZero - device data not allocated\n"; + return 1; + } + // TODO: We should not need to access raw data in this class + mem_.setZeroArrayOnDevice(&d_data_[j * n_size_], n_size_); + cpu_updated_[j] = false; + gpu_updated_[j] = true; + break; + } + return 0; + } + + /** + * @brief set vector data to a given constant. + * + * In case of multivectors, entire multivector is set to the constant. + * + * @param[in] C - Constant value to set + * @param[in] memspace - Memory space of the data to be set (HOST or DEVICE) + * + */ + template + int Vector::setToConst(ScalarT C, memory::MemorySpace memspace) + { + using namespace memory; + switch (memspace) + { + case HOST: + if (h_data_ == nullptr) + { + out::error() << "Vector::setToConst - host data not allocated\n"; + return 1; + } + mem_.setArrayToConstOnHost(h_data_, C, n_size_ * k_); + setHostUpdated(true); + setDeviceUpdated(false); + break; + case DEVICE: + if (d_data_ == nullptr) + { + out::error() << "Vector::setToConst - device data not allocated\n"; + return 1; + } + mem_.setArrayToConstOnDevice(d_data_, C, n_size_ * k_); + setHostUpdated(false); + setDeviceUpdated(true); + break; + } + return 0; + } + + /** + * @brief set the data of a single vector in a multivector to a given constant. + * + * @param[in] j - Index of a vector in the multivector + * @param[in] C - Constant value to set + * @param[in] memspace - Memory space of the data to be set (HOST or DEVICE) + * + * @pre `j` < `k_`, i.e., `j` is smaller than the number of vectors. + */ + template + int Vector::setToConst(IdxT j, ScalarT C, memory::MemorySpace memspace) + { + using namespace memory; + switch (memspace) + { + case HOST: + if (h_data_ == nullptr) + { + out::error() << "Vector::setToConst - host data not allocated\n"; + return 1; + } + mem_.setArrayToConstOnHost(&h_data_[n_size_ * j], C, n_size_); + cpu_updated_[j] = true; + gpu_updated_[j] = false; + break; + case DEVICE: + if (d_data_ == nullptr) + { + out::error() << "Vector::setToConst - device data not allocated\n"; + return 1; + } + mem_.setArrayToConstOnDevice(&d_data_[n_size_ * j], C, n_size_); + cpu_updated_[j] = false; + gpu_updated_[j] = true; + break; + } + return 0; + } + + /** + * @brief Resize vector to `new_n_size`. + * + * Use for vectors and multivectors that change size throughout computation. + * + * @note Vector needs to have capacity set to maximum expected size. + * @warning This method is not to be used in vectors who do not own their + * data. + * + * @param[in] new_n_size - New vector length + * + * @return 0 if successful, 1 otherwise. + * + * @pre `new_n_size` <= `n_capacity_` + */ + template + int Vector::resize(IdxT new_n_size) + { + assert(owns_cpu_data_ && owns_gpu_data_ + && "Cannot resize if vector is not owning the data."); + + if (new_n_size > n_capacity_) + { + out::error() << "Vector::resize - requested size " << new_n_size + << " exceeds capacity " << n_capacity_ << "\n"; + return 1; + } + else + { + n_size_ = new_n_size; + return 0; + } + } + + /** + * @brief Copy HOST or DEVICE data of a single vector in a multivector to _dest_. + * + * Supports cross-space copies, e.g. vector _i_ from HOST to DEVICE. + * + * @param[out] dest - Destination array + * @param[in] i - Index of a vector in the multivector + * @param[in] memspaceSrc - Memory space of the source data (HOST or DEVICE) + * @param[in] memspaceDst - Memory space of the destination (HOST or DEVICE) + * + * @return 0 if successful, 1 otherwise. + * + * @pre `i` < `k_`, i.e., `i` is smaller than the number of vectors. + * @pre _dest_ is allocated with at least _n_ elements in _memspaceDst_. + * @post All elements of vector _i_ are copied to _dest_. + */ + template + int Vector::copyToExternal(ScalarT* dest, IdxT i, memory::MemorySpace memspaceSrc, memory::MemorySpace memspaceDst) + { + using namespace memory; + if (i >= k_) + { + out::error() << "Vector::copyToExternal - vector index " << i + << " out of range, multivector has only " << k_ << " vectors\n"; + return 1; + } + if (dest == nullptr) + { + out::error() << "Vector::copyToExternal - destination pointer for vector " << i << " is null\n"; + return 1; + } + ScalarT* data = getData(i, memspaceSrc); + if (data == nullptr) + { + out::error() << "Vector::copyToExternal - source data for vector " << i << " is null or stale\n"; + return 1; + } + switch (memspaceSrc) + { + case HOST: + switch (memspaceDst) + { + case HOST: + mem_.copyArrayHostToHost(dest, data, n_size_); + break; + case DEVICE: + mem_.copyArrayHostToDevice(dest, data, n_size_); + break; + } + break; + case DEVICE: + switch (memspaceDst) + { + case HOST: + mem_.copyArrayDeviceToHost(dest, data, n_size_); + break; + case DEVICE: + mem_.copyArrayDeviceToDevice(dest, data, n_size_); + break; + } + break; + } + return 0; + } + + /** + * @brief copy HOST or DEVICE data of multivector to _dest_. + * + * This function allows to copy data between different memory spaces in one call. + * For example, you can copy data of multivector from HOST to DEVICE, or from DEVICE to HOST. + * + * @param[out] dest - Pointer to the memory to which data is copied + * @param[in] memspaceSrc - Memory space (HOST or DEVICE) of the data to be copied + * @param[in] memspaceDst - Memory space (HOST or DEVICE) to which data is copied + * + * @return 0 if successful, 1 otherwise. + * + * @pre _dest_ is allocated, and the size of _dest_ is at least _n_ * _k_ (total length of all vectors in the multivector). + * @pre _dest_ is allocated in memspaceOutDst memory space. + * @post All elements of all vectors in multivector are copied to the array _dest_. + */ + template + int Vector::copyToExternal(ScalarT* dest, memory::MemorySpace memspaceSrc, memory::MemorySpace memspaceDst) + { + using namespace memory; + ScalarT* data = this->getData(memspaceSrc); + // Check that the source data is not null and up to date + if (data == nullptr) + { + out::error() << "Vector::copyToExternal - source data is null or stale\n"; + return 1; + } + // Check that the destination memory space is allocated + if (dest == nullptr) + { + out::error() << "Vector::copyToExternal - destination pointer is null\n"; + return 1; + } + switch (memspaceSrc) + { + case HOST: + if (!cpu_updated_[0]) + { + out::error() << "Vector::copyToExternal - source data is stale\n"; + return 1; + } + switch (memspaceDst) + { + case HOST: + mem_.copyArrayHostToHost(dest, data, n_size_ * k_); + break; + case DEVICE: + mem_.copyArrayHostToDevice(dest, data, n_size_ * k_); + break; + } + break; + case DEVICE: + if (!gpu_updated_[0]) + { + out::error() << "Vector::copyToExternal - source data is stale\n"; + return 1; + } + switch (memspaceDst) + { + case HOST: + mem_.copyArrayDeviceToHost(dest, data, n_size_ * k_); + break; + case DEVICE: + mem_.copyArrayDeviceToDevice(dest, data, n_size_ * k_); + break; + } + break; + } + return 0; + } + + // + // Private methods + // + + template + void Vector::setHostUpdated(bool is_updated) + { + std::fill(cpu_updated_, cpu_updated_ + k_, is_updated); + } + + template + void Vector::setDeviceUpdated(bool is_updated) + { + std::fill(gpu_updated_, gpu_updated_ + k_, is_updated); + } + + // template class Vector; + template class Vector; + // template class Vector; + + } // namespace LinearAlgebra +} // namespace GridKit diff --git a/GridKit/LinearAlgebra/Vector/Vector.hpp b/GridKit/LinearAlgebra/Vector/Vector.hpp new file mode 100644 index 000000000..2690fb0d9 --- /dev/null +++ b/GridKit/LinearAlgebra/Vector/Vector.hpp @@ -0,0 +1,97 @@ +#pragma once +#include + +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + /** + * @brief This class implements vectors (dense arrays) and multivectors and + * some basic utilities (get size, allocate, set data, get data, etc). + * + * + * What you need to know: + * - Multivectors are stored in one array, organized column-wise. A vector + * is a multivector of size 1. + * - There is a mirroring memory approach: the class has DEVICE and HOST + * data pointers. If needed, only one (or none) can be used or allocated. + * Unless triggered directly or by other function, the data is NOT + * automatically updated between HOST and DEVICE. + * - Constructor DOES NOT allocate memory. This has to be done separately. + * - You can get (and set) "raw" data easily, if needed. + * - There is memory ownership utility - vector can own memory (separate + * flags for HOST and DEVICE) or not, depending on how it is used. + * + * @author Kasia Swirydowicz + * @author Slaven Peles + */ + template + class Vector + { + public: + Vector() = default; + Vector(IdxT n); + Vector(IdxT n, IdxT k); + ~Vector(); + + Vector(const Vector&) = delete; + Vector(Vector&&) = delete; + Vector& operator=(const Vector&) = delete; + Vector& operator=(Vector&&) = delete; + + int copyFromExternal(const ScalarT* source, + memory::MemorySpace memspaceIn = memory::HOST, + memory::MemorySpace memspaceOut = memory::HOST); + int copyFromExternal(const Vector& source, + memory::MemorySpace memspaceIn = memory::HOST, + memory::MemorySpace memspaceOut = memory::HOST); + + ScalarT* getData(memory::MemorySpace memspace = memory::HOST); + ScalarT* getData(IdxT i, memory::MemorySpace memspace = memory::HOST); + const ScalarT* getData(memory::MemorySpace memspace = memory::HOST) const; + const ScalarT* getData(IdxT i, memory::MemorySpace memspace = memory::HOST) const; + + IdxT getCapacity() const; + IdxT getSize() const; + IdxT getNumVectors() const; + + int setDataUpdated(memory::MemorySpace memspace = memory::HOST); + int setDataUpdated(IdxT j, memory::MemorySpace memspace = memory::HOST); + int setData(ScalarT* data, memory::MemorySpace memspace = memory::HOST); + int allocate(memory::MemorySpace memspace = memory::HOST); + int setToZero(memory::MemorySpace memspace = memory::HOST); + int setToZero(IdxT i, memory::MemorySpace memspace = memory::HOST); + int setToConst(ScalarT C, memory::MemorySpace memspace = memory::HOST); + int setToConst(IdxT i, ScalarT C, memory::MemorySpace memspace = memory::HOST); + int syncData(memory::MemorySpace memspaceOut = memory::HOST); + int syncData(IdxT j, memory::MemorySpace memspaceOut = memory::HOST); + int resize(IdxT new_n_current); + int copyToExternal(ScalarT* dest, + IdxT i, + memory::MemorySpace memspaceSrc = memory::HOST, + memory::MemorySpace memspaceDst = memory::HOST); + int copyToExternal(ScalarT* dest, + memory::MemorySpace memspaceSrc = memory::HOST, + memory::MemorySpace memspaceDst = memory::HOST); + + private: + void setHostUpdated(bool is_updated); + void setDeviceUpdated(bool is_updated); + + IdxT n_capacity_{0}; ///< vector capacity + IdxT k_{0}; ///< number of vectors in multivector + IdxT n_size_{0}; ///< actual size of the vector + ScalarT* d_data_{nullptr}; ///< DEVICE data array + ScalarT* h_data_{nullptr}; ///< HOST data array + bool* gpu_updated_{nullptr}; ///< DEVICE data flags (updated or not) + bool* cpu_updated_{nullptr}; ///< HOST data flags (updated or not) + + 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 + }; + } // namespace LinearAlgebra +} // namespace GridKit diff --git a/tests/UnitTests/LinearAlgebra/CMakeLists.txt b/tests/UnitTests/LinearAlgebra/CMakeLists.txt index 58e28ddcb..801a7b3b7 100644 --- a/tests/UnitTests/LinearAlgebra/CMakeLists.txt +++ b/tests/UnitTests/LinearAlgebra/CMakeLists.txt @@ -1 +1,2 @@ add_subdirectory(SparseMatrix) +add_subdirectory(Vector) diff --git a/tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt b/tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt new file mode 100644 index 000000000..c97de6e54 --- /dev/null +++ b/tests/UnitTests/LinearAlgebra/Vector/CMakeLists.txt @@ -0,0 +1,7 @@ +add_executable(test_vector runVectorTests.cpp) +target_link_libraries(test_vector GridKit::dense_vector GridKit::testing) + +add_test(NAME VectorTest COMMAND $) + +install(TARGETS test_vector + RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp b/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp new file mode 100644 index 000000000..87d0e7c50 --- /dev/null +++ b/tests/UnitTests/LinearAlgebra/Vector/VectorTests.hpp @@ -0,0 +1,436 @@ +#pragma once +#include +#include +#include +#include +#include +#include + +#include +// #include +// #include +// #include +#include +#include + +namespace GridKit +{ + namespace Testing + { + using namespace LinearAlgebra; + + /** + * @class Tests for vector operations. + */ + template + class VectorTests + { + public: + VectorTests(memory::MemorySpace memspace = memory::HOST) + : memspace_(memspace) + { + } + + virtual ~VectorTests() + { + } + + /** + * @brief Test vector constructor with specified size and number of vectors. + * + * @param[in] N Number of elements in the vector. + * @param[in] k Number of vectors in the multivector. + * @return TestOutcome indicating success or failure of the test. + */ + TestOutcome vectorConstructor(IdxT N, IdxT k) + { + TestStatus status; + status = true; + + Vector x(N, k); + + if (x.getCapacity() != N) + { + std::cout << "The capacity of the vector is " << x.getCapacity() + << ", expected: " << N << "\n"; + status *= false; + } + + if (x.getSize() != N) + { + std::cout << "The size of the vector is " << x.getSize() + << ", expected: " << N << "\n"; + status *= false; + } + + if (x.getNumVectors() != k) + { + std::cout << "The number of vectors in the multivector is " << x.getNumVectors() + << ", expected: " << k << "\n"; + status *= false; + } + + return status.report(__func__); + } + + /** + * @brief Test vector constructor with specified size and default number of vectors (1). + * + * @param[in] N Number of elements in the vector. + * @return TestOutcome indicating success or failure of the test. + */ + TestOutcome vectorConstructor(IdxT N) + { + return vectorConstructor(N, 1); + } + + /** + * @brief Test resizing a vector to a new size. + * + * @param[in] N Current size of the vector. + * @param[in] new_N New size to which the vector should be resized. + * @return TestOutcome indicating success or failure of the test. + */ + TestOutcome resize(IdxT N, IdxT new_N) + { + TestStatus status; + status = true; + + Vector x(N); + x.allocate(memspace_); + + x.resize(new_N); + + if (x.getSize() != new_N) + { + std::cout << "The size of the vector after resizing is " << x.getSize() + << ", expected: " << new_N << "\n"; + status *= false; + } + + return status.report(__func__); + } + + /** + * @brief Test setting data in a vector from array. + * + * @param[in] N Number of elements in the vector. + * @return TestOutcome indicating success or failure of the test. + */ + TestOutcome setData(IdxT N) + { + TestStatus status; + status = true; + + Vector x(N); + + ScalarT* data = new ScalarT[N]; + for (IdxT i = 0; i < N; ++i) + { + data[i] = 0.1 * (ScalarT) i; + } + x.setData(data, memspace_); + + const ScalarT* x_data = x.getData(memspace_); + + if (x_data == nullptr) + { + std::cout << "The data pointer is null after setting.\n"; + status *= false; + } + else + { + for (IdxT i = 0; i < N; ++i) + { + if (!isEqual(x_data[i], data[i])) + { + std::cout << "The data in the vector is incorrect at index " << i + << ", expected: " << data[i] + << ", got: " << x_data[i] << "\n"; + status *= false; + break; + } + } + } + + delete[] data; + return status.report(__func__); + } + + /** + * @brief Test copying data between vector-array and vector-vector. + * + * This creates an array, copies it to a vector in the current memory space, then + * copies it to another vector in the same memory space, and finally back to a third on the + * HOST. Then, it verifies the content of the final vector. This test only passes if all + * copies are successful. + * + * @param[in] N Number of elements in the vector. + * @return TestOutcome indicating success or failure of the test. + */ + TestOutcome copyFromExternal(IdxT N) + { + TestStatus status; + status = true; + + Vector x(N); + ScalarT* data = new ScalarT[N]; + for (IdxT i = 0; i < N; ++i) + { + data[i] = 0.1 * (ScalarT) i; + } + + // array -> memspace + x.allocate(memspace_); + x.copyFromExternal(data, memory::HOST, memspace_); + + // memspace -> memspace + Vector y(N); + y.allocate(memspace_); + y.copyFromExternal(x, memspace_, memspace_); + + // memspace -> host + Vector z(N); + z.allocate(memory::HOST); + z.copyFromExternal(y, memspace_, memory::HOST); + + const ScalarT* z_data = z.getData(memory::HOST); + + if (z_data == nullptr) + { + std::cout << "The data pointer is null after copying from vector.\n"; + status *= false; + } + else + { + for (IdxT i = 0; i < N; ++i) + { + if (!isEqual(z_data[i], data[i])) + { + std::cout << "The data in the copied vector is incorrect at index " << i + << ", expected: " << data[i] + << ", got: " << z_data[i] << "\n"; + status *= false; + break; + } + } + } + + delete[] data; + return status.report(__func__); + } + + /** + * @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 + * to HOST for verification. + * + * @param[in] N Number of elements in the vector. + * @return TestOutcome indicating success or failure of the test. + */ + TestOutcome copyToExternal(IdxT N) + { + TestStatus status = true; + + Vector x(N); + ScalarT* data = new ScalarT[N]; + for (int i = 0; i < N; ++i) + { + data[i] = 0.1 * (ScalarT) i; + } + + x.allocate(memspace_); + x.copyFromExternal(data, memory::HOST, memspace_); + + // Copy data to an array on current memspace + ScalarT* dest = new ScalarT[N]; + // second argument is source, third is destination + x.copyToExternal(dest, memspace_, memspace_); + + // Copy to host to verify + ScalarT* dest_h = new ScalarT[N]; + if (memspace_ == memory::DEVICE) + { + mh_.copyArrayDeviceToHost(dest_h, dest, N); + dest = dest_h; + } + else + { + // If we are on HOST, we can use dest directly + delete[] dest_h; + dest_h = dest; + } + + // Verify the copied data + for (int i = 0; i < N; ++i) + { + if (!isEqual(dest_h[i], data[i])) + { + std::cout << "The data in the destination array is incorrect at index " << i + << ", expected: " << data[i] + << ", got: " << dest_h[i] << "\n"; + status *= false; + break; + } + } + + delete[] data; + delete[] dest; + return status.report(__func__); + } + + /** + * @brief Test setting all elements of a vector to a constant value. + * + * @param[in] N Number of elements in the vector. + * @return TestOutcome indicating success or failure of the test. + */ + TestOutcome setToConst(IdxT N) + { + constexpr ScalarT ONE = 1.0; + constexpr ScalarT ZERO = 0.0; + + TestStatus success = true; + + IdxT vector_size = N; + IdxT number_vectors = 3; + + Vector x(vector_size, number_vectors); + x.allocate(memspace_); + if (memspace_ == memory::DEVICE) + x.allocate(memory::HOST); + + x.setToZero(memspace_); + success *= verifyAnswer(x, ZERO); + + x.setToConst(1, ONE, memspace_); // set vector 1 to ones + success *= verifyAnswer(vector_size, x.getData(0, memspace_), ZERO); + success *= verifyAnswer(vector_size, x.getData(1, memspace_), ONE); + success *= verifyAnswer(vector_size, x.getData(2, memspace_), ZERO); + + x.setToConst(ONE, memspace_); + success *= verifyAnswer(x, ONE); + + x.setToZero(1, memspace_); // set vector 1 to zeros + success *= verifyAnswer(vector_size, x.getData(0, memspace_), ONE); + success *= verifyAnswer(vector_size, x.getData(1, memspace_), ZERO); + success *= verifyAnswer(vector_size, x.getData(2, memspace_), ONE); + + return success.report(__func__); + } + + /** + * @brief Test syncing data between HOST and DEVICE memory spaces. + * + * Creates a vector allocated in the specified memory space, then sync + * to the other memory space and verify that the data is synced correctly. + * + * @param[in] N Number of elements in the vector. + * @return TestOutcome returns a report on the test. + */ + TestOutcome syncData(IdxT N = 4) + { + constexpr ScalarT ONE = 1.0; + constexpr ScalarT ZERO = 0.0; + + TestStatus success; + success = true; + + if (memspace_ == memory::HOST) + { + return success.report(__func__); + } + + IdxT vector_size = N; + IdxT number_vectors = 3; + + Vector x(vector_size, number_vectors); + x.allocate(memspace_); + x.allocate(memory::HOST); + + // Set all vectors in x on device to ones + x.setToConst(ONE, memspace_); + // Sync host (all ones on the host, as well) + x.syncData(memory::HOST); + // Set vector 1 to all zeros on host + x.setToZero(1, memory::HOST); + // Sync vector 1 on device + x.syncData(1, memspace_); + + // Check what we have on device now is correct + success *= verifyAnswer(vector_size, x.getData(0, memspace_), ONE); + success *= verifyAnswer(vector_size, x.getData(1, memspace_), ZERO); + success *= verifyAnswer(vector_size, x.getData(2, memspace_), ONE); + + return success.report(__func__); + } + + private: + memory::MemorySpace memspace_{memory::HOST}; + MemoryHandler mh_; + + /// Check if vector elements are set to the same number + 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) + { + // std::cout << x->getData("cpu")[i] << "\n"; + 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; + } + + /// Check if an array elements are set to the same number + bool verifyAnswer(IdxT size, const ScalarT* data, ScalarT answer) + { + bool success = true; + ScalarT* x = nullptr; + + // If the data is on device copy it to the host + if (memspace_ == memory::DEVICE) + { + mh_.allocateArrayOnHost(&x, size); + mh_.copyArrayDeviceToHost(x, data, size); + // Set `data` to point to the host copy + data = x; + } + + for (size_t i = 0; i < static_cast(size); ++i) + { + if (!isEqual(data[i], answer)) + { + std::cout << std::setprecision(16); + success = false; + std::cout << "Solution vector element x[" << i << "] = " << data[i] + << ", expected: " << answer << "\n"; + break; + } + } + + if (memspace_ == memory::DEVICE) + { + mh_.deleteOnHost(x); + data = nullptr; + } + + return success; + } + }; // class + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp b/tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp new file mode 100644 index 000000000..06336eb2d --- /dev/null +++ b/tests/UnitTests/LinearAlgebra/Vector/runVectorTests.cpp @@ -0,0 +1,52 @@ +#include +#include +#include + +#include "VectorTests.hpp" + +int main(int, char**) +{ + GridKit::Testing::TestingResults result; + + { + GridKit::Testing::VectorTests test; + + std::cout << "Running vector tests on CPU:\n"; + result += test.vectorConstructor(50, 5); + result += test.vectorConstructor(50); + + result += test.setData(50); + + result += test.copyFromExternal(50); + // result += test.copyToExternal(50); + + result += test.resize(100, 50); + + result += test.setToConst(50); + result += test.setToConst(50); + } + +#ifdef GRIDKIT_USE_GPU + { + GridKit::Testing::VectorTests test(GridKit::memory::DEVICE); + + std::cout << "Running Testing on GPU:\n"; + result += test.vectorConstructor(50, 5); + result += test.vectorConstructor(50); + + result += test.setData(50); + + result += test.copyFromExternal(50); + // result += test.copyToExternal(50); + + result += test.resize(100, 50); + + result += test.setToConst(50); + result += test.setToConst(50); + result += test.syncData(50); + result += test.syncData(50); + } +#endif + + return result.summary(); +}