From e27ec57c7280a80cf7d0e1f3a1cdf399040dd51a Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Sat, 22 Jun 2024 21:29:45 -0400 Subject: [PATCH 1/5] Make coo2csr a standalone function. --- resolve/matrix/CMakeLists.txt | 2 + resolve/matrix/MatrixHandler.cpp | 283 ++++++++++++++++--------------- resolve/matrix/utilities.cpp | 163 ++++++++++++++++++ resolve/matrix/utilities.hpp | 17 ++ 4 files changed, 325 insertions(+), 140 deletions(-) create mode 100644 resolve/matrix/utilities.cpp create mode 100644 resolve/matrix/utilities.hpp diff --git a/resolve/matrix/CMakeLists.txt b/resolve/matrix/CMakeLists.txt index ad38dd5d..d6159742 100644 --- a/resolve/matrix/CMakeLists.txt +++ b/resolve/matrix/CMakeLists.txt @@ -15,6 +15,7 @@ set(Matrix_SRC Coo.cpp MatrixHandler.cpp MatrixHandlerCpu.cpp + utilities.h ) # C++ code that depends on CUDA SDK libraries @@ -35,6 +36,7 @@ set(Matrix_HEADER_INSTALL Csr.hpp Csc.hpp MatrixHandler.hpp + utilities.cpp ) # Build shared library ReSolve::matrix diff --git a/resolve/matrix/MatrixHandler.cpp b/resolve/matrix/MatrixHandler.cpp index 3fc9a492..2a6b46fb 100644 --- a/resolve/matrix/MatrixHandler.cpp +++ b/resolve/matrix/MatrixHandler.cpp @@ -116,148 +116,151 @@ namespace ReSolve { */ int MatrixHandler::coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, memory::MemorySpace memspace) { - //count nnzs first - index_type nnz_unpacked = 0; - index_type nnz = A_coo->getNnz(); - index_type n = A_coo->getNumRows(); - bool symmetric = A_coo->symmetric(); - bool expanded = A_coo->expanded(); - - index_type* nnz_counts = new index_type[n]; - std::fill_n(nnz_counts, n, 0); - index_type* coo_rows = A_coo->getRowData(memory::HOST); - index_type* coo_cols = A_coo->getColData(memory::HOST); - real_type* coo_vals = A_coo->getValues( memory::HOST); - - index_type* diag_control = new index_type[n]; //for DEDUPLICATION of the diagonal - std::fill_n(diag_control, n, 0); - index_type nnz_unpacked_no_duplicates = 0; - index_type nnz_no_duplicates = nnz; - - - //maybe check if they exist? - for (index_type i = 0; i < nnz; ++i) - { - nnz_counts[coo_rows[i]]++; - nnz_unpacked++; - nnz_unpacked_no_duplicates++; - if ((coo_rows[i] != coo_cols[i])&& (symmetric) && (!expanded)) - { - nnz_counts[coo_cols[i]]++; - nnz_unpacked++; - nnz_unpacked_no_duplicates++; - } - if (coo_rows[i] == coo_cols[i]){ - if (diag_control[coo_rows[i]] > 0) { - //duplicate - nnz_unpacked_no_duplicates--; - nnz_no_duplicates--; - } - diag_control[coo_rows[i]]++; - } - } - A_csr->setExpanded(true); - A_csr->setNnzExpanded(nnz_unpacked_no_duplicates); - index_type* csr_ia = new index_type[n+1]; - std::fill_n(csr_ia, n + 1, 0); - index_type* csr_ja = new index_type[nnz_unpacked]; - real_type* csr_a = new real_type[nnz_unpacked]; - index_type* nnz_shifts = new index_type[n]; - std::fill_n(nnz_shifts, n , 0); - - IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; - - csr_ia[0] = 0; - - for (index_type i = 1; i < n + 1; ++i){ - csr_ia[i] = csr_ia[i - 1] + nnz_counts[i - 1] - (diag_control[i-1] - 1); - } - - int r, start; - - - for (index_type i = 0; i < nnz; ++i){ - //which row - r = coo_rows[i]; - start = csr_ia[r]; - - if ((start + nnz_shifts[r]) > nnz_unpacked) { - out::warning() << "index out of bounds (case 1) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; - } - if ((r == coo_cols[i]) && (diag_control[r] > 1)) {//diagonal, and there are duplicates - bool already_there = false; - for (index_type j = start; j < start + nnz_shifts[r]; ++j) - { - index_type c = tmp[j].getIdx(); - if (c == r) { - real_type val = tmp[j].getValue(); - val += coo_vals[i]; - tmp[j].setValue(val); - already_there = true; - out::warning() << " duplicate found, row " << c << " adding in place " << j << " current value: " << val << std::endl; - } - } - if (!already_there){ // first time this duplicates appears - - tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); - tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); - - nnz_shifts[r]++; - } - } else {//not diagonal - tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); - tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); - nnz_shifts[r]++; - - if ((coo_rows[i] != coo_cols[i]) && (symmetric == 1)) - { - r = coo_cols[i]; - start = csr_ia[r]; - - if ((start + nnz_shifts[r]) > nnz_unpacked) - out::warning() << "index out of bounds (case 2) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; - tmp[start + nnz_shifts[r]].setIdx(coo_rows[i]); - tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); - nnz_shifts[r]++; - } - } - } - //now sort whatever is inside rows - - for (int i = 0; i < n; ++i) - { - - //now sorting (and adding 1) - int colStart = csr_ia[i]; - int colEnd = csr_ia[i + 1]; - int length = colEnd - colStart; - std::sort(&tmp[colStart],&tmp[colStart] + length); - } - - for (index_type i = 0; i < nnz_unpacked; ++i) - { - csr_ja[i] = tmp[i].getIdx(); - csr_a[i] = tmp[i].getValue(); - } -#if 0 - for (int i = 0; isetNnz(nnz_no_duplicates); + matrix::coo2csr(A_coo, A_csr); A_csr->updateData(csr_ia, csr_ja, csr_a, memory::HOST, memspace); - delete [] nnz_counts; - delete [] tmp; - delete [] nnz_shifts; - delete [] csr_ia; - delete [] csr_ja; - delete [] csr_a; - delete [] diag_control; +// //count nnzs first +// index_type nnz_unpacked = 0; +// index_type nnz = A_coo->getNnz(); +// index_type n = A_coo->getNumRows(); +// bool symmetric = A_coo->symmetric(); +// bool expanded = A_coo->expanded(); + +// index_type* nnz_counts = new index_type[n]; +// std::fill_n(nnz_counts, n, 0); +// index_type* coo_rows = A_coo->getRowData(memory::HOST); +// index_type* coo_cols = A_coo->getColData(memory::HOST); +// real_type* coo_vals = A_coo->getValues( memory::HOST); + +// index_type* diag_control = new index_type[n]; //for DEDUPLICATION of the diagonal +// std::fill_n(diag_control, n, 0); +// index_type nnz_unpacked_no_duplicates = 0; +// index_type nnz_no_duplicates = nnz; + + +// //maybe check if they exist? +// for (index_type i = 0; i < nnz; ++i) +// { +// nnz_counts[coo_rows[i]]++; +// nnz_unpacked++; +// nnz_unpacked_no_duplicates++; +// if ((coo_rows[i] != coo_cols[i])&& (symmetric) && (!expanded)) +// { +// nnz_counts[coo_cols[i]]++; +// nnz_unpacked++; +// nnz_unpacked_no_duplicates++; +// } +// if (coo_rows[i] == coo_cols[i]){ +// if (diag_control[coo_rows[i]] > 0) { +// //duplicate +// nnz_unpacked_no_duplicates--; +// nnz_no_duplicates--; +// } +// diag_control[coo_rows[i]]++; +// } +// } +// A_csr->setExpanded(true); +// A_csr->setNnzExpanded(nnz_unpacked_no_duplicates); +// index_type* csr_ia = new index_type[n+1]; +// std::fill_n(csr_ia, n + 1, 0); +// index_type* csr_ja = new index_type[nnz_unpacked]; +// real_type* csr_a = new real_type[nnz_unpacked]; +// index_type* nnz_shifts = new index_type[n]; +// std::fill_n(nnz_shifts, n , 0); + +// IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; + +// csr_ia[0] = 0; + +// for (index_type i = 1; i < n + 1; ++i){ +// csr_ia[i] = csr_ia[i - 1] + nnz_counts[i - 1] - (diag_control[i-1] - 1); +// } + +// int r, start; + + +// for (index_type i = 0; i < nnz; ++i){ +// //which row +// r = coo_rows[i]; +// start = csr_ia[r]; + +// if ((start + nnz_shifts[r]) > nnz_unpacked) { +// out::warning() << "index out of bounds (case 1) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; +// } +// if ((r == coo_cols[i]) && (diag_control[r] > 1)) {//diagonal, and there are duplicates +// bool already_there = false; +// for (index_type j = start; j < start + nnz_shifts[r]; ++j) +// { +// index_type c = tmp[j].getIdx(); +// if (c == r) { +// real_type val = tmp[j].getValue(); +// val += coo_vals[i]; +// tmp[j].setValue(val); +// already_there = true; +// out::warning() << " duplicate found, row " << c << " adding in place " << j << " current value: " << val << std::endl; +// } +// } +// if (!already_there){ // first time this duplicates appears + +// tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); +// tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); + +// nnz_shifts[r]++; +// } +// } else {//not diagonal +// tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); +// tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); +// nnz_shifts[r]++; + +// if ((coo_rows[i] != coo_cols[i]) && (symmetric == 1)) +// { +// r = coo_cols[i]; +// start = csr_ia[r]; + +// if ((start + nnz_shifts[r]) > nnz_unpacked) +// out::warning() << "index out of bounds (case 2) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; +// tmp[start + nnz_shifts[r]].setIdx(coo_rows[i]); +// tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); +// nnz_shifts[r]++; +// } +// } +// } +// //now sort whatever is inside rows + +// for (int i = 0; i < n; ++i) +// { + +// //now sorting (and adding 1) +// int colStart = csr_ia[i]; +// int colEnd = csr_ia[i + 1]; +// int length = colEnd - colStart; +// std::sort(&tmp[colStart],&tmp[colStart] + length); +// } + +// for (index_type i = 0; i < nnz_unpacked; ++i) +// { +// csr_ja[i] = tmp[i].getIdx(); +// csr_a[i] = tmp[i].getValue(); +// } +// #if 0 +// for (int i = 0; isetNnz(nnz_no_duplicates); + // A_csr->updateData(csr_ia, csr_ja, csr_a, memory::HOST, memspace); + + // delete [] nnz_counts; + // delete [] tmp; + // delete [] nnz_shifts; + // delete [] csr_ia; + // delete [] csr_ja; + // delete [] csr_a; + // delete [] diag_control; return 0; } diff --git a/resolve/matrix/utilities.cpp b/resolve/matrix/utilities.cpp new file mode 100644 index 00000000..d23850bf --- /dev/null +++ b/resolve/matrix/utilities.cpp @@ -0,0 +1,163 @@ +#include +#include +#include +#include "utilities.hpp" + +namespace ReSolve +{ + namespace matrix + { + /** + * @brief Creates a CSR from a COO matrix. + * + * @param[in] A_coo + * @param[out] A_csr + * @return int - Error code, 0 if successful. + * + * @pre `A_coo` is a valid sparse matrix in unorderes COO format. + * Duplicates are allowed. Up-to-date values and indices must be + * on host. + * + * @post `A_csr` is representing the same matrix as `A_coo` but in CSR + * format. `A_csr` is allocated and stored on host. + * + * @invariant `A_coo` is not changed. + */ + int coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr) + { + //count nnzs first + index_type nnz_unpacked = 0; + index_type nnz = A_coo->getNnz(); + index_type n = A_coo->getNumRows(); + bool symmetric = A_coo->symmetric(); + bool expanded = A_coo->expanded(); + + index_type* nnz_counts = new index_type[n]; + std::fill_n(nnz_counts, n, 0); + index_type* coo_rows = A_coo->getRowData(memory::HOST); + index_type* coo_cols = A_coo->getColData(memory::HOST); + real_type* coo_vals = A_coo->getValues( memory::HOST); + + index_type* diag_control = new index_type[n]; // for DEDUPLICATION of the diagonal + std::fill_n(diag_control, n, 0); + index_type nnz_unpacked_no_duplicates = 0; + index_type nnz_no_duplicates = nnz; + + // maybe check if they exist? + for (index_type i = 0; i < nnz; ++i) { + nnz_counts[coo_rows[i]]++; + nnz_unpacked++; + nnz_unpacked_no_duplicates++; + if ((coo_rows[i] != coo_cols[i]) && (symmetric) && (!expanded)) + { + nnz_counts[coo_cols[i]]++; + nnz_unpacked++; + nnz_unpacked_no_duplicates++; + } + if (coo_rows[i] == coo_cols[i]) { + if (diag_control[coo_rows[i]] > 0) { + // duplicate + nnz_unpacked_no_duplicates--; + nnz_no_duplicates--; + } + diag_control[coo_rows[i]]++; + } + } + A_csr->setExpanded(true); + A_csr->setNnzExpanded(nnz_unpacked_no_duplicates); + index_type* csr_ia = new index_type[n+1]; + std::fill_n(csr_ia, n + 1, 0); + index_type* csr_ja = new index_type[nnz_unpacked]; + real_type* csr_a = new real_type[nnz_unpacked]; + index_type* nnz_shifts = new index_type[n]; + std::fill_n(nnz_shifts, n , 0); + + IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; + + csr_ia[0] = 0; + + for (index_type i = 1; i < n + 1; ++i) { + csr_ia[i] = csr_ia[i - 1] + nnz_counts[i - 1] - (diag_control[i-1] - 1); + } + + int r, start; + + + for (index_type i = 0; i < nnz; ++i){ + //which row + r = coo_rows[i]; + start = csr_ia[r]; + + if ((start + nnz_shifts[r]) > nnz_unpacked) { + out::warning() << "index out of bounds (case 1) start: " << start + << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; + } + if ((r == coo_cols[i]) && (diag_control[r] > 1)) {//diagonal, and there are duplicates + bool already_there = false; + for (index_type j = start; j < start + nnz_shifts[r]; ++j) + { + index_type c = tmp[j].getIdx(); + if (c == r) { + real_type val = tmp[j].getValue(); + val += coo_vals[i]; + tmp[j].setValue(val); + already_there = true; + out::warning() << " duplicate found, row " << c << " adding in place " << j << " current value: " << val << std::endl; + } + } + if (!already_there) { // first time this duplicates appears + + tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); + tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); + + nnz_shifts[r]++; + } + } else { // non-diagonal + tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); + tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); + nnz_shifts[r]++; + + if ((coo_rows[i] != coo_cols[i]) && (symmetric == 1)) + { + r = coo_cols[i]; + start = csr_ia[r]; + + if ((start + nnz_shifts[r]) > nnz_unpacked) + out::warning() << "index out of bounds (case 2) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; + tmp[start + nnz_shifts[r]].setIdx(coo_rows[i]); + tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); + nnz_shifts[r]++; + } + } + } + + //now sort whatever is inside rows + for (int i = 0; i < n; ++i) { + //now sorting (and adding 1) + int colStart = csr_ia[i]; + int colEnd = csr_ia[i + 1]; + int length = colEnd - colStart; + std::sort(&tmp[colStart], &tmp[colStart] + length); + } + + for (index_type i = 0; i < nnz_unpacked; ++i) + { + csr_ja[i] = tmp[i].getIdx(); + csr_a[i] = tmp[i].getValue(); + } + A_csr->setNnz(nnz_no_duplicates); + A_csr->updateData(csr_ia, csr_ja, csr_a, memory::HOST, memspace); + + delete [] nnz_counts; + delete [] tmp; + delete [] nnz_shifts; + delete [] csr_ia; + delete [] csr_ja; + delete [] csr_a; + delete [] diag_control; + + return 0; + + } + } +} \ No newline at end of file diff --git a/resolve/matrix/utilities.hpp b/resolve/matrix/utilities.hpp new file mode 100644 index 00000000..20582eaa --- /dev/null +++ b/resolve/matrix/utilities.hpp @@ -0,0 +1,17 @@ +#pragma once + +namespace ReSolve +{ + namespace matrix + { + // Forward declarations + class Coo; + class Csr; + + /// @brief + /// @param A_coo + /// @param A_csr + /// @return + int coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr); + } +} From d49f5611fd0bfe3fd1ee58354726a3485243f31f Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Mon, 24 Jun 2024 15:43:44 -0400 Subject: [PATCH 2/5] Remove duplicated coo2csr code. --- resolve/matrix/CMakeLists.txt | 4 +- resolve/matrix/Csr.cpp | 144 +---------------------------- resolve/matrix/Csr.hpp | 2 - resolve/matrix/MatrixHandler.cpp | 150 +------------------------------ resolve/matrix/utilities.cpp | 53 ++++++----- resolve/matrix/utilities.hpp | 9 +- 6 files changed, 41 insertions(+), 321 deletions(-) diff --git a/resolve/matrix/CMakeLists.txt b/resolve/matrix/CMakeLists.txt index d6159742..ba2935ab 100644 --- a/resolve/matrix/CMakeLists.txt +++ b/resolve/matrix/CMakeLists.txt @@ -15,7 +15,7 @@ set(Matrix_SRC Coo.cpp MatrixHandler.cpp MatrixHandlerCpu.cpp - utilities.h + utilities.cpp ) # C++ code that depends on CUDA SDK libraries @@ -36,7 +36,7 @@ set(Matrix_HEADER_INSTALL Csr.hpp Csc.hpp MatrixHandler.hpp - utilities.cpp + utilities.hpp ) # Build shared library ReSolve::matrix diff --git a/resolve/matrix/Csr.cpp b/resolve/matrix/Csr.cpp index 325bb76b..308f1a29 100644 --- a/resolve/matrix/Csr.cpp +++ b/resolve/matrix/Csr.cpp @@ -7,6 +7,7 @@ #include "Coo.hpp" #include #include +#include namespace ReSolve { @@ -35,7 +36,7 @@ namespace ReSolve A_coo->symmetric(), A_coo->expanded()) { - coo2csr(A_coo, memspace); + matrix::coo2csr(A_coo, this, memspace); } /** @@ -374,149 +375,10 @@ namespace ReSolve assert(nnz_ == A_coo->getNnz()); assert(is_symmetric_ == A_coo->symmetric()); // <- Do we need to check for this? - return coo2csr(A_coo, memspaceOut); + return matrix::coo2csr(A_coo, this, memspaceOut); } - int matrix::Csr::coo2csr(matrix::Coo* A_coo, memory::MemorySpace memspace) - { - //count nnzs first - index_type nnz_unpacked = 0; - index_type nnz = A_coo->getNnz(); - index_type n = A_coo->getNumRows(); - bool symmetric = A_coo->symmetric(); - bool expanded = A_coo->expanded(); - - index_type* nnz_counts = new index_type[n]; - std::fill_n(nnz_counts, n, 0); - index_type* coo_rows = A_coo->getRowData(memory::HOST); - index_type* coo_cols = A_coo->getColData(memory::HOST); - real_type* coo_vals = A_coo->getValues( memory::HOST); - - index_type* diag_control = new index_type[n]; //for DEDUPLICATION of the diagonal - std::fill_n(diag_control, n, 0); - index_type nnz_unpacked_no_duplicates = 0; - index_type nnz_no_duplicates = nnz; - - - //maybe check if they exist? - for (index_type i = 0; i < nnz; ++i) - { - nnz_counts[coo_rows[i]]++; - nnz_unpacked++; - nnz_unpacked_no_duplicates++; - if ((coo_rows[i] != coo_cols[i])&& (symmetric) && (!expanded)) - { - nnz_counts[coo_cols[i]]++; - nnz_unpacked++; - nnz_unpacked_no_duplicates++; - } - if (coo_rows[i] == coo_cols[i]){ - if (diag_control[coo_rows[i]] > 0) { - //duplicate - nnz_unpacked_no_duplicates--; - nnz_no_duplicates--; - } - diag_control[coo_rows[i]]++; - } - } - this->setExpanded(true); - this->setNnzExpanded(nnz_unpacked_no_duplicates); - index_type* csr_ia = new index_type[n+1]; - std::fill_n(csr_ia, n + 1, 0); - index_type* csr_ja = new index_type[nnz_unpacked]; - real_type* csr_a = new real_type[nnz_unpacked]; - index_type* nnz_shifts = new index_type[n]; - std::fill_n(nnz_shifts, n , 0); - - IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; - - csr_ia[0] = 0; - - for (index_type i = 1; i < n + 1; ++i){ - csr_ia[i] = csr_ia[i - 1] + nnz_counts[i - 1] - (diag_control[i-1] - 1); - } - - int r, start; - - - for (index_type i = 0; i < nnz; ++i){ - //which row - r = coo_rows[i]; - start = csr_ia[r]; - - if ((start + nnz_shifts[r]) > nnz_unpacked) { - out::warning() << "index out of bounds (case 1) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; - } - if ((r == coo_cols[i]) && (diag_control[r] > 1)) {//diagonal, and there are duplicates - bool already_there = false; - for (index_type j = start; j < start + nnz_shifts[r]; ++j) - { - index_type c = tmp[j].getIdx(); - if (c == r) { - real_type val = tmp[j].getValue(); - val += coo_vals[i]; - tmp[j].setValue(val); - already_there = true; - out::warning() << " duplicate found, row " << c << " adding in place " << j << " current value: " << val << std::endl; - } - } - if (!already_there){ // first time this duplicates appears - - tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); - tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); - - nnz_shifts[r]++; - } - } else {//not diagonal - tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); - tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); - nnz_shifts[r]++; - - if ((coo_rows[i] != coo_cols[i]) && (symmetric == 1)) - { - r = coo_cols[i]; - start = csr_ia[r]; - - if ((start + nnz_shifts[r]) > nnz_unpacked) - out::warning() << "index out of bounds (case 2) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; - tmp[start + nnz_shifts[r]].setIdx(coo_rows[i]); - tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); - nnz_shifts[r]++; - } - } - } - //now sort whatever is inside rows - - for (int i = 0; i < n; ++i) - { - //now sorting (and adding 1) - int colStart = csr_ia[i]; - int colEnd = csr_ia[i + 1]; - int length = colEnd - colStart; - std::sort(&tmp[colStart],&tmp[colStart] + length); - } - - for (index_type i = 0; i < nnz_unpacked; ++i) - { - csr_ja[i] = tmp[i].getIdx(); - csr_a[i] = tmp[i].getValue(); - } - - this->setNnz(nnz_no_duplicates); - this->updateData(csr_ia, csr_ja, csr_a, memory::HOST, memspace); - - delete [] nnz_counts; - delete [] tmp; - delete [] nnz_shifts; - delete [] csr_ia; - delete [] csr_ja; - delete [] csr_a; - delete [] diag_control; - - return 0; - } - /** * @brief Prints matrix data. * diff --git a/resolve/matrix/Csr.hpp b/resolve/matrix/Csr.hpp index 7841a7fb..a3896a1d 100644 --- a/resolve/matrix/Csr.hpp +++ b/resolve/matrix/Csr.hpp @@ -49,8 +49,6 @@ namespace ReSolve { namespace matrix { int updateFromCoo(matrix::Coo* mat, memory::MemorySpace memspaceOut); - private: - int coo2csr(matrix::Coo* mat, memory::MemorySpace memspace); }; }} // namespace ReSolve::matrix diff --git a/resolve/matrix/MatrixHandler.cpp b/resolve/matrix/MatrixHandler.cpp index 2a6b46fb..127f03b2 100644 --- a/resolve/matrix/MatrixHandler.cpp +++ b/resolve/matrix/MatrixHandler.cpp @@ -6,7 +6,7 @@ #include #include #include -#include +#include #include "MatrixHandler.hpp" #include "MatrixHandlerCpu.hpp" @@ -116,153 +116,7 @@ namespace ReSolve { */ int MatrixHandler::coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, memory::MemorySpace memspace) { - matrix::coo2csr(A_coo, A_csr); - A_csr->updateData(csr_ia, csr_ja, csr_a, memory::HOST, memspace); - -// //count nnzs first -// index_type nnz_unpacked = 0; -// index_type nnz = A_coo->getNnz(); -// index_type n = A_coo->getNumRows(); -// bool symmetric = A_coo->symmetric(); -// bool expanded = A_coo->expanded(); - -// index_type* nnz_counts = new index_type[n]; -// std::fill_n(nnz_counts, n, 0); -// index_type* coo_rows = A_coo->getRowData(memory::HOST); -// index_type* coo_cols = A_coo->getColData(memory::HOST); -// real_type* coo_vals = A_coo->getValues( memory::HOST); - -// index_type* diag_control = new index_type[n]; //for DEDUPLICATION of the diagonal -// std::fill_n(diag_control, n, 0); -// index_type nnz_unpacked_no_duplicates = 0; -// index_type nnz_no_duplicates = nnz; - - -// //maybe check if they exist? -// for (index_type i = 0; i < nnz; ++i) -// { -// nnz_counts[coo_rows[i]]++; -// nnz_unpacked++; -// nnz_unpacked_no_duplicates++; -// if ((coo_rows[i] != coo_cols[i])&& (symmetric) && (!expanded)) -// { -// nnz_counts[coo_cols[i]]++; -// nnz_unpacked++; -// nnz_unpacked_no_duplicates++; -// } -// if (coo_rows[i] == coo_cols[i]){ -// if (diag_control[coo_rows[i]] > 0) { -// //duplicate -// nnz_unpacked_no_duplicates--; -// nnz_no_duplicates--; -// } -// diag_control[coo_rows[i]]++; -// } -// } -// A_csr->setExpanded(true); -// A_csr->setNnzExpanded(nnz_unpacked_no_duplicates); -// index_type* csr_ia = new index_type[n+1]; -// std::fill_n(csr_ia, n + 1, 0); -// index_type* csr_ja = new index_type[nnz_unpacked]; -// real_type* csr_a = new real_type[nnz_unpacked]; -// index_type* nnz_shifts = new index_type[n]; -// std::fill_n(nnz_shifts, n , 0); - -// IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; - -// csr_ia[0] = 0; - -// for (index_type i = 1; i < n + 1; ++i){ -// csr_ia[i] = csr_ia[i - 1] + nnz_counts[i - 1] - (diag_control[i-1] - 1); -// } - -// int r, start; - - -// for (index_type i = 0; i < nnz; ++i){ -// //which row -// r = coo_rows[i]; -// start = csr_ia[r]; - -// if ((start + nnz_shifts[r]) > nnz_unpacked) { -// out::warning() << "index out of bounds (case 1) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; -// } -// if ((r == coo_cols[i]) && (diag_control[r] > 1)) {//diagonal, and there are duplicates -// bool already_there = false; -// for (index_type j = start; j < start + nnz_shifts[r]; ++j) -// { -// index_type c = tmp[j].getIdx(); -// if (c == r) { -// real_type val = tmp[j].getValue(); -// val += coo_vals[i]; -// tmp[j].setValue(val); -// already_there = true; -// out::warning() << " duplicate found, row " << c << " adding in place " << j << " current value: " << val << std::endl; -// } -// } -// if (!already_there){ // first time this duplicates appears - -// tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); -// tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); - -// nnz_shifts[r]++; -// } -// } else {//not diagonal -// tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]); -// tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); -// nnz_shifts[r]++; - -// if ((coo_rows[i] != coo_cols[i]) && (symmetric == 1)) -// { -// r = coo_cols[i]; -// start = csr_ia[r]; - -// if ((start + nnz_shifts[r]) > nnz_unpacked) -// out::warning() << "index out of bounds (case 2) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; -// tmp[start + nnz_shifts[r]].setIdx(coo_rows[i]); -// tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); -// nnz_shifts[r]++; -// } -// } -// } -// //now sort whatever is inside rows - -// for (int i = 0; i < n; ++i) -// { - -// //now sorting (and adding 1) -// int colStart = csr_ia[i]; -// int colEnd = csr_ia[i + 1]; -// int length = colEnd - colStart; -// std::sort(&tmp[colStart],&tmp[colStart] + length); -// } - -// for (index_type i = 0; i < nnz_unpacked; ++i) -// { -// csr_ja[i] = tmp[i].getIdx(); -// csr_a[i] = tmp[i].getValue(); -// } -// #if 0 -// for (int i = 0; isetNnz(nnz_no_duplicates); - // A_csr->updateData(csr_ia, csr_ja, csr_a, memory::HOST, memspace); - - // delete [] nnz_counts; - // delete [] tmp; - // delete [] nnz_shifts; - // delete [] csr_ia; - // delete [] csr_ja; - // delete [] csr_a; - // delete [] diag_control; - - return 0; + return matrix::coo2csr(A_coo, A_csr, memspace); } /** diff --git a/resolve/matrix/utilities.cpp b/resolve/matrix/utilities.cpp index d23850bf..a7e1a3e8 100644 --- a/resolve/matrix/utilities.cpp +++ b/resolve/matrix/utilities.cpp @@ -1,10 +1,13 @@ #include #include #include +#include +#include #include "utilities.hpp" namespace ReSolve { + using out = io::Logger; namespace matrix { /** @@ -14,46 +17,47 @@ namespace ReSolve * @param[out] A_csr * @return int - Error code, 0 if successful. * - * @pre `A_coo` is a valid sparse matrix in unorderes COO format. + * @pre `A_coo` is a valid sparse matrix in unordered COO format. * Duplicates are allowed. Up-to-date values and indices must be * on host. * - * @post `A_csr` is representing the same matrix as `A_coo` but in CSR - * format. `A_csr` is allocated and stored on host. + * @post `A_csr` is representing the same matrix as `A_coo` but in + * _general_ CSR format. `A_csr` is allocated and stored on host. * * @invariant `A_coo` is not changed. */ - int coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr) + int coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, memory::MemorySpace memspace) { //count nnzs first index_type nnz_unpacked = 0; - index_type nnz = A_coo->getNnz(); - index_type n = A_coo->getNumRows(); - bool symmetric = A_coo->symmetric(); - bool expanded = A_coo->expanded(); + const index_type nnz = A_coo->getNnz(); + const index_type n = A_coo->getNumRows(); + const bool symmetric = A_coo->symmetric(); + const bool expanded = A_coo->expanded(); - index_type* nnz_counts = new index_type[n]; + index_type* nnz_counts = new index_type[n]; ///< Number of elements per row std::fill_n(nnz_counts, n, 0); - index_type* coo_rows = A_coo->getRowData(memory::HOST); - index_type* coo_cols = A_coo->getColData(memory::HOST); - real_type* coo_vals = A_coo->getValues( memory::HOST); + const index_type* coo_rows = A_coo->getRowData(memory::HOST); + const index_type* coo_cols = A_coo->getColData(memory::HOST); + const real_type* coo_vals = A_coo->getValues( memory::HOST); index_type* diag_control = new index_type[n]; // for DEDUPLICATION of the diagonal std::fill_n(diag_control, n, 0); index_type nnz_unpacked_no_duplicates = 0; index_type nnz_no_duplicates = nnz; - // maybe check if they exist? + // Count matrix elements for (index_type i = 0; i < nnz; ++i) { nnz_counts[coo_rows[i]]++; nnz_unpacked++; nnz_unpacked_no_duplicates++; - if ((coo_rows[i] != coo_cols[i]) && (symmetric) && (!expanded)) - { + // Count elements added after expansion + if ((coo_rows[i] != coo_cols[i]) && (symmetric) && (!expanded)) { nnz_counts[coo_cols[i]]++; nnz_unpacked++; nnz_unpacked_no_duplicates++; } + // Bookkeeping of diagonal elements that were counted if (coo_rows[i] == coo_cols[i]) { if (diag_control[coo_rows[i]] > 0) { // duplicate @@ -65,6 +69,8 @@ namespace ReSolve } A_csr->setExpanded(true); A_csr->setNnzExpanded(nnz_unpacked_no_duplicates); + + // Allocate matrix format conversion workspace index_type* csr_ia = new index_type[n+1]; std::fill_n(csr_ia, n + 1, 0); index_type* csr_ja = new index_type[nnz_unpacked]; @@ -74,16 +80,15 @@ namespace ReSolve IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; + // Set CSR row pointers csr_ia[0] = 0; - for (index_type i = 1; i < n + 1; ++i) { csr_ia[i] = csr_ia[i - 1] + nnz_counts[i - 1] - (diag_control[i-1] - 1); } int r, start; - - for (index_type i = 0; i < nnz; ++i){ + for (index_type i = 0; i < nnz; ++i) { //which row r = coo_rows[i]; start = csr_ia[r]; @@ -102,7 +107,8 @@ namespace ReSolve val += coo_vals[i]; tmp[j].setValue(val); already_there = true; - out::warning() << " duplicate found, row " << c << " adding in place " << j << " current value: " << val << std::endl; + out::warning() << " duplicate found, row " << c + << " adding in place " << j << " current value: " << val << std::endl; } } if (!already_there) { // first time this duplicates appears @@ -117,13 +123,14 @@ namespace ReSolve tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); nnz_shifts[r]++; - if ((coo_rows[i] != coo_cols[i]) && (symmetric == 1)) - { + if ((coo_rows[i] != coo_cols[i]) && (symmetric == 1)) { r = coo_cols[i]; start = csr_ia[r]; - if ((start + nnz_shifts[r]) > nnz_unpacked) - out::warning() << "index out of bounds (case 2) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; + if ((start + nnz_shifts[r]) > nnz_unpacked) { + out::warning() << "index out of bounds (case 2) start: " << start + << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl; + } tmp[start + nnz_shifts[r]].setIdx(coo_rows[i]); tmp[start + nnz_shifts[r]].setValue(coo_vals[i]); nnz_shifts[r]++; diff --git a/resolve/matrix/utilities.hpp b/resolve/matrix/utilities.hpp index 20582eaa..ecaee79e 100644 --- a/resolve/matrix/utilities.hpp +++ b/resolve/matrix/utilities.hpp @@ -1,5 +1,7 @@ #pragma once +#include + namespace ReSolve { namespace matrix @@ -8,10 +10,7 @@ namespace ReSolve class Coo; class Csr; - /// @brief - /// @param A_coo - /// @param A_csr - /// @return - int coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr); + /// @brief Converts symmetric or general COO to general CSR matrix + int coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, memory::MemorySpace memspace); } } From 52fa85db2b49e2b3d09691591a62b6786099b74d Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Tue, 25 Jun 2024 09:08:29 -0400 Subject: [PATCH 3/5] Fix missing header for sort function. --- resolve/matrix/utilities.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/resolve/matrix/utilities.cpp b/resolve/matrix/utilities.cpp index a7e1a3e8..7ad866f5 100644 --- a/resolve/matrix/utilities.cpp +++ b/resolve/matrix/utilities.cpp @@ -1,3 +1,4 @@ +#include #include #include #include From f430b29e43179e123ab0664b6219359abb6cbec2 Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Tue, 2 Jul 2024 18:25:47 -0400 Subject: [PATCH 4/5] Rename matrix/utilities.* files --- resolve/matrix/{utilities.cpp => Utilities.cpp} | 0 resolve/matrix/{utilities.hpp => Utilities.hpp} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename resolve/matrix/{utilities.cpp => Utilities.cpp} (100%) rename resolve/matrix/{utilities.hpp => Utilities.hpp} (100%) diff --git a/resolve/matrix/utilities.cpp b/resolve/matrix/Utilities.cpp similarity index 100% rename from resolve/matrix/utilities.cpp rename to resolve/matrix/Utilities.cpp diff --git a/resolve/matrix/utilities.hpp b/resolve/matrix/Utilities.hpp similarity index 100% rename from resolve/matrix/utilities.hpp rename to resolve/matrix/Utilities.hpp From 31046e45294dd46dce23da46aba43702244a9c43 Mon Sep 17 00:00:00 2001 From: Slaven Peles Date: Tue, 2 Jul 2024 18:58:06 -0400 Subject: [PATCH 5/5] Fix errors caught by CI. --- resolve/matrix/CMakeLists.txt | 4 +-- resolve/matrix/Csr.cpp | 3 +- resolve/matrix/MatrixHandler.cpp | 2 +- resolve/matrix/Utilities.cpp | 39 ++++++++++++++++++-- resolve/utilities/misc/IndexValuePair.hpp | 43 ----------------------- 5 files changed, 41 insertions(+), 50 deletions(-) delete mode 100644 resolve/utilities/misc/IndexValuePair.hpp diff --git a/resolve/matrix/CMakeLists.txt b/resolve/matrix/CMakeLists.txt index ba2935ab..6a92975e 100644 --- a/resolve/matrix/CMakeLists.txt +++ b/resolve/matrix/CMakeLists.txt @@ -15,7 +15,7 @@ set(Matrix_SRC Coo.cpp MatrixHandler.cpp MatrixHandlerCpu.cpp - utilities.cpp + Utilities.cpp ) # C++ code that depends on CUDA SDK libraries @@ -36,7 +36,7 @@ set(Matrix_HEADER_INSTALL Csr.hpp Csc.hpp MatrixHandler.hpp - utilities.hpp + Utilities.hpp ) # Build shared library ReSolve::matrix diff --git a/resolve/matrix/Csr.cpp b/resolve/matrix/Csr.cpp index 308f1a29..268898bc 100644 --- a/resolve/matrix/Csr.cpp +++ b/resolve/matrix/Csr.cpp @@ -5,9 +5,8 @@ #include "Csr.hpp" #include "Coo.hpp" -#include #include -#include +#include namespace ReSolve { diff --git a/resolve/matrix/MatrixHandler.cpp b/resolve/matrix/MatrixHandler.cpp index 127f03b2..b61a93b5 100644 --- a/resolve/matrix/MatrixHandler.cpp +++ b/resolve/matrix/MatrixHandler.cpp @@ -6,7 +6,7 @@ #include #include #include -#include +#include #include "MatrixHandler.hpp" #include "MatrixHandlerCpu.hpp" diff --git a/resolve/matrix/Utilities.cpp b/resolve/matrix/Utilities.cpp index 7ad866f5..bd5a6f9b 100644 --- a/resolve/matrix/Utilities.cpp +++ b/resolve/matrix/Utilities.cpp @@ -2,12 +2,47 @@ #include #include #include -#include #include -#include "utilities.hpp" +#include "Utilities.hpp" namespace ReSolve { + /// @brief Helper class for COO matrix sorting + class IndexValuePair + { + public: + IndexValuePair() : idx_(0), value_(0.0) + {} + ~IndexValuePair() + {} + void setIdx (index_type new_idx) + { + idx_ = new_idx; + } + void setValue (real_type new_value) + { + value_ = new_value; + } + + index_type getIdx() + { + return idx_; + } + real_type getValue() + { + return value_; + } + + bool operator < (const IndexValuePair& str) const + { + return (idx_ < str.idx_); + } + + private: + index_type idx_; + real_type value_; + }; + using out = io::Logger; namespace matrix { diff --git a/resolve/utilities/misc/IndexValuePair.hpp b/resolve/utilities/misc/IndexValuePair.hpp deleted file mode 100644 index b09ccf97..00000000 --- a/resolve/utilities/misc/IndexValuePair.hpp +++ /dev/null @@ -1,43 +0,0 @@ -#pragma once - - -namespace ReSolve { - - /// @brief Helper class for COO matrix sorting - class IndexValuePair - { - public: - IndexValuePair() : idx_(0), value_(0.0) - {} - ~IndexValuePair() - {} - void setIdx (index_type new_idx) - { - idx_ = new_idx; - } - void setValue (real_type new_value) - { - value_ = new_value; - } - - index_type getIdx() - { - return idx_; - } - real_type getValue() - { - return value_; - } - - bool operator < (const IndexValuePair& str) const - { - return (idx_ < str.idx_); - } - - private: - index_type idx_; - real_type value_; - }; - -} // namespace ReSolve -