diff --git a/CHANGELOG.md b/CHANGELOG.md index eb388b870..e2753c5c6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -44,6 +44,7 @@ - Minor performance improvements to residual evaluation in PowerElectronics module. - Added full support for sparse Jacobians obtained with Enzyme in PhasorDynamics. - Added `Node` class to the PowerElectronics module to separate nodes from circuit components. +- Refactor Jacobian assembly in `PowerElectronics` module to reuse the CSR pattern. ## v0.1 - Refactored code to support adding different model families. diff --git a/GridKit/LinearAlgebra/CMakeLists.txt b/GridKit/LinearAlgebra/CMakeLists.txt index d290617be..3d324d485 100644 --- a/GridKit/LinearAlgebra/CMakeLists.txt +++ b/GridKit/LinearAlgebra/CMakeLists.txt @@ -1,3 +1,9 @@ add_subdirectory(SparseMatrix) add_subdirectory(DenseMatrix) + +install( + FILES + ${CMAKE_CURRENT_SOURCE_DIR}/MemoryUtils.hpp + DESTINATION + include/GridKit/LinearAlgebra) diff --git a/GridKit/LinearAlgebra/SparseMatrix/CMakeLists.txt b/GridKit/LinearAlgebra/SparseMatrix/CMakeLists.txt index 56d5adc13..14c11bebf 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/CMakeLists.txt +++ b/GridKit/LinearAlgebra/SparseMatrix/CMakeLists.txt @@ -1,6 +1,8 @@ -gridkit_add_library(CsrMatrix +gridkit_add_library(SparseMatrix SOURCES + CooMatrix.cpp CsrMatrix.cpp HEADERS COO_Matrix.hpp + CooMatrix.hpp CsrMatrix.hpp) diff --git a/GridKit/LinearAlgebra/SparseMatrix/CooMatrix.cpp b/GridKit/LinearAlgebra/SparseMatrix/CooMatrix.cpp new file mode 100644 index 000000000..f99521df3 --- /dev/null +++ b/GridKit/LinearAlgebra/SparseMatrix/CooMatrix.cpp @@ -0,0 +1,606 @@ + +#include "CooMatrix.hpp" + +#include +#include +#include +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + template + CooMatrix::CooMatrix() + { + } + + /** + * @brief basic constructor. It DOES NOT allocate any memory! + * + * @param[in] n - number of rows + * @param[in] m - number of columns + * @param[in] nnz_ - number of non-zeros + */ + template + CooMatrix::CooMatrix(IdxT n, + IdxT m, + IdxT nnz) + : n_{n}, + m_{m}, + nnz_{nnz} + { + setNotUpdated(); + + // set everything to nullptr + h_row_data_ = nullptr; + h_col_data_ = nullptr; + h_val_data_ = nullptr; + + d_row_data_ = nullptr; + d_col_data_ = nullptr; + d_val_data_ = nullptr; + + owns_cpu_sparsity_pattern_ = false; + owns_cpu_values_ = false; + + owns_gpu_sparsity_pattern_ = false; + owns_gpu_values_ = false; + } + + /** + * @brief Hijacking constructor + * + * @param[in] n + * @param[in] m + * @param[in] nnz_ + * @param[in,out] rows + * @param[in,out] cols + * @param[in,out] vals + * @param[in] memspace_src + * @param[in] memspace_dst + */ + template + CooMatrix::CooMatrix(IdxT n, + IdxT m, + IdxT nnz_, + IdxT** rows, + IdxT** cols, + RealT** vals, + memory::MemorySpace memspace_src, + memory::MemorySpace memspace_dst) + : CooMatrix(n, m, nnz_) + { + int control = -1; + if ((memspace_src == memory::HOST) && (memspace_dst == memory::HOST)) + { + control = 0; + } + if ((memspace_src == memory::HOST) && (memspace_dst == memory::DEVICE)) + { + control = 1; + } + if ((memspace_src == memory::DEVICE) && (memspace_dst == memory::HOST)) + { + control = 2; + } + if ((memspace_src == memory::DEVICE) && (memspace_dst == memory::DEVICE)) + { + control = 3; + } + + switch (control) + { + case 0: // cpu->cpu + // Set host data + h_row_data_ = *rows; + h_col_data_ = *cols; + h_val_data_ = *vals; + h_data_updated_ = true; + owns_cpu_sparsity_pattern_ = true; + owns_cpu_values_ = true; + // Set device data to null + if (d_row_data_ || d_col_data_ || d_val_data_) + { + std::cerr << "Device data unexpectedly allocated. " + << "Possible bug in matrix::Sparse class.\n"; + } + d_row_data_ = nullptr; + d_col_data_ = nullptr; + d_val_data_ = nullptr; + d_data_updated_ = false; + owns_gpu_sparsity_pattern_ = false; + owns_gpu_values_ = false; + // Hijack data from the source + *rows = nullptr; + *cols = nullptr; + *vals = nullptr; + break; + case 2: // gpu->cpu + // Set device data and copy it to host + d_row_data_ = *rows; + d_col_data_ = *cols; + d_val_data_ = *vals; + d_data_updated_ = true; + owns_gpu_sparsity_pattern_ = true; + owns_gpu_values_ = true; + syncData(memspace_dst); + // Hijack data from the source + *rows = nullptr; + *cols = nullptr; + *vals = nullptr; + break; + case 1: // cpu->gpu + // Set host data and copy it to device + h_row_data_ = *rows; + h_col_data_ = *cols; + h_val_data_ = *vals; + h_data_updated_ = true; + owns_cpu_sparsity_pattern_ = true; + owns_cpu_values_ = true; + syncData(memspace_dst); + + // Hijack data from the source + *rows = nullptr; + *cols = nullptr; + *vals = nullptr; + break; + case 3: // gpu->gpu + // Set device data + d_row_data_ = *rows; + d_col_data_ = *cols; + d_val_data_ = *vals; + d_data_updated_ = true; + owns_gpu_sparsity_pattern_ = true; + owns_gpu_values_ = true; + // Set host data to null + if (h_row_data_ || h_col_data_ || h_val_data_) + { + std::cerr << "Host data unexpectedly allocated. " + << "Possible bug in matrix::Sparse class.\n"; + } + h_row_data_ = nullptr; + h_col_data_ = nullptr; + h_val_data_ = nullptr; + h_data_updated_ = false; + owns_cpu_sparsity_pattern_ = false; + owns_cpu_values_ = false; + // Hijack data from the source + *rows = nullptr; + *cols = nullptr; + *vals = nullptr; + break; + default: + std::cerr << "CooMatrix constructor failed! " + << "Possible bug in memory spaces setting.\n"; + break; + } + } + + template + CooMatrix::~CooMatrix() + { + this->destroyMatrixData(memory::HOST); + this->destroyMatrixData(memory::DEVICE); + delete[] map_to_sorted_; + delete[] map_to_dedup_; + } + + /** + * @brief set the matrix update flags to false (for both HOST and DEVICE). + */ + template + void CooMatrix::setNotUpdated() + { + h_data_updated_ = false; + d_data_updated_ = false; + } + + /** + * @brief get number of matrix rows + * + * @return number of matrix rows. + */ + template + IdxT CooMatrix::getNumRows() const + { + return this->n_; + } + + /** + * @brief get number of matrix columns + * + * @return number of matrix columns. + */ + template + IdxT CooMatrix::getNumColumns() const + { + return this->m_; + } + + /** + * @brief get number of non-zeros in the matrix. + * + * @return number of non-zeros. + */ + template + IdxT CooMatrix::getNnz() const + { + return this->nnz_; + } + + /** + * @brief Set the pointers for matrix row, column, value data. + * + * Useful if interfacing with other codes - this function only assigns + * pointers, but it does not allocate nor copy anything. The data ownership + * flags would be set to false (default). + * + * @param[in] row_data - pointer to row data (array of integers) + * @param[in] col_data - pointer to column data (array of integers) + * @param[in] val_data - pointer to value data (array of real numbers) + * @param[in] memspace - memory space (HOST or DEVICE) of incoming data + * + * @return 0 if successful, 1 if not. + */ + template + int CooMatrix::setDataPointers(IdxT* row_data, + IdxT* col_data, + RealT* val_data, + memory::MemorySpace memspace) + { + using namespace memory; + + setNotUpdated(); + + switch (memspace) + { + case HOST: + if (owns_cpu_sparsity_pattern_ && (h_row_data_ || h_col_data_)) + { + std::cerr << "Trying to set matrix host data, but the data already set!\n"; + std::cerr << "Ignoring setDataPointers function call ...\n"; + return 1; + } + if (owns_cpu_values_ && h_val_data_) + { + std::cerr << "Trying to set matrix host values, but the values already set!\n"; + std::cerr << "Ignoring setValuesPointer function call ...\n"; + return 1; + } + h_row_data_ = row_data; + h_col_data_ = col_data; + h_val_data_ = val_data; + h_data_updated_ = true; + owns_cpu_sparsity_pattern_ = false; + owns_cpu_values_ = false; + break; + case DEVICE: + if (owns_gpu_sparsity_pattern_ && (d_row_data_ || d_col_data_)) + { + std::cerr << "Trying to set matrix host data, but the data already set!\n"; + std::cerr << "Ignoring setDataPointers function call ...\n"; + return 1; + } + if (owns_gpu_values_ && d_val_data_) + { + std::cerr << "Trying to set matrix device values, but the values already set!\n"; + std::cerr << "Ignoring setValuesPointer function call ...\n"; + return 1; + } + d_row_data_ = row_data; + d_col_data_ = col_data; + d_val_data_ = val_data; + d_data_updated_ = true; + owns_gpu_sparsity_pattern_ = false; + owns_gpu_values_ = false; + break; + } + return 0; + } + + /** + * @brief destroy matrix data (HOST or DEVICE) if the matrix owns it + * (will attempt to destroy all three arrays). + * + * @param[in] memspace - memory space (HOST or DEVICE) of incoming data + * + * @return 0 if successful, -1 if not. + * + */ + template + int CooMatrix::destroyMatrixData(memory::MemorySpace memspace) + { + using namespace memory; + switch (memspace) + { + case HOST: + if (owns_cpu_sparsity_pattern_) + { + delete[] h_row_data_; + delete[] h_col_data_; + h_row_data_ = nullptr; + h_col_data_ = nullptr; + } + if (owns_cpu_values_) + { + delete[] h_val_data_; + h_val_data_ = nullptr; + } + return 0; + case DEVICE: + if (owns_gpu_sparsity_pattern_) + { + mem_.deleteOnDevice(d_row_data_); + mem_.deleteOnDevice(d_col_data_); + d_row_data_ = nullptr; + d_col_data_ = nullptr; + } + if (owns_gpu_values_) + { + mem_.deleteOnDevice(d_val_data_); + d_val_data_ = nullptr; + } + return 0; + default: + return -1; + } + } + + /** + * @brief Sort COO entries by (row, col), merge duplicate + * entries, and build CSR row pointers. + * + * @note This function operates on host (CPU) data only. + * + * @return Pointer to CSR row pointer array + */ + template + IdxT* CooMatrix::getCsrRowData() + { + if (!h_data_updated_) + { + std::cerr << "CooMatrix::getCsrRowData requires up-to-date host data, but host is out of date!\n"; + assert(h_data_updated_); + return nullptr; + } + + map_to_sorted_ = new IdxT[static_cast(nnz_)]; + + for (IdxT i = 0; i < nnz_; i++) + { + map_to_sorted_[i] = i; + } + + // Sort by (row, col) + std::sort(map_to_sorted_, map_to_sorted_ + nnz_, [this](IdxT a, IdxT b) + { if (h_row_data_[a] != h_row_data_[b]) return h_row_data_[a] < h_row_data_[b]; + return h_col_data_[a] < h_col_data_[b]; }); + + IdxT* row_data_tmp = new IdxT[static_cast(nnz_)]; + IdxT* col_data_tmp = new IdxT[static_cast(nnz_)]; + RealT* val_data_tmp = new RealT[static_cast(nnz_)]; + + for (IdxT i = 0; i < nnz_; i++) + { + row_data_tmp[i] = h_row_data_[map_to_sorted_[i]]; + col_data_tmp[i] = h_col_data_[map_to_sorted_[i]]; + val_data_tmp[i] = h_val_data_[map_to_sorted_[i]]; + } + for (IdxT i = 0; i < nnz_; i++) + { + h_row_data_[i] = row_data_tmp[i]; + h_col_data_[i] = col_data_tmp[i]; + h_val_data_[i] = val_data_tmp[i]; + } + delete[] row_data_tmp; + delete[] col_data_tmp; + delete[] val_data_tmp; + + map_to_dedup_ = new IdxT[static_cast(nnz_)]; + IdxT* csr_row_data = new IdxT[static_cast(n_ + 1)](); + + map_to_dedup_[0] = 0; + csr_row_data[h_row_data_[0] + 1]++; + + // Deduplicate sorted entries + IdxT w = 0; + for (IdxT i = 1; i < nnz_; i++) + { + if (h_row_data_[i] == h_row_data_[w] && h_col_data_[i] == h_col_data_[w]) + { + h_val_data_[w] += h_val_data_[i]; + map_to_dedup_[i] = w; + } + else + { + ++w; + h_row_data_[w] = h_row_data_[i]; + h_col_data_[w] = h_col_data_[i]; + h_val_data_[w] = h_val_data_[i]; + map_to_dedup_[i] = w; + csr_row_data[h_row_data_[w] + 1]++; + } + } + // nnz after deduplication + nnz_ = w + 1; + + for (IdxT i = 0; i < n_; i++) + { + csr_row_data[i + 1] += csr_row_data[i]; + } + + return csr_row_data; + } + + template + const IdxT* CooMatrix::getMapToSorted() const + { + return map_to_sorted_; + } + + template + const IdxT* CooMatrix::getMapToDeduplicated() const + { + return map_to_dedup_; + } + + template + IdxT* CooMatrix::getRowData(memory::MemorySpace memspace) + { + using namespace memory; + + switch (memspace) + { + case HOST: + return this->h_row_data_; + case DEVICE: + return this->d_row_data_; + default: + return nullptr; + } + } + + template + IdxT* CooMatrix::getColData(memory::MemorySpace memspace) + { + using namespace memory; + + switch (memspace) + { + case HOST: + return this->h_col_data_; + case DEVICE: + return this->d_col_data_; + default: + return nullptr; + } + } + + template + RealT* CooMatrix::getValues(memory::MemorySpace memspace) + { + using namespace memory; + + switch (memspace) + { + case HOST: + return this->h_val_data_; + case DEVICE: + return this->d_val_data_; + default: + return nullptr; + } + } + + /** + * @brief Sync data in memspace with the updated memory space. + * + * @param memspace - memory space to be synced up (HOST or DEVICE) + * @return int - 0 if successful, error code otherwise + * + * @pre The memory space other than `memspace` must be up-to-date. Otherwise, + * this function will return an error. + * + * @see CooMatrix::setUpdated + */ + template + int CooMatrix::syncData(memory::MemorySpace memspace) + { + using namespace memory; + + switch (memspace) + { + case HOST: + // check if we need to copy or not + assert(((h_row_data_ == nullptr) == (h_col_data_ == nullptr)) && "In CooMatrix::syncData one of host row or column data is null!\n"); + + if (h_data_updated_) + { + std::cerr << "CooMatrix::syncData is trying to sync host, but host already up to date!\n"; + assert(!h_data_updated_); + return 1; + } + if (!d_data_updated_) + { + std::cerr << "CooMatrix::syncData is trying to sync host with device, but device is out of date!\n" + << "See CooMatrix::syncData documentation\n."; + assert(d_data_updated_); + } + if ((h_row_data_ == nullptr) && (h_col_data_ == nullptr)) + { + h_row_data_ = new IdxT[static_cast(n_ + 1)]; + h_col_data_ = new IdxT[static_cast(nnz_)]; + owns_cpu_sparsity_pattern_ = true; + } + if (h_val_data_ == nullptr) + { + h_val_data_ = new RealT[static_cast(nnz_)]; + owns_cpu_values_ = true; + } + mem_.copyArrayDeviceToHost(h_row_data_, d_row_data_, n_ + 1); + mem_.copyArrayDeviceToHost(h_col_data_, d_col_data_, nnz_); + mem_.copyArrayDeviceToHost(h_val_data_, d_val_data_, nnz_); + h_data_updated_ = true; + return 0; + case DEVICE: + assert(((d_row_data_ == nullptr) == (d_col_data_ == nullptr)) && "In CooMatrix::syncData one of device row or column data is null!\n"); + + if (d_data_updated_) + { + std::cerr << "CooMatrix::syncData is trying to sync device, but device already up to date!\n"; + assert(!d_data_updated_); + return 1; + } + if (!h_data_updated_) + { + std::cerr << "CooMatrix::syncData is trying to sync device with host, but host is out of date!\n" + << "See CooMatrix::syncData documentation\n."; + assert(h_data_updated_); + } + if ((d_row_data_ == nullptr) && (d_col_data_ == nullptr)) + { + mem_.allocateArrayOnDevice(&d_row_data_, n_ + 1); + mem_.allocateArrayOnDevice(&d_col_data_, nnz_); + owns_gpu_sparsity_pattern_ = true; + } + if (d_val_data_ == nullptr) + { + mem_.allocateArrayOnDevice(&d_val_data_, nnz_); + owns_gpu_values_ = true; + } + mem_.copyArrayHostToDevice(d_row_data_, h_row_data_, n_ + 1); + mem_.copyArrayHostToDevice(d_col_data_, h_col_data_, nnz_); + mem_.copyArrayHostToDevice(d_val_data_, h_val_data_, nnz_); + d_data_updated_ = true; + return 0; + default: + return 1; + } // switch + } + + /** + * @brief Prints matrix data. + * + * @param out - Output stream where the matrix data is printed + */ + template + void CooMatrix::print(std::ostream& out, IdxT indexing_base) + { + out << std::scientific << std::setprecision(std::numeric_limits::digits10); + for (IdxT i = 0; i < n_; ++i) + { + for (IdxT j = h_row_data_[i]; j < h_row_data_[i + 1]; ++j) + { + out << i + indexing_base << " " + << h_col_data_[j] + indexing_base << " " + << h_val_data_[j] << "\n"; + } + } + } + + template class CooMatrix; + template class CooMatrix; + template class CooMatrix; + } // namespace LinearAlgebra +} // namespace GridKit diff --git a/GridKit/LinearAlgebra/SparseMatrix/CooMatrix.hpp b/GridKit/LinearAlgebra/SparseMatrix/CooMatrix.hpp new file mode 100644 index 000000000..b167d2760 --- /dev/null +++ b/GridKit/LinearAlgebra/SparseMatrix/CooMatrix.hpp @@ -0,0 +1,90 @@ +#pragma once + +#include +#include +#include + +#include + +namespace GridKit +{ + namespace LinearAlgebra + { + + template + class CooMatrix + { + public: + CooMatrix(); + + CooMatrix(IdxT n, IdxT m, IdxT nnz); + + CooMatrix(IdxT n, + IdxT m, + IdxT nnz, + IdxT** rows, + IdxT** cols, + RealT** vals, + memory::MemorySpace memspace_src = memory::HOST, + memory::MemorySpace memspace_dst = memory::HOST); + + ~CooMatrix(); + + // accessors + IdxT getNumRows() const; + IdxT getNumColumns() const; + IdxT getNnz() const; + + const IdxT* getMapToSorted() const; + const IdxT* getMapToDeduplicated() const; + + IdxT* getRowData(memory::MemorySpace memspace = memory::HOST); + IdxT* getColData(memory::MemorySpace memspace = memory::HOST); + RealT* getValues(memory::MemorySpace memspace = memory::HOST); + + int setDataPointers(IdxT* row_data, + IdxT* col_data, + RealT* val_data, + memory::MemorySpace memspace); + + int destroyMatrixData(memory::MemorySpace memspace); + + void print(std::ostream& file_out = std::cout, IdxT indexing_base = 0); + + int syncData(memory::MemorySpace memspace_out); + + IdxT* getCsrRowData(); + + private: + IdxT n_{0}; ///< number of rows + IdxT m_{0}; ///< number of columns + IdxT nnz_{0}; ///< number of non-zeros + + // host data + IdxT* h_row_data_{nullptr}; ///< row data (HOST) + IdxT* h_col_data_{nullptr}; ///< column data (HOST) + RealT* h_val_data_{nullptr}; ///< value data (HOST) + bool h_data_updated_{false}; ///< HOST update flag + + // gpu data + IdxT* d_row_data_{nullptr}; ///< row data (DEVICE) + IdxT* d_col_data_{nullptr}; ///< column data (DEVICE) + RealT* d_val_data_{nullptr}; ///< value data (DEVICE) + bool d_data_updated_{false}; ///< DEVICE update flag + + void setNotUpdated(); + + // Data ownership flags + bool owns_cpu_sparsity_pattern_{false}; ///< for row/col data + bool owns_cpu_values_{false}; ///< for nonzero values + + bool owns_gpu_sparsity_pattern_{false}; ///< for row/col data + bool owns_gpu_values_{false}; ///< for nonzero values + + 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 + }; + } // namespace LinearAlgebra +} // namespace GridKit diff --git a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp index a1296c570..87fa7433d 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp +++ b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.cpp @@ -820,5 +820,6 @@ namespace GridKit template class CsrMatrix; template class CsrMatrix; + template class CsrMatrix; } // namespace LinearAlgebra } // namespace GridKit diff --git a/GridKit/Model/Evaluator.hpp b/GridKit/Model/Evaluator.hpp index 93c812ffc..ffec10cf2 100644 --- a/GridKit/Model/Evaluator.hpp +++ b/GridKit/Model/Evaluator.hpp @@ -4,6 +4,7 @@ #include #include +#include #include #include @@ -19,8 +20,9 @@ namespace GridKit class Evaluator { public: - using RealT = typename GridKit::ScalarTraits::RealT; - using MatrixT = GridKit::LinearAlgebra::COO_Matrix; //\todo Use CsrMatrix + using RealT = typename GridKit::ScalarTraits::RealT; + using MatrixT = GridKit::LinearAlgebra::COO_Matrix; //\todo Use CsrMatrix + using CsrMatrix = GridKit::LinearAlgebra::CsrMatrix; Evaluator() { @@ -82,6 +84,16 @@ namespace GridKit { } + /** + * @brief Return a pointer to the CSR Jacobian + * + * @todo Remove this and use CsrMatirx for jac_ + */ + virtual CsrMatrix* getCsrJacobian() const + { + return nullptr; + } + /** * @brief Is the Jacobian defined. Used in IDA to determine wether DQ is used or not * diff --git a/GridKit/Model/PowerElectronics/CircuitComponent.hpp b/GridKit/Model/PowerElectronics/CircuitComponent.hpp index cbf8a6135..d83b7a078 100644 --- a/GridKit/Model/PowerElectronics/CircuitComponent.hpp +++ b/GridKit/Model/PowerElectronics/CircuitComponent.hpp @@ -18,8 +18,9 @@ namespace GridKit class CircuitComponent : public Model::Evaluator { public: - using RealT = typename Model::Evaluator::RealT; - using MatrixT = typename Model::Evaluator::MatrixT; + using RealT = typename Model::Evaluator::RealT; + using MatrixT = typename Model::Evaluator::MatrixT; + using CsrMatrix = typename Model::Evaluator::CsrMatrix; CircuitComponent() = default; diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index f5796fc88..634409227 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -2,14 +2,18 @@ #pragma once +#include #include #include #include #include #include +#include +#include +#include #include -#include +#include #include namespace GridKit @@ -60,7 +64,10 @@ namespace GridKit class PowerElectronicsModel : public CircuitComponent { using RealT = typename CircuitComponent::RealT; + using MatrixT = typename CircuitComponent::MatrixT; + using CsrMatrix = typename CircuitComponent::CsrMatrix; using component_type = CircuitComponent; + using node_type = CircuitNode; using CircuitComponent::size_; using CircuitComponent::nnz_; @@ -123,7 +130,11 @@ namespace GridKit virtual ~PowerElectronicsModel() { for (auto comp : this->components_) + { delete comp; + } + delete csr_jac_; + delete[] map_to_csr_; } /** @@ -162,12 +173,13 @@ namespace GridKit } /** - * @brief Allocate the vector data with size amount - * @todo Add capability to go through component model connection to get the size of the actual vector + * @brief Allocate system vectors and construct the system CSR Jacobian * - * @param[in] s size of the vector + * @param[in] s size of the vector (total number of unknowns) * * @post System model vectors allocated with size s + * @post CSR Jacobian sparsity pattern is computed + * @post COO->CSR mapping is computed * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ @@ -185,6 +197,81 @@ namespace GridKit yp_.resize(size_); f_.resize(size_); + // Evaluate component Jacobians to get sparsity + distributeVectors(); + for (const auto& component : components_) + { + component->evaluateJacobian(); + } + + // Count the number of non-zeros + IdxT nnz_dup = 0; + for (const auto& component : components_) + { + std::tuple&, std::vector&, std::vector&> entries = component->getJacobian().getEntries(); + const auto& [r, c, v] = entries; + + for (IdxT i = 0; i < static_cast(r.size()); ++i) + { + if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) + { + ++nnz_dup; + } + } + } + + // Allocate COO triplet arrays (we own these until we hand off to CsrMatrix) + IdxT* rows_dup = new IdxT[nnz_dup]; + IdxT* cols_dup = new IdxT[nnz_dup]; + RealT* vals_dup = new RealT[nnz_dup]; + + IdxT counter = 0; + for (const auto& component : components_) + { + std::tuple&, std::vector&, std::vector&> entries = component->getJacobian().getEntries(); + const auto& [r, c, v] = entries; + + for (IdxT i = 0; i < static_cast(r.size()); ++i) + { + if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) + { + rows_dup[counter] = component->getNodeConnection(r[i]); + cols_dup[counter] = component->getNodeConnection(c[i]); + vals_dup[counter] = v[i]; + counter++; + } + } + } + + // Build the system COO Jacobian + LinearAlgebra::CooMatrix jac(size_, size_, nnz_dup, &rows_dup, &cols_dup, &vals_dup); + + // Populate CSR data with sort and deduplicate + IdxT* row_ptrs = jac.getCsrRowData(); + + // Deduplicated nnz + nnz_ = jac.getNnz(); + + // Allocate cols/vals with deduplicated nnz + IdxT* cols = new IdxT[nnz_]; + RealT* vals = new RealT[nnz_]; + + std::copy(jac.getColData(), jac.getColData() + nnz_, cols); + std::copy(jac.getValues(), jac.getValues() + nnz_, vals); + + // Create the CSR Jacobian + csr_jac_ = new CsrMatrix(size_, size_, nnz_, &row_ptrs, &cols, &vals); + + const IdxT* map_to_sorted = jac.getMapToSorted(); + const IdxT* map_to_dedup = jac.getMapToDeduplicated(); + + // Build a mappping from original COO index to CSR index + map_to_csr_ = new IdxT[nnz_dup]; + for (IdxT i = 0; i < nnz_dup; ++i) + { + map_to_csr_[map_to_sorted[i]] = map_to_dedup[i]; + } + return 0; } @@ -237,7 +324,7 @@ namespace GridKit return 0; } - int tagDifferentiable() + int tagDifferentiable() final { return 0; } @@ -282,46 +369,44 @@ namespace GridKit } /** - * @brief Creates the Sparse COO Jacobian representing \alpha dF/dy' + dF/dy + * @brief Creates the system Jacobian representing \alpha dF/dy' + dF/dy + * + * Updates the CSR Jacobian values using the per-component mappings + * computed during allocate(). * * @return int 0 if successful, positive if there's a recoverable error, negative if unrecoverable */ int evaluateJacobian() final { - jac_.zeroMatrix(); distributeVectors(); - // Evaluate component jacs + // Zero out values + RealT* vals = csr_jac_->getValues(); + for (IdxT i = 0; i < csr_jac_->getNnz(); ++i) + { + vals[i] = 0.0; + } + + // Update CSR values from component Jacobians + IdxT counter = 0; for (const auto& component : components_) { component->evaluateJacobian(); - // get references to local jacobian - std::tuple&, std::vector&, std::vector&> tpm = component->getJacobian().getEntries(); - const auto& [r, c, v] = tpm; + std::tuple&, std::vector&, std::vector&> entries = component->getJacobian().getEntries(); + const auto& [r, c, v] = entries; - // Create copies of data to handle groundings - std::vector rgr; - std::vector cgr; - std::vector vgr; - for (IdxT i = 0; i < static_cast(r.size()); i++) + for (IdxT i = 0; i < static_cast(r.size()); ++i) { if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) { - rgr.push_back(component->getNodeConnection(r[i])); - cgr.push_back(component->getNodeConnection(c[i])); - vgr.push_back(v[i]); + vals[map_to_csr_[counter]] += v[i]; + ++counter; } } - - // AXPY to Global Jacobian - // elementwise jac_(rgr, cgr) += vgr - jac_.axpy(1.0, rgr, cgr, vgr); } - // jac_.printMatrixMarket("ScaleMicrogrid_Jacobian_N2_number" + std::to_string(jac_call_count_) + ".mtx", "Jacobian N2 number " + std::to_string(jac_call_count_)); jac_call_count_++; - return 0; } @@ -330,7 +415,6 @@ namespace GridKit */ int evaluateIntegrand() final { - return 0; } @@ -398,12 +482,16 @@ namespace GridKit * @param[in] filename * @param[in] title */ - void printJacobianMatrixMarket(std::string filename, std::string title) { jac_.printMatrixMarket(filename, title); } + CsrMatrix* getCsrJacobian() const override + { + return csr_jac_; + } + void addComponent(component_type* component) { components_.push_back(component); @@ -414,6 +502,9 @@ namespace GridKit std::vector components_; + IdxT* map_to_csr_{nullptr}; + CsrMatrix* csr_jac_{nullptr}; + int jac_call_count_{0}; bool use_jac_; diff --git a/GridKit/Solver/Dynamic/CMakeLists.txt b/GridKit/Solver/Dynamic/CMakeLists.txt index 0cceed1fd..4e4abd1e4 100644 --- a/GridKit/Solver/Dynamic/CMakeLists.txt +++ b/GridKit/Solver/Dynamic/CMakeLists.txt @@ -20,7 +20,8 @@ if (GRIDKIT_ENABLE_SUNDIALS_SPARSE) PUBLIC SUNDIALS::idas PUBLIC SUNDIALS::sunlinsolklu PUBLIC GridKit::definitions - PUBLIC GridKit::utilities_logger) + PUBLIC GridKit::utilities_logger + PUBLIC GridKit::SparseMatrix) else() gridkit_add_library(solvers_dyn SOURCES @@ -31,5 +32,6 @@ else() PUBLIC SUNDIALS::nvecserial PUBLIC SUNDIALS::idas PUBLIC GridKit::definitions - PUBLIC GridKit::utilities_logger) + PUBLIC GridKit::utilities_logger + PUBLIC GridKit::SparseMatrix) endif() diff --git a/GridKit/Solver/Dynamic/Ida.cpp b/GridKit/Solver/Dynamic/Ida.cpp index aec2d9b0c..cac99195c 100644 --- a/GridKit/Solver/Dynamic/Ida.cpp +++ b/GridKit/Solver/Dynamic/Ida.cpp @@ -168,8 +168,20 @@ namespace AnalysisManager { int retval = 0; - sunindextype n = static_cast(model_->size()); - sunindextype nnz = static_cast((model_->getJacobian()).nnz()); + sunindextype n = static_cast(model_->size()); + sunindextype nnz; + + /// @todo Remove this flag and require all apps to configure CsrJac + bool use_csr_jac = (model_->getCsrJacobian() != nullptr); + + if (use_csr_jac) + { + nnz = static_cast(model_->getCsrJacobian()->getNnz()); + } + else + { + nnz = static_cast((model_->getJacobian()).nnz()); + } JacobianMat_ = SUNSparseMatrix(n, n, @@ -184,7 +196,14 @@ namespace AnalysisManager retval = IDASetLinearSolver(solver_, linearSolver_, JacobianMat_); checkOutput(retval, "IDASetLinearSolver"); - retval = IDASetJacFn(solver_, this->Jac); + if (use_csr_jac) + { + retval = IDASetJacFn(solver_, this->CsrJac); + } + else + { + retval = IDASetJacFn(solver_, this->Jac); + } checkOutput(retval, "IDASetJacFn"); return retval; @@ -804,6 +823,52 @@ namespace AnalysisManager return 0; } + /** + * @brief Jacobian evaluation + * + * @note The model Jacobian is stored in CSR format. + * + * @tparam ScalarT + * @tparam IdxT + * + * @todo Update PhasorDynamics to use this implementation. + */ + template + int Ida::CsrJac(RealT t, RealT cj, N_Vector yy, N_Vector yp, N_Vector, SUNMatrix J, void* user_data, N_Vector, N_Vector, N_Vector) + { + GridKit::Model::Evaluator* model = static_cast*>(user_data); + + model->updateTime(t, cj); + copyVec(yy, model->y()); + copyVec(yp, model->yp()); + + model->evaluateJacobian(); + + using CsrMatrix = GridKit::LinearAlgebra::CsrMatrix; + CsrMatrix* Jac = model->getCsrJacobian(); + + SUNMatZero(J); + + sunindextype* sun_row_ptrs = SUNSparseMatrix_IndexPointers(J); + sunindextype* sun_cols = SUNSparseMatrix_IndexValues(J); + RealT* sun_vals = SUNSparseMatrix_Data(J); + + IdxT n = Jac->getNumRows(); + IdxT nnz = Jac->getNnz(); + + // Get reference to the jacobian entries + IdxT* row_ptrs = Jac->getRowData(); + IdxT* cols = Jac->getColData(); + RealT* vals = Jac->getValues(); + + // Copy data from model jac to sundials + std::copy(row_ptrs, row_ptrs + n + 1, sun_row_ptrs); + std::copy(cols, cols + nnz, sun_cols); + std::copy(vals, vals + nnz, sun_vals); + + return 0; + } + /** * @brief Integrand evaluation * @@ -972,7 +1037,7 @@ namespace AnalysisManager /** * @brief Accumulate another stats object into this one, allowing for stats to be kept - * across multiple simulations with IDA. + * across multiple simulations with IDA */ IdaStats& IdaStats::operator+=(const IdaStats& other) { diff --git a/GridKit/Solver/Dynamic/Ida.hpp b/GridKit/Solver/Dynamic/Ida.hpp index 20ff56aa8..0c4c09bf1 100644 --- a/GridKit/Solver/Dynamic/Ida.hpp +++ b/GridKit/Solver/Dynamic/Ida.hpp @@ -151,6 +151,17 @@ namespace AnalysisManager N_Vector tmp2, N_Vector tmp3); + static int CsrJac(RealT t, + RealT cj, + N_Vector yy, + N_Vector yp, + N_Vector resvec, + SUNMatrix J, + void* user_data, + N_Vector tmp1, + N_Vector tmp2, + N_Vector tmp3); + static int Integrand(RealT t, N_Vector yy, N_Vector yp, diff --git a/examples/LinearAlgebra/SparseTest/CMakeLists.txt b/examples/LinearAlgebra/SparseTest/CMakeLists.txt index 3ef32ea66..c4c664d5b 100644 --- a/examples/LinearAlgebra/SparseTest/CMakeLists.txt +++ b/examples/LinearAlgebra/SparseTest/CMakeLists.txt @@ -1,7 +1,7 @@ add_executable(spmattest SparseTest.cpp) -target_link_libraries(spmattest GridKit::CsrMatrix) +target_link_libraries(spmattest GridKit::SparseMatrix) add_test(NAME SparseMatrixTest COMMAND spmattest) install(TARGETS spmattest RUNTIME DESTINATION bin) diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp index 44804aa7a..c8f970cc9 100644 --- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp +++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp @@ -97,7 +97,7 @@ int main(int /* argc */, char const** /* argv */) sysmodel.updateTime(0.0, 1.0); sysmodel.evaluateJacobian(); std::cout << "Intial Jacobian with alpha = 1:\n"; - sysmodel.getJacobian().printMatrix(); + sysmodel.getCsrJacobian()->print(); // Create numerical integrator and configure it for the generator model AnalysisManager::Sundials::Ida idas(&sysmodel); diff --git a/tests/UnitTests/LinearAlgebra/SparseMatrix/CMakeLists.txt b/tests/UnitTests/LinearAlgebra/SparseMatrix/CMakeLists.txt index f94f273bd..25f371579 100644 --- a/tests/UnitTests/LinearAlgebra/SparseMatrix/CMakeLists.txt +++ b/tests/UnitTests/LinearAlgebra/SparseMatrix/CMakeLists.txt @@ -1,8 +1,12 @@ -add_executable(test_sparse runSparseTests.cpp) -target_link_libraries(test_sparse GridKit::CsrMatrix) - -add_test(NAME SparseTest COMMAND $) +add_executable(test_sparse_csr runSparseTests.cpp) +target_link_libraries(test_sparse_csr GridKit::SparseMatrix) +add_test(NAME SparseTest COMMAND $) -install(TARGETS test_sparse +add_executable(test_sparse_coo runSparseCooTests.cpp) +target_link_libraries(test_sparse_coo GridKit::SparseMatrix) + +add_test(NAME SparseCooTest COMMAND $) + +install(TARGETS test_sparse_csr test_sparse_coo RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCooTests.hpp b/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCooTests.hpp new file mode 100644 index 000000000..d784029a6 --- /dev/null +++ b/tests/UnitTests/LinearAlgebra/SparseMatrix/SparseCooTests.hpp @@ -0,0 +1,178 @@ +#include +#include + +namespace GridKit +{ + namespace Testing + { + using namespace LinearAlgebra; + + template + class SparseCooTests + { + + using CooMatrix = LinearAlgebra::CooMatrix; + + public: + SparseCooTests(memory::MemorySpace memspace = memory::HOST) + : memspace_(memspace) + { + } + + ~SparseCooTests() + { + } + + /** + * @brief Constructor test for COO sparse matrix + * + * Constructs a COO matrix with the given parameters and confirms that the + * parameters are stored correctly. + * + * @param[in] n - number of rows + * @param[in] m - number of columns + * @param[in] nnz - number of non-zeros + * + * @return TestOutcome indicating success or failure of the test + */ + TestOutcome constructor(IdxT n, IdxT m, IdxT nnz) + { + TestStatus success; + success = true; + + CooMatrix A(n, m, nnz); + + if (A.getNumRows() != n || A.getNumColumns() != m || A.getNnz() != nnz) + { + std::cout << "Matrix dimensions do not match expected values.\n"; + success = false; + } + + return success.report(__func__); + } + + /** + * @brief Test set data pointers for the COO sparse matrix + * + * Sets the row, column, and value data pointers for the COO matrix and + * checks if the pointers are set correctly. + * + * @param[in] n - number of rows + * @param[in] m - number of columns + * @param[in] row_density - number of non-zeros per row + * + * @return TestOutcome indicating success or failure of the test + */ + TestOutcome setDataPointers(IdxT n, IdxT m, IdxT row_density) + { + TestStatus success; + success = true; + + IdxT nnz = n * row_density; + + CooMatrix A(n, m, nnz); + + IdxT* h_row_data = new IdxT[static_cast(nnz)]; + for (IdxT i = 0; i < nnz; ++i) + { + h_row_data[i] = i / row_density; + } + IdxT* h_col_data = new IdxT[static_cast(nnz)]; + for (IdxT i = 0; i < nnz; ++i) + { + h_col_data[i] = i % m; + } + ScalarT* h_val_data = new ScalarT[static_cast(nnz)]; + for (IdxT i = 0; i < nnz; ++i) + { + h_val_data[i] = static_cast(i + 1); + } + + if (A.setDataPointers(h_row_data, h_col_data, h_val_data, memory::HOST) != 0) + { + std::cout << "Failed to set data pointers.\n"; + success = false; + } + else if (A.getRowData(memory::HOST) != h_row_data || A.getColData(memory::HOST) != h_col_data || A.getValues(memory::HOST) != h_val_data) + { + std::cout << "Data pointers do not point to expected values.\n"; + success = false; + } + + delete[] h_row_data; + delete[] h_col_data; + delete[] h_val_data; + + return success.report(__func__); + } + + /** + * @brief Test getCsrRowData on a small COO matrix with duplicates + * + * Builds a 3x3 COO matrix with 5 entries (including one duplicate pair), + * calls getCsrRowData(), and verifies that the returned CSR row pointers + * and the updated nnz are correct. + * + * @return TestOutcome indicating success or failure of the test + */ + TestOutcome getCsrRowData() + { + TestStatus success; + success = true; + + const IdxT n = 3; + const IdxT m = 3; + const IdxT nnz_dup = 5; + + CooMatrix A(n, m, nnz_dup); + + IdxT* rows = new IdxT[nnz_dup]{2, 0, 0, 1, 1}; + IdxT* cols = new IdxT[nnz_dup]{2, 0, 1, 1, 1}; + ScalarT* vals = new ScalarT[nnz_dup]{5, 1, 2, 3, 4}; + + if (A.setDataPointers(rows, cols, vals, memory::HOST) != 0) + { + std::cout << "Failed to set data pointers.\n"; + success = false; + } + else + { + IdxT* row_ptrs = A.getCsrRowData(); + + const IdxT expected_nnz = 4; + const IdxT expected_row_ptrs[4] = {0, 2, 3, 4}; + + if (A.getNnz() != expected_nnz) + { + std::cout << "nnz after deduplication is " << A.getNnz() + << ", expected " << expected_nnz << ".\n"; + success = false; + } + + for (IdxT i = 0; i <= n; ++i) + { + if (row_ptrs[i] != expected_row_ptrs[i]) + { + std::cout << "row_ptrs[" << i << "] = " << row_ptrs[i] + << ", expected " << expected_row_ptrs[i] << ".\n"; + success = false; + break; + } + } + + delete[] row_ptrs; + } + + delete[] rows; + delete[] cols; + delete[] vals; + + return success.report(__func__); + } + + private: + memory::MemorySpace memspace_; + MemoryHandler mem_; + }; // class SparseCooTests + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/LinearAlgebra/SparseMatrix/runSparseCooTests.cpp b/tests/UnitTests/LinearAlgebra/SparseMatrix/runSparseCooTests.cpp new file mode 100644 index 000000000..7d74afa69 --- /dev/null +++ b/tests/UnitTests/LinearAlgebra/SparseMatrix/runSparseCooTests.cpp @@ -0,0 +1,38 @@ +#include "SparseCooTests.hpp" + +using namespace GridKit; +using namespace LinearAlgebra; +using namespace Testing; + +/** + * @brief Run COO sparse matrix tests with a given backend + * + * @param[in] backend - name of the hardware backend + * @param[in] memspace - memory space for the tests + * @param[out] result - test results + */ +template +void runTests(const std::string& backend, memory::MemorySpace memspace, TestingResults& result) +{ + std::cout << "Running tests on " << backend << ":\n"; + + SparseCooTests test(memspace); + + result += test.constructor(50, 50, 2); + result += test.constructor(50, 100, 2); + + result += test.setDataPointers(50, 50, 2); + result += test.setDataPointers(50, 100, 2); + + result += test.getCsrRowData(); +} + +int main(int, char**) +{ + TestingResults result; + runTests("CPU", memory::HOST, result); + runTests("CPU", memory::HOST, result); + runTests("CPU", memory::HOST, result); + + return result.summary(); +} diff --git a/tests/UnitTests/LinearAlgebra/SparseMatrix/runSparseTests.cpp b/tests/UnitTests/LinearAlgebra/SparseMatrix/runSparseTests.cpp index 44ec49195..393e9ad0b 100644 --- a/tests/UnitTests/LinearAlgebra/SparseMatrix/runSparseTests.cpp +++ b/tests/UnitTests/LinearAlgebra/SparseMatrix/runSparseTests.cpp @@ -43,6 +43,7 @@ int main(int, char**) TestingResults result; runTests("CPU", memory::HOST, result); runTests("CPU", memory::HOST, result); + runTests("CPU", memory::HOST, result); return result.summary(); }