diff --git a/resolve/matrix/Coo.cpp b/resolve/matrix/Coo.cpp index c697d4c3..3555d68b 100644 --- a/resolve/matrix/Coo.cpp +++ b/resolve/matrix/Coo.cpp @@ -56,7 +56,8 @@ namespace ReSolve h_col_data_ = *cols; h_val_data_ = *vals; h_data_updated_ = true; - owns_cpu_vals_ = owns_cpu_data_ = true; + owns_cpu_vals_ = true; + owns_cpu_data_ = true; // Set device data to null if (d_row_data_ || d_col_data_ || d_val_data_) { out::error() << "Device data unexpectedly allocated. " @@ -66,7 +67,8 @@ namespace ReSolve d_col_data_ = nullptr; d_val_data_ = nullptr; d_data_updated_ = false; - owns_gpu_vals_ = owns_gpu_data_ = false; + owns_gpu_vals_ = true; + owns_gpu_data_ = false; // Hijack data from the source *rows = nullptr; *cols = nullptr; @@ -78,7 +80,8 @@ namespace ReSolve d_col_data_ = *cols; d_val_data_ = *vals; d_data_updated_ = true; - owns_gpu_vals_ = owns_gpu_data_ = true; + owns_gpu_vals_ = true; + owns_gpu_data_ = true; copyData(memspaceDst); // Hijack data from the source *rows = nullptr; @@ -91,7 +94,8 @@ namespace ReSolve h_col_data_ = *cols; h_val_data_ = *vals; h_data_updated_ = true; - owns_cpu_vals_ = owns_cpu_data_ = true; + owns_cpu_vals_ = true; + owns_cpu_data_ = true; copyData(memspaceDst); // Hijack data from the source @@ -105,7 +109,8 @@ namespace ReSolve d_col_data_ = *cols; d_val_data_ = *vals; d_data_updated_ = true; - owns_gpu_vals_ = owns_gpu_data_ = true; + owns_gpu_vals_ = true; + owns_gpu_data_ = true; // Set host data to null if (h_row_data_ || h_col_data_ || h_val_data_) { out::error() << "Host data unexpectedly allocated. " @@ -115,7 +120,8 @@ namespace ReSolve h_col_data_ = nullptr; h_val_data_ = nullptr; h_data_updated_ = false; - owns_cpu_vals_ = owns_cpu_data_ = false; + owns_cpu_vals_ = true; + owns_cpu_data_ = false; // Hijack data from the source *rows = nullptr; *cols = nullptr; diff --git a/resolve/matrix/Csr.cpp b/resolve/matrix/Csr.cpp index 9a5fe717..78c6502e 100644 --- a/resolve/matrix/Csr.cpp +++ b/resolve/matrix/Csr.cpp @@ -8,7 +8,7 @@ #include #include -namespace ReSolve +namespace ReSolve { using out = io::Logger; @@ -19,9 +19,9 @@ namespace ReSolve matrix::Csr::Csr(index_type n, index_type m, index_type nnz) : Sparse(n, m, nnz) { } - - matrix::Csr::Csr(index_type n, - index_type m, + + matrix::Csr::Csr(index_type n, + index_type m, index_type nnz, bool symmetric, bool expanded) : Sparse(n, m, nnz, symmetric, expanded) @@ -40,20 +40,20 @@ namespace ReSolve /** * @brief Hijacking constructor - * - * @param[in] n - * @param[in] m - * @param[in] nnz - * @param[in] symmetric - * @param[in] expanded - * @param[in,out] rows - * @param[in,out] cols - * @param[in,out] vals - * @param[in] memspaceSrc - * @param[in] memspaceDst + * + * @param[in] n + * @param[in] m + * @param[in] nnz + * @param[in] symmetric + * @param[in] expanded + * @param[in,out] rows + * @param[in,out] cols + * @param[in,out] vals + * @param[in] memspaceSrc + * @param[in] memspaceDst */ - matrix::Csr::Csr(index_type n, - index_type m, + matrix::Csr::Csr(index_type n, + index_type m, index_type nnz, bool symmetric, bool expanded, @@ -79,7 +79,8 @@ namespace ReSolve h_col_data_ = *cols; h_val_data_ = *vals; h_data_updated_ = true; - owns_cpu_vals_ = owns_cpu_data_ = true; + owns_cpu_vals_ = true; + owns_cpu_data_ = true; // Set device data to null if (d_row_data_ || d_col_data_ || d_val_data_) { out::error() << "Device data unexpectedly allocated. " @@ -89,7 +90,8 @@ namespace ReSolve d_col_data_ = nullptr; d_val_data_ = nullptr; d_data_updated_ = false; - owns_gpu_vals_ = owns_gpu_data_ = false; + owns_gpu_vals_ = true; + owns_gpu_data_ = false; // Hijack data from the source *rows = nullptr; *cols = nullptr; @@ -101,7 +103,8 @@ namespace ReSolve d_col_data_ = *cols; d_val_data_ = *vals; d_data_updated_ = true; - owns_gpu_vals_ = owns_gpu_data_ = true; + owns_gpu_vals_ = true; + owns_gpu_data_ = true; copyData(memspaceDst); // Hijack data from the source *rows = nullptr; @@ -114,7 +117,8 @@ namespace ReSolve h_col_data_ = *cols; h_val_data_ = *vals; h_data_updated_ = true; - owns_cpu_vals_ = owns_cpu_data_ = true; + owns_cpu_vals_ = true; + owns_cpu_data_ = true; copyData(memspaceDst); // Hijack data from the source @@ -128,7 +132,8 @@ namespace ReSolve d_col_data_ = *cols; d_val_data_ = *vals; d_data_updated_ = true; - owns_gpu_vals_ = owns_gpu_data_ = true; + owns_gpu_vals_ = true; + owns_gpu_data_ = true; // Set host data to null if (h_row_data_ || h_col_data_ || h_val_data_) { out::error() << "Host data unexpectedly allocated. " @@ -138,7 +143,8 @@ namespace ReSolve h_col_data_ = nullptr; h_val_data_ = nullptr; h_data_updated_ = false; - owns_cpu_vals_ = owns_cpu_data_ = false; + owns_cpu_vals_ = true; + owns_cpu_data_ = false; // Hijack data from the source *rows = nullptr; *cols = nullptr; @@ -218,7 +224,7 @@ namespace ReSolve this->h_row_data_ = new index_type[n_ + 1]; this->h_col_data_ = new index_type[nnz_current]; owns_cpu_data_ = true; - } + } if (h_val_data_ == nullptr) { this->h_val_data_ = new real_type[nnz_current]; owns_cpu_vals_ = true; @@ -231,18 +237,18 @@ namespace ReSolve out::error() << "In Csr::updateData one of device row or column data is null!\n"; } if ((d_row_data_ == nullptr) && (d_col_data_ == nullptr)) { - mem_.allocateArrayOnDevice(&d_row_data_, n_ + 1); + mem_.allocateArrayOnDevice(&d_row_data_, n_ + 1); mem_.allocateArrayOnDevice(&d_col_data_, nnz_current); owns_gpu_vals_ = true; } if (d_val_data_ == nullptr) { - mem_.allocateArrayOnDevice(&d_val_data_, nnz_current); + mem_.allocateArrayOnDevice(&d_val_data_, nnz_current); owns_gpu_data_ = true; } } - //copy + //copy switch(control) { case 0: //cpu->cpu mem_.copyArrayHostToHost(h_row_data_, row_data, n_ + 1); @@ -272,7 +278,7 @@ namespace ReSolve return -1; } return 0; - } + } int matrix::Csr::updateData(index_type* row_data, index_type* col_data, real_type* val_data, index_type new_nnz, memory::MemorySpace memspaceIn, memory::MemorySpace memspaceOut) { @@ -280,7 +286,7 @@ namespace ReSolve this->nnz_ = new_nnz; int i = this->updateData(row_data, col_data, val_data, memspaceIn, memspaceOut); return i; - } + } int matrix::Csr::allocateMatrixData(memory::MemorySpace memspace) { @@ -290,23 +296,23 @@ namespace ReSolve if (memspace == memory::HOST) { this->h_row_data_ = new index_type[n_ + 1]; - std::fill(h_row_data_, h_row_data_ + n_ + 1, 0); + std::fill(h_row_data_, h_row_data_ + n_ + 1, 0); this->h_col_data_ = new index_type[nnz_current]; - std::fill(h_col_data_, h_col_data_ + nnz_current, 0); + std::fill(h_col_data_, h_col_data_ + nnz_current, 0); this->h_val_data_ = new real_type[nnz_current]; - std::fill(h_val_data_, h_val_data_ + nnz_current, 0.0); + std::fill(h_val_data_, h_val_data_ + nnz_current, 0.0); owns_cpu_data_ = true; owns_cpu_vals_ = true; - return 0; + return 0; } if (memspace == memory::DEVICE) { - mem_.allocateArrayOnDevice(&d_row_data_, n_ + 1); - mem_.allocateArrayOnDevice(&d_col_data_, nnz_current); - mem_.allocateArrayOnDevice(&d_val_data_, nnz_current); + mem_.allocateArrayOnDevice(&d_row_data_, n_ + 1); + mem_.allocateArrayOnDevice(&d_col_data_, nnz_current); + mem_.allocateArrayOnDevice(&d_val_data_, nnz_current); owns_gpu_data_ = true; owns_gpu_vals_ = true; - return 0; + return 0; } return -1; } @@ -329,11 +335,11 @@ namespace ReSolve } if ((h_row_data_ == nullptr) && (h_col_data_ == nullptr)) { h_row_data_ = new index_type[n_ + 1]; - h_col_data_ = new index_type[nnz_current]; + h_col_data_ = new index_type[nnz_current]; owns_cpu_data_ = true; } if (h_val_data_ == nullptr) { - h_val_data_ = new real_type[nnz_current]; + h_val_data_ = new real_type[nnz_current]; owns_cpu_vals_ = true; } mem_.copyArrayDeviceToHost(h_row_data_, d_row_data_, n_ + 1); @@ -348,12 +354,12 @@ namespace ReSolve out::error() << "In Csr::copyData one of device row or column data is null!\n"; } if ((d_row_data_ == nullptr) && (d_col_data_ == nullptr)) { - mem_.allocateArrayOnDevice(&d_row_data_, n_ + 1); - mem_.allocateArrayOnDevice(&d_col_data_, nnz_current); + mem_.allocateArrayOnDevice(&d_row_data_, n_ + 1); + mem_.allocateArrayOnDevice(&d_col_data_, nnz_current); owns_gpu_data_ = true; } if (d_val_data_ == nullptr) { - mem_.allocateArrayOnDevice(&d_val_data_, nnz_current); + mem_.allocateArrayOnDevice(&d_val_data_, nnz_current); owns_gpu_vals_ = true; } mem_.copyArrayHostToDevice(d_row_data_, h_row_data_, n_ + 1); @@ -429,7 +435,7 @@ namespace ReSolve index_type* nnz_shifts = new index_type[n]; std::fill_n(nnz_shifts, n , 0); - IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; + IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; csr_ia[0] = 0; @@ -449,7 +455,7 @@ namespace ReSolve 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; + bool already_there = false; for (index_type j = start; j < start + nnz_shifts[r]; ++j) { index_type c = tmp[j].getIdx(); @@ -459,8 +465,8 @@ namespace ReSolve 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]); @@ -512,14 +518,14 @@ namespace ReSolve delete [] csr_ia; delete [] csr_ja; delete [] csr_a; - delete [] diag_control; + delete [] diag_control; return 0; } /** * @brief Prints matrix data. - * + * * @param out - Output stream where the matrix data is printed */ void matrix::Csr::print(std::ostream& out) @@ -527,11 +533,10 @@ namespace ReSolve out << std::scientific << std::setprecision(std::numeric_limits::digits10); for(index_type i = 0; i < n_; ++i) { for (index_type j = h_row_data_[i]; j < h_row_data_[i+1]; ++j) { - out << i << " " + out << i << " " << h_col_data_[j] << " " << h_val_data_[j] << "\n"; } } } -} // namespace ReSolve - +} // namespace ReSolve diff --git a/resolve/matrix/Utilities.cpp b/resolve/matrix/Utilities.cpp index ba301076..78186d2c 100644 --- a/resolve/matrix/Utilities.cpp +++ b/resolve/matrix/Utilities.cpp @@ -10,7 +10,9 @@ namespace ReSolve { /** - * @invariant The input matrix must be deduplicated, otherwise the result is undefined + * @pre The input matrix must be deduplicated + * @pre If the input matrix is symmetric and unexpanded, it is either upper + * or lower triangular */ int matrix::expand(matrix::Coo& A) { @@ -75,7 +77,9 @@ namespace ReSolve } /** - * @invariant The input matrix must be deduplicated, otherwise the result is undefined + * @pre The input matrix must be deduplicated + * @pre If the input matrix is symmetric and unexpanded, it is either upper + * or lower triangular */ int matrix::expand(matrix::Csr& A) { @@ -149,7 +153,9 @@ namespace ReSolve } /** - * @invariant The input matrix must be deduplicated, otherwise the result is undefined + * @pre The input matrix must be deduplicated + * @pre If the input matrix is symmetric and unexpanded, it is either upper + * or lower triangular */ int matrix::expand(matrix::Csc& A) { diff --git a/tests/functionality/testLUSOL.cpp b/tests/functionality/testLUSOL.cpp index 17f1609f..b0e7050a 100644 --- a/tests/functionality/testLUSOL.cpp +++ b/tests/functionality/testLUSOL.cpp @@ -64,7 +64,7 @@ int main(int argc, char* argv[]) return 1; } - std::unique_ptr rhs(ReSolve::io::readRhsFromFile(rhs_file)); + real_type* rhs = ReSolve::io::readRhsFromFile(rhs_file); rhs_file.close(); std::unique_ptr x(new real_type[A->getNumRows()]); @@ -72,7 +72,7 @@ int main(int argc, char* argv[]) std::unique_ptr vec_x(new vector_type(A->getNumRows())); std::unique_ptr vec_r(new vector_type(A->getNumRows())); - vec_rhs->update(rhs.get(), ReSolve::memory::HOST, ReSolve::memory::HOST); + vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); vec_x->allocate(ReSolve::memory::HOST); @@ -99,7 +99,7 @@ int main(int argc, char* argv[]) std::fill_n(x_data, A->getNumRows(), 1.0); vec_test->setData(x_data, ReSolve::memory::HOST); - vec_r->update(rhs.get(), ReSolve::memory::HOST, ReSolve::memory::HOST); + vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_diff->update(x_data, ReSolve::memory::HOST, ReSolve::memory::HOST); matrix_handler->setValuesChanged(true, ReSolve::memory::HOST); @@ -117,7 +117,7 @@ int main(int argc, char* argv[]) real_type normDiffMatrix = sqrt(vector_handler->dot(vec_diff.get(), vec_diff.get(), ReSolve::memory::HOST)); // compute the residual using exact solution - vec_r->update(rhs.get(), ReSolve::memory::HOST, ReSolve::memory::HOST); + vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); error_sum += specializedMatvec(A.get(), vec_test.get(), vec_r.get(), @@ -126,7 +126,7 @@ int main(int argc, char* argv[]) real_type exactSol_normRmatrix = sqrt(vector_handler->dot(vec_r.get(), vec_r.get(), ReSolve::memory::HOST)); // evaluate the residual ON THE CPU using COMPUTED solution - vec_r->update(rhs.get(), ReSolve::memory::HOST, ReSolve::memory::HOST); + vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); error_sum += specializedMatvec(A.get(), vec_x.get(), @@ -144,6 +144,7 @@ int main(int argc, char* argv[]) std::cout << "\t ||x-x_true||_2/||x_true||_2 : " << normDiffMatrix / normXtrue << " (scaled solution error)\n"; std::cout << "\t ||b-A*x_exact||_2 : " << exactSol_normRmatrix << " (control; residual norm with exact solution)\n\n"; + delete[] rhs; delete[] x_data; real_type scaled_residual_norm_one = normRmatrix / normB; @@ -164,7 +165,7 @@ int main(int argc, char* argv[]) return 1; } - rhs = std::unique_ptr(ReSolve::io::readRhsFromFile(rhs_file)); + rhs = ReSolve::io::readRhsFromFile(rhs_file); rhs_file.close(); x = std::unique_ptr(new real_type[A->getNumRows()]); @@ -172,7 +173,7 @@ int main(int argc, char* argv[]) vec_x = std::unique_ptr(new vector_type(A->getNumRows())); vec_r = std::unique_ptr(new vector_type(A->getNumRows())); - vec_rhs->update(rhs.get(), ReSolve::memory::HOST, ReSolve::memory::HOST); + vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_rhs->setDataUpdated(ReSolve::memory::HOST); vec_x->allocate(ReSolve::memory::HOST); @@ -199,7 +200,7 @@ int main(int argc, char* argv[]) std::fill_n(x_data, A->getNumRows(), 1.0); vec_test->setData(x_data, ReSolve::memory::HOST); - vec_r->update(rhs.get(), ReSolve::memory::HOST, ReSolve::memory::HOST); + vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); vec_diff->update(x_data, ReSolve::memory::HOST, ReSolve::memory::HOST); matrix_handler->setValuesChanged(true, ReSolve::memory::HOST); @@ -217,7 +218,7 @@ int main(int argc, char* argv[]) normDiffMatrix = sqrt(vector_handler->dot(vec_diff.get(), vec_diff.get(), ReSolve::memory::HOST)); // compute the residual using exact solution - vec_r->update(rhs.get(), ReSolve::memory::HOST, ReSolve::memory::HOST); + vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); error_sum += specializedMatvec(A.get(), vec_test.get(), vec_r.get(), @@ -226,7 +227,7 @@ int main(int argc, char* argv[]) exactSol_normRmatrix = sqrt(vector_handler->dot(vec_r.get(), vec_r.get(), ReSolve::memory::HOST)); // evaluate the residual ON THE CPU using COMPUTED solution - vec_r->update(rhs.get(), ReSolve::memory::HOST, ReSolve::memory::HOST); + vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); error_sum += specializedMatvec(A.get(), vec_x.get(), @@ -244,6 +245,7 @@ int main(int argc, char* argv[]) std::cout << "\t ||x-x_true||_2/||x_true||_2 : " << normDiffMatrix / normXtrue << " (scaled solution error)\n"; std::cout << "\t ||b-A*x_exact||_2 : " << exactSol_normRmatrix << " (control; residual norm with exact solution)\n\n"; + delete[] rhs; delete[] x_data; real_type scaled_residual_norm_two = normRmatrix / normB; diff --git a/tests/unit/matrix/MatrixExpansionTests.hpp b/tests/unit/matrix/MatrixExpansionTests.hpp index 5772d731..168eaa8e 100644 --- a/tests/unit/matrix/MatrixExpansionTests.hpp +++ b/tests/unit/matrix/MatrixExpansionTests.hpp @@ -16,6 +16,9 @@ namespace ReSolve { namespace tests { + /** + * @class Unit tests for `ReSolve::matrix::expand` + */ class MatrixExpansionTests : TestBase { public: @@ -27,44 +30,50 @@ namespace ReSolve { } - TestOutcome cooMatrix5x5() + /** + * @brief Correctness test of `ReSolve::matrix::expand` on a COO + * matrix + */ + TestOutcome cooMatrix() { TestStatus status; std::unique_ptr A = buildCooMatrix5x5(); matrix::expand(*A); - A->print(); - - status *= validateAnswer(*A, target_triples_5x5_); + status *= verifyAnswer(*A, target_triples_5x5_); return status.report(__func__); } - TestOutcome csrMatrix5x5() + /** + * @brief Correctness test of `ReSolve::matrix::expand` on a CSR + * matrix + */ + TestOutcome csrMatrix() { TestStatus status; std::unique_ptr A = buildCsrMatrix5x5(); matrix::expand(*A); - A->print(); - - status *= validateAnswer(*A, target_triples_5x5_); + status *= verifyAnswer(*A, target_triples_5x5_); return status.report(__func__); } - TestOutcome cscMatrix5x5() + /** + * @brief Correctness test of `ReSolve::matrix::expand` on a CSC + * matrix + */ + TestOutcome cscMatrix() { TestStatus status; std::unique_ptr A = buildCscMatrix5x5(); matrix::expand(*A); - A->print(); - - status *= validateAnswer(*A, target_triples_5x5_); + status *= verifyAnswer(*A, target_triples_5x5_); return status.report(__func__); } @@ -77,8 +86,8 @@ namespace ReSolve {0, 4, 3.0}, {4, 0, 3.0}}; - bool validateAnswer(matrix::Coo& A, - std::vector> target) + bool verifyAnswer(matrix::Coo& A, + std::vector> target) { std::shared_ptr i(new index_type(0)); std::shared_ptr nnz(new index_type(A.getNnz())); @@ -90,7 +99,7 @@ namespace ReSolve return false; } - return validateAnswer( + return verifyAnswer( [=]() -> std::tuple, bool> { if (*i == *nnz) { return {{0, 0, 0}, false}; @@ -102,7 +111,7 @@ namespace ReSolve target); } - bool validateAnswer(matrix::Csr& A, + bool verifyAnswer(matrix::Csr& A, std::vector> target) { std::shared_ptr i(new index_type(0)), j(new index_type(0)); @@ -115,7 +124,7 @@ namespace ReSolve return false; } - return validateAnswer( + return verifyAnswer( [=]() -> std::tuple, bool> { if (*j == *nnz) { return {{0, 0, 0}, false}; @@ -131,7 +140,7 @@ namespace ReSolve target); } - bool validateAnswer(matrix::Csc& A, + bool verifyAnswer(matrix::Csc& A, std::vector> target) { std::shared_ptr i(new index_type(0)), j(new index_type(0)); @@ -144,7 +153,7 @@ namespace ReSolve return false; } - return validateAnswer( + return verifyAnswer( [=]() -> std::tuple, bool> { if (*j == *nnz) { return {{0, 0, 0}, false}; @@ -160,7 +169,7 @@ namespace ReSolve target); } - bool validateAnswer(std::function, bool>()> f, + bool verifyAnswer(std::function, bool>()> f, std::vector> target) { bool ok; diff --git a/tests/unit/matrix/runMatrixExpansionTests.cpp b/tests/unit/matrix/runMatrixExpansionTests.cpp index 33b9c783..59979ef8 100644 --- a/tests/unit/matrix/runMatrixExpansionTests.cpp +++ b/tests/unit/matrix/runMatrixExpansionTests.cpp @@ -15,9 +15,9 @@ int main(int, char**) std::cout << "Running tests on CPU:\n"; ReSolve::tests::MatrixExpansionTests test; - result += test.cooMatrix5x5(); - result += test.csrMatrix5x5(); - result += test.cscMatrix5x5(); + result += test.cooMatrix(); + result += test.csrMatrix(); + result += test.cscMatrix(); std::cout << "\n"; }