From 4af64c54d68738911b6d48531e583ce2dc17e4f6 Mon Sep 17 00:00:00 2001 From: superwhiskers Date: Fri, 5 Jul 2024 14:12:41 -0400 Subject: [PATCH] clean-slate coo2csr implementation + tests --- resolve/matrix/CMakeLists.txt | 4 +- resolve/matrix/Csr.cpp | 2 +- resolve/matrix/MatrixHandler.cpp | 2 +- resolve/matrix/Utilities.cpp | 225 ++++++++++++++++++ .../matrix/{utilities.hpp => Utilities.hpp} | 0 resolve/matrix/utilities.cpp | 171 ------------- tests/unit/matrix/CMakeLists.txt | 7 +- tests/unit/matrix/MatrixConversionTests.hpp | 213 +++++++++++++++++ .../unit/matrix/runMatrixConversionTests.cpp | 18 ++ 9 files changed, 466 insertions(+), 176 deletions(-) create mode 100644 resolve/matrix/Utilities.cpp rename resolve/matrix/{utilities.hpp => Utilities.hpp} (100%) delete mode 100644 resolve/matrix/utilities.cpp create mode 100644 tests/unit/matrix/MatrixConversionTests.hpp create mode 100644 tests/unit/matrix/runMatrixConversionTests.cpp 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..4735e68a 100644 --- a/resolve/matrix/Csr.cpp +++ b/resolve/matrix/Csr.cpp @@ -7,7 +7,7 @@ #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 new file mode 100644 index 00000000..50e9852c --- /dev/null +++ b/resolve/matrix/Utilities.cpp @@ -0,0 +1,225 @@ +#include +#include +#include + +#include "Utilities.hpp" + +#include +#include +#include +#include +#include + +namespace ReSolve +{ + using out = io::Logger; + namespace matrix + { + /** + * @brief Creates a CSR from a COO matrix. + * + * @param[in] A + * @param[out] B + * @return int - Error code, 0 if successful. + * + * @pre `A` is a valid sparse matrix in unordered COO format. Duplicates are allowed. + * Up-to-date values and indices must be on the host. + * + * @post `B` represents the same matrix as `A` but is in the CSR format. `B` is + * allocated and stored on the host. + * + * @invariant `A` is not changed. + */ + int coo2csr(matrix::Coo* A, matrix::Csr* B, memory::MemorySpace memspace) + { + // NOTE: a note on deduplication: + // currently, this function allocates more memory than necessary to contain + // the matrix. this is because it's impossible to cleanly check the amount + // of space each row needs ahead of time without storing a full dense matrix + // with the dimensions of A + // + // additionally, this necessitates that we shrink the rows to remove this + // unused space, which is costly and necessitates a lot of shifting + // + // one potential solution to this that i thought of was to store an array of + // bloom filters (maybe u64s or u128s) of the size of the number of columns, + // each filter indicating if there is a value stored at that column in a row + // + // if, during the preprocessing step, we encounter a potential duplicate, we + // backtrack to see if there actually is one. given that the number of + // duplicates is probably small, this shouldn't carry too much of a + // performance penalty + // + // if performance becomes a problem, this could be coupled with arrays + // maintaining the last and/or first seen indices in the coo arrays at which + // a row had an associated value, to speed up the backtracking process + // + // this would be applied during the first allocation phase when an element + // for a row is encountered, prior to the increment of its size in + // `new_rows` + + index_type* rows = A->getRowData(memory::HOST); + index_type* columns = A->getColData(memory::HOST); + real_type* values = A->getValues(memory::HOST); + + if (rows == nullptr || columns == nullptr || values == nullptr) { + return 0; + } + + index_type nnz = A->getNnz(); + index_type n_rows = A->getNumRows(); + index_type* new_rows = new index_type[n_rows + 1]; + std::fill_n(new_rows, n_rows + 1, 0); + + // NOTE: this is the only auxiliary storage buffer used by this conversion + // function. it is first used to track the number of values on the + // diagonal (if the matrix is symmetric and unexpanded), then it is + // used to track the amount of spaces used in each row's value and + // column data + std::unique_ptr used(new index_type[n_rows]); + std::fill_n(used.get(), n_rows, 0); + + // allocation stage + + // NOTE: so aside from tracking the number of mirrored values for allocation + // purposes, we also have to compute the amount of nonzeroes in each row + // since we're dealing with a COO input matrix. we store these values in + // row + 1 slots in new_rows, to avoid allocating extra space :) + + if (!A->symmetric() || A->expanded()) { + // NOTE: this branch is special cased because there is a bunch of extra + // bookkeeping done if a matrix is symmetric and unexpanded + for (index_type i = 0; i < nnz; i++) { + new_rows[rows[i] + 1]++; + } + } else { + bool upper_triangular = false; + for (index_type i = 0; i < nnz; i++) { + new_rows[rows[i] + 1]++; + if (rows[i] != columns[i]) { + used[columns[i]]++; + + if (rows[i] > columns[i] && upper_triangular) { + assert(false && "a matrix indicated to be symmetric triangular was not actually symmetric triangular"); + return -1; + } + upper_triangular = rows[i] < columns[i]; + } + } + } + + for (index_type row = 0; row < n_rows; row++) { + new_rows[row + 1] += new_rows[row] + used[row]; + used[row] = 0; + } + + index_type* new_columns = new index_type[new_rows[n_rows]]; + // TODO: we need to find a better way to fix this. what this fixes is + // that otherwise, the column values are all zero by default and + // this seems to mess with the way the insertion position is + // selected if the column of the pair we're looking to insert + // is zero + std::fill_n(new_columns, new_rows[n_rows], -1); + real_type* new_values = new real_type[new_rows[n_rows]]; + + // fill stage + + for (index_type i = 0; i < nnz; i++) { + index_type insertion_pos = + static_cast( + std::lower_bound(&new_columns[new_rows[rows[i]]], + &new_columns[new_rows[rows[i]] + used[rows[i]]], + columns[i]) - + new_columns); + + // NOTE: the stuff beyond here is basically insertion sort. i'm not + // sure if it's the best way of going about sorting the + // column-value pairs, but it seemed to me to be the most + // natural way of going about this. the only potential + // benefit to using something else (that i can forsee) would + // be that on really large matrices with many nonzeroes in a + // row, something like mergesort might be better or maybe + // a reimpl of plain std::sort. the advantage offered by + // insertion sort here is that we can do it while we fill + // the row, as opposed to doing sorting in a postprocessing + // step as was done prior + + if (new_columns[insertion_pos] == columns[i]) { + new_values[insertion_pos] = values[i]; + } else { + for (index_type offset = new_rows[rows[i]] + used[rows[i]]++; + offset > insertion_pos; + offset--) { + std::swap(new_columns[offset], new_columns[offset - 1]); + std::swap(new_values[offset], new_values[offset - 1]); + } + + new_columns[insertion_pos] = columns[i]; + new_values[insertion_pos] = values[i]; + } + + if ((A->symmetric() && !A->expanded()) && (columns[i] != rows[i])) { + index_type mirrored_insertion_pos = + static_cast( + std::lower_bound(&new_columns[new_rows[columns[i]]], + &new_columns[new_rows[columns[i]] + used[columns[i]]], + rows[i]) - + new_columns); + + if (new_columns[mirrored_insertion_pos] == rows[i]) { + new_values[mirrored_insertion_pos] = values[i]; + } else { + for (index_type offset = new_rows[columns[i]] + used[columns[i]]++; + offset > mirrored_insertion_pos; + offset--) { + std::swap(new_columns[offset], new_columns[offset - 1]); + std::swap(new_values[offset], new_values[offset - 1]); + } + + new_columns[mirrored_insertion_pos] = rows[i]; + new_values[mirrored_insertion_pos] = values[i]; + } + } + } + + // backshifting stage + + for (index_type row = 0; row < n_rows - 1; row++) { + index_type row_nnz = new_rows[row + 1] - new_rows[row]; + if (used[row] != row_nnz) { + index_type correction = row_nnz - used[row]; + + for (index_type corrected_row = row + 1; + corrected_row < n_rows; + corrected_row++) { + for (index_type offset = new_rows[corrected_row]; + offset < new_rows[corrected_row + 1]; + offset++) { + new_columns[offset - correction] = new_columns[offset]; + new_values[offset - correction] = new_values[offset]; + } + + new_rows[corrected_row] -= correction; + } + + new_rows[n_rows] -= correction; + } + } + + if (B->destroyMatrixData(memory::HOST) != 0 || + B->setMatrixData(new_rows, new_columns, new_values, memory::HOST) != 0) { + assert(false && "invalid state after coo -> csr conversion"); + return -1; + } + + // TODO: set data ownership / value ownership here. we'd be leaking + // memory here otherwise + + B->setSymmetric(A->symmetric()); + B->setNnz(new_rows[n_rows]); + B->setExpanded(true); + + return 0; + } + } // namespace matrix +} // namespace ReSolve 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 diff --git a/resolve/matrix/utilities.cpp b/resolve/matrix/utilities.cpp deleted file mode 100644 index 7ad866f5..00000000 --- a/resolve/matrix/utilities.cpp +++ /dev/null @@ -1,171 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include "utilities.hpp" - -namespace ReSolve -{ - using out = io::Logger; - 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 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 - * _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, memory::MemorySpace memspace) - { - //count nnzs first - index_type nnz_unpacked = 0; - 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]; ///< Number of elements per row - std::fill_n(nnz_counts, n, 0); - 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; - - // Count matrix elements - for (index_type i = 0; i < nnz; ++i) { - nnz_counts[coo_rows[i]]++; - nnz_unpacked++; - nnz_unpacked_no_duplicates++; - // 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 - nnz_unpacked_no_duplicates--; - nnz_no_duplicates--; - } - diag_control[coo_rows[i]]++; - } - } - 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]; - 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]; - - // 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) { - //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/tests/unit/matrix/CMakeLists.txt b/tests/unit/matrix/CMakeLists.txt index 5b939056..d313248a 100644 --- a/tests/unit/matrix/CMakeLists.txt +++ b/tests/unit/matrix/CMakeLists.txt @@ -18,6 +18,10 @@ target_link_libraries(runMatrixHandlerTests.exe PRIVATE ReSolve resolve_matrix) add_executable(runMatrixFactorizationTests.exe runMatrixFactorizationTests.cpp) target_link_libraries(runMatrixFactorizationTests.exe PRIVATE ReSolve resolve_matrix) +# Build matrix conversion tests +add_executable(runMatrixConversionTests.exe runMatrixConversionTests.cpp) +target_link_libraries(runMatrixConversionTests.exe PRIVATE ReSolve resolve_matrix) + # Build LUSOL-related tests if(RESOLVE_USE_LUSOL) add_executable(runLUSOLTests.exe runLUSOLTests.cpp) @@ -25,7 +29,7 @@ if(RESOLVE_USE_LUSOL) endif() # Install tests -set(installable_tests runMatrixIoTests.exe runMatrixHandlerTests.exe runMatrixFactorizationTests.exe) +set(installable_tests runMatrixIoTests.exe runMatrixHandlerTests.exe runMatrixFactorizationTests.exe runMatrixConversionTests.exe) if(RESOLVE_USE_LUSOL) list(APPEND installable_tests runLUSOLTests.exe) endif() @@ -35,6 +39,7 @@ install(TARGETS ${installable_tests} add_test(NAME matrix_test COMMAND $) add_test(NAME matrix_handler_test COMMAND $) add_test(NAME matrix_factorization_test COMMAND $) +add_test(NAME matrix_conversion_test COMMAND $) if(RESOLVE_USE_LUSOL) add_test(NAME lusol_test COMMAND $) endif() diff --git a/tests/unit/matrix/MatrixConversionTests.hpp b/tests/unit/matrix/MatrixConversionTests.hpp new file mode 100644 index 00000000..9fe5cc50 --- /dev/null +++ b/tests/unit/matrix/MatrixConversionTests.hpp @@ -0,0 +1,213 @@ +#pragma once + +#include +#include +#include +#include + +using index_type = ReSolve::index_type; +using real_type = ReSolve::real_type; + +namespace ReSolve +{ + namespace tests + { + /** + * @class Unit tests for matrix conversions + */ + class MatrixConversionTests : TestBase + { + public: + MatrixConversionTests() + { + } + + virtual ~MatrixConversionTests() + { + } + + TestOutcome simpleUpperUnexpandedSymmetricMatrix() + { + TestStatus status; + + matrix::Coo A(simple_upper_unexpanded_symmetric_n_, + simple_upper_unexpanded_symmetric_m_, + simple_upper_unexpanded_symmetric_nnz_, + true, + false); + A.updateData(simple_upper_unexpanded_symmetric_i_, + simple_upper_unexpanded_symmetric_j_, + simple_upper_unexpanded_symmetric_a_, + memory::HOST, + memory::HOST); + ReSolve::matrix::Csr B(A.getNumRows(), A.getNumColumns(), 0); + + status *= ReSolve::matrix::coo2csr(&A, &B, memory::HOST) == 0; + status *= this->verifyAnswer(&B, + simple_symmetric_expected_n_, + simple_symmetric_expected_m_, + simple_symmetric_expected_nnz_, + simple_symmetric_expected_i_, + simple_symmetric_expected_j_, + simple_symmetric_expected_a_); + + return status.report(__func__); + } + + TestOutcome simpleLowerUnexpandedSymmetricMatrix() + { + TestStatus status; + + matrix::Coo A(simple_lower_unexpanded_symmetric_n_, + simple_lower_unexpanded_symmetric_m_, + simple_lower_unexpanded_symmetric_nnz_, + true, + false); + A.updateData(simple_lower_unexpanded_symmetric_i_, + simple_lower_unexpanded_symmetric_j_, + simple_lower_unexpanded_symmetric_a_, + memory::HOST, + memory::HOST); + ReSolve::matrix::Csr B(A.getNumRows(), A.getNumColumns(), 0); + + status *= ReSolve::matrix::coo2csr(&A, &B, memory::HOST) == 0; + status *= this->verifyAnswer(&B, + simple_symmetric_expected_n_, + simple_symmetric_expected_m_, + simple_symmetric_expected_nnz_, + simple_symmetric_expected_i_, + simple_symmetric_expected_j_, + simple_symmetric_expected_a_); + + return status.report(__func__); + } + + TestOutcome simpleMainDiagonalOnlyMatrix() + { + TestStatus status; + + matrix::Coo A(simple_main_diagonal_only_n_, + simple_main_diagonal_only_m_, + simple_main_diagonal_only_nnz_, + true, + false); + A.updateData(simple_main_diagonal_only_i_j_, + simple_main_diagonal_only_i_j_, + simple_main_diagonal_only_a_, + memory::HOST, + memory::HOST); + ReSolve::matrix::Csr B(A.getNumRows(), A.getNumColumns(), 0); + + status *= ReSolve::matrix::coo2csr(&A, &B, memory::HOST) == 0; + status *= this->verifyAnswer(&B, + simple_main_diagonal_only_n_, + simple_main_diagonal_only_m_, + simple_main_diagonal_only_nnz_, + simple_main_diagonal_only_i_j_, + simple_main_diagonal_only_i_j_, + simple_main_diagonal_only_a_); + + return status.report(__func__); + } + + TestOutcome simpleFullAsymmetricMatrix() + { + TestStatus status; + + matrix::Coo A(simple_asymmetric_n_, + simple_asymmetric_m_, + simple_asymmetric_nnz_); + A.updateData(simple_asymmetric_i_, + simple_asymmetric_j_, + simple_asymmetric_a_, + memory::HOST, + memory::HOST); + ReSolve::matrix::Csr B(A.getNumRows(), A.getNumColumns(), 0); + + status *= ReSolve::matrix::coo2csr(&A, &B, memory::HOST) == 0; + status *= this->verifyAnswer(&B, + simple_asymmetric_expected_n_, + simple_asymmetric_expected_m_, + simple_asymmetric_expected_nnz_, + simple_asymmetric_expected_i_, + simple_asymmetric_expected_j_, + simple_asymmetric_expected_a_); + + return status.report(__func__); + } + + private: + const index_type simple_symmetric_expected_n_ = 5; + const index_type simple_symmetric_expected_m_ = 5; + const index_type simple_symmetric_expected_nnz_ = 8; + index_type simple_symmetric_expected_i_[8] = {0, 1, 1, 1, 2, 3, 3, 4}; + index_type simple_symmetric_expected_j_[8] = {0, 1, 2, 3, 1, 1, 4, 3}; + real_type simple_symmetric_expected_a_[8] = {2.0, 4.0, 6.0, 7.0, 6.0, 7.0, 8.0, 8.0}; + + const index_type simple_upper_unexpanded_symmetric_n_ = 5; + const index_type simple_upper_unexpanded_symmetric_m_ = 5; + const index_type simple_upper_unexpanded_symmetric_nnz_ = 8; + index_type simple_upper_unexpanded_symmetric_i_[8] = {0, 0, 1, 1, 1, 1, 1, 3}; + index_type simple_upper_unexpanded_symmetric_j_[8] = {0, 0, 1, 1, 2, 2, 3, 4}; + real_type simple_upper_unexpanded_symmetric_a_[8] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0}; + + const index_type simple_main_diagonal_only_n_ = 5; + const index_type simple_main_diagonal_only_m_ = 5; + const index_type simple_main_diagonal_only_nnz_ = 3; + index_type simple_main_diagonal_only_i_j_[3] = {1, 3, 4}; + real_type simple_main_diagonal_only_a_[3] = {1.0, 2.0, 3.0}; + + const index_type simple_lower_unexpanded_symmetric_n_ = 5; + const index_type simple_lower_unexpanded_symmetric_m_ = 5; + const index_type simple_lower_unexpanded_symmetric_nnz_ = 8; + index_type simple_lower_unexpanded_symmetric_i_[8] = {0, 0, 1, 1, 2, 2, 3, 4}; + index_type simple_lower_unexpanded_symmetric_j_[8] = {0, 0, 1, 1, 1, 1, 1, 3}; + real_type simple_lower_unexpanded_symmetric_a_[8] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0}; + + const index_type simple_asymmetric_n_ = 5; + const index_type simple_asymmetric_m_ = 5; + const index_type simple_asymmetric_nnz_ = 10; + index_type simple_asymmetric_i_[10] = {0, 1, 1, 2, 3, 3, 3, 4, 1, 1}; + index_type simple_asymmetric_j_[10] = {0, 1, 3, 1, 1, 4, 4, 3, 2, 2}; + real_type simple_asymmetric_a_[10] = {2.0, 4.0, 7.0, 9.0, 6.0, 7.0, 8.0, 8.0, 5.0, 6.0}; + + const index_type simple_asymmetric_expected_n_ = 5; + const index_type simple_asymmetric_expected_m_ = 5; + const index_type simple_asymmetric_expected_nnz_ = 8; + index_type simple_asymmetric_expected_i_[8] = {0, 1, 1, 1, 2, 3, 3, 4}; + index_type simple_asymmetric_expected_j_[8] = {0, 1, 2, 3, 1, 1, 4, 3}; + real_type simple_asymmetric_expected_a_[8] = {2.0, 4.0, 6.0, 7.0, 9.0, 6.0, 8.0, 8.0}; + + bool verifyAnswer(matrix::Csr* A, + const index_type& n, + const index_type& m, + const index_type& nnz, + index_type* is, + index_type* js, + real_type* as) + { + if (n != A->getNumRows() || m != A->getNumColumns() || nnz != A->getNnz()) { + return false; + } + + index_type* rows = A->getRowData(memory::HOST); + index_type* columns = A->getColData(memory::HOST); + real_type* values = A->getValues(memory::HOST); + + index_type answer_offset = 0; + for (index_type i = 0; i < A->getNumRows(); i++) { + for (index_type offset = rows[i]; offset < rows[i + 1]; offset++) { + if (i != is[answer_offset] || + columns[offset] != js[answer_offset] || + !isEqual(values[offset], as[answer_offset])) { + return false; + } + answer_offset++; + } + } + + return true; + } + }; + } // namespace tests +} // namespace ReSolve diff --git a/tests/unit/matrix/runMatrixConversionTests.cpp b/tests/unit/matrix/runMatrixConversionTests.cpp new file mode 100644 index 00000000..d8e9c156 --- /dev/null +++ b/tests/unit/matrix/runMatrixConversionTests.cpp @@ -0,0 +1,18 @@ +#include + +#include "MatrixConversionTests.hpp" + +int main() +{ + ReSolve::tests::TestingResults result; + ReSolve::tests::MatrixConversionTests test; + + result += test.simpleUpperUnexpandedSymmetricMatrix(); + result += test.simpleLowerUnexpandedSymmetricMatrix(); + result += test.simpleMainDiagonalOnlyMatrix(); + result += test.simpleFullAsymmetricMatrix(); + + std::cout << std::endl; + + return result.summary(); +}