diff --git a/examples/ExampleHelper.hpp b/examples/ExampleHelper.hpp index ec812c263..b3c1590bb 100644 --- a/examples/ExampleHelper.hpp +++ b/examples/ExampleHelper.hpp @@ -142,14 +142,10 @@ namespace ReSolve ReSolve::vector::Vector* r, ReSolve::vector::Vector* x) { + assert(res_ != nullptr && "resetSystem should be called after setSystem"); A_ = A; r_ = r; x_ = x; - if (res_ == nullptr) - { - res_ = new ReSolve::vector::Vector(A->getNumRows()); - } - computeNorms(); } diff --git a/examples/gluRefactor.cpp b/examples/gluRefactor.cpp index 322b392ba..3b1cbc6c9 100644 --- a/examples/gluRefactor.cpp +++ b/examples/gluRefactor.cpp @@ -264,7 +264,14 @@ int gluRefactor(int argc, char* argv[]) RESOLVE_RANGE_POP("Triangular solve"); // Print summary of the results - helper.resetSystem(A, vec_rhs, vec_x); + if (i == 0) + { + helper.setSystem(A, vec_rhs, vec_x); + } + else + { + helper.resetSystem(A, vec_rhs, vec_x); + } helper.printSummary(); } // for (int i = 0; i < num_systems; ++i) diff --git a/examples/gpuRefactor.cpp b/examples/gpuRefactor.cpp index 90f5ee06e..39f0dcc01 100644 --- a/examples/gpuRefactor.cpp +++ b/examples/gpuRefactor.cpp @@ -260,7 +260,7 @@ int gpuRefactor(int argc, char* argv[]) std::cout << "KLU solve status: " << status << std::endl; // Print summary of results - helper.resetSystem(A, vec_rhs, vec_x); + helper.setSystem(A, vec_rhs, vec_x); helper.printShortSummary(); if (i == 1) diff --git a/examples/kluFactor.cpp b/examples/kluFactor.cpp index 78ad54d50..e119af39d 100644 --- a/examples/kluFactor.cpp +++ b/examples/kluFactor.cpp @@ -184,8 +184,14 @@ int main(int argc, char* argv[]) status = KLU.solve(vec_rhs, vec_x); std::cout << "KLU solve status: " << status << std::endl; - - helper.resetSystem(A, vec_rhs, vec_x); + if (i == 0) + { + helper.setSystem(A, vec_rhs, vec_x); + } + else + { + helper.resetSystem(A, vec_rhs, vec_x); + } helper.printShortSummary(); if (is_iterative_refinement) { diff --git a/examples/kluRefactor.cpp b/examples/kluRefactor.cpp index ddf188d4a..8325b5e12 100644 --- a/examples/kluRefactor.cpp +++ b/examples/kluRefactor.cpp @@ -191,7 +191,14 @@ int main(int argc, char* argv[]) status = KLU->solve(vec_rhs, vec_x); std::cout << "KLU solve status: " << status << std::endl; - helper.resetSystem(A, vec_rhs, vec_x); + if (i == 0) + { + helper.setSystem(A, vec_rhs, vec_x); + } + else + { + helper.resetSystem(A, vec_rhs, vec_x); + } helper.printShortSummary(); if (is_iterative_refinement) { diff --git a/examples/randGmres.cpp b/examples/randGmres.cpp index 629cc2649..440b4031c 100644 --- a/examples/randGmres.cpp +++ b/examples/randGmres.cpp @@ -183,7 +183,7 @@ int runGmresExample(int argc, char* argv[]) FGMRES.solve(vec_rhs, vec_x); // Print summary of results - helper.resetSystem(A, vec_rhs, vec_x); + helper.setSystem(A, vec_rhs, vec_x); std::cout << "\nRandomized GMRES result on " << hwbackend << "\n"; std::cout << "---------------------------------\n"; helper.printIrSummary(&FGMRES); diff --git a/examples/sysRefactor.cpp b/examples/sysRefactor.cpp index 3a095213c..509a8efbc 100644 --- a/examples/sysRefactor.cpp +++ b/examples/sysRefactor.cpp @@ -320,7 +320,14 @@ int sysRefactor(int argc, char* argv[]) std::cout << "Triangular solve status: " << status << std::endl; // Print summary of results - helper.resetSystem(A, vec_rhs, vec_x); + if (i == 0) + { + helper.setSystem(A, vec_rhs, vec_x); + } + else + { + helper.resetSystem(A, vec_rhs, vec_x); + } helper.printShortSummary(); if ((i > 1) && is_iterative_refinement) { diff --git a/resolve/LinSolverDirect.cpp b/resolve/LinSolverDirect.cpp index d82e1046c..1f41cdb19 100644 --- a/resolve/LinSolverDirect.cpp +++ b/resolve/LinSolverDirect.cpp @@ -55,6 +55,7 @@ namespace ReSolve { if (A == nullptr) { + out::error() << "LinSolverDirect::setup: input matrix A is nullptr" << std::endl; return 1; } A_ = A; diff --git a/resolve/LinSolverDirectCuSolverRf.cpp b/resolve/LinSolverDirectCuSolverRf.cpp index d4f94429d..72b66d310 100644 --- a/resolve/LinSolverDirectCuSolverRf.cpp +++ b/resolve/LinSolverDirectCuSolverRf.cpp @@ -50,6 +50,7 @@ namespace ReSolve * * Sets up the cuSolverRf factorization for the given matrix A and its * L and U factors. The permutation vectors P and Q are also set up. + * This function should not be called more than once for the same object. * * @param[in] A - pointer to the matrix A * @param[in] L - pointer to the lower triangular factor L in CSR @@ -75,31 +76,41 @@ namespace ReSolve this->A_ = A; index_type n = A_->getNumRows(); - // Remember - P and Q are generally CPU variables! - // Factorization data is stored in the handle. - // If function is called again, destroy the old handle to get rid of old data. if (setup_completed_) { - cusolverRfDestroy(handle_cusolverrf_); - cusolverRfCreate(&handle_cusolverrf_); + out::error() << "Trying to setup LinSolverDirectCuSolverRf, but the setup has been already done! " << std::endl; + return 1; } if (d_P_ == nullptr) { mem_.allocateArrayOnDevice(&d_P_, n); } + else + { + out::error() << "Trying to allocate permutation vector P in " << __func__ << " in LinSolverDirectCuSolverRf, but the permutation vector P is already allocated! " << std::endl; + return 1; + } if (d_Q_ == nullptr) { mem_.allocateArrayOnDevice(&d_Q_, n); } - - if (d_T_ != nullptr) + else { - mem_.deleteOnDevice(d_T_); + out::error() << "Trying to allocate permutation vector Q in " << __func__ << " in LinSolverDirectCuSolverRf, but the permutation vector Q is already allocated! " << std::endl; + return 1; } - mem_.allocateArrayOnDevice(&d_T_, n); + if (d_T_ == nullptr) + { + mem_.allocateArrayOnDevice(&d_T_, n); + } + else + { + out::error() << "Trying to allocate temporary vector T in " << __func__ << " in LinSolverDirectCuSolverRf, but the temporary vector T is already allocated! " << std::endl; + return 1; + } mem_.copyArrayHostToDevice(d_P_, P, n); mem_.copyArrayHostToDevice(d_Q_, Q, n); diff --git a/resolve/LinSolverDirectKLU.cpp b/resolve/LinSolverDirectKLU.cpp index ecde2f48f..8a1f2ee93 100644 --- a/resolve/LinSolverDirectKLU.cpp +++ b/resolve/LinSolverDirectKLU.cpp @@ -127,7 +127,7 @@ namespace ReSolve if (Symbolic_ == nullptr) { out::error() << "Symbolic_ factorization failed with Common_.status = " - << Common_.status << "\n"; + << Common_.status << std::endl; return 1; } return 0; @@ -167,6 +167,8 @@ namespace ReSolve if (Numeric_ == nullptr) { + out::error() << "Numeric_ factorization failed with Common_.status = " + << Common_.status << std::endl; return 1; } else diff --git a/resolve/LinSolverDirectRocSolverRf.cpp b/resolve/LinSolverDirectRocSolverRf.cpp index 4f0ebc50e..c2481d6dd 100644 --- a/resolve/LinSolverDirectRocSolverRf.cpp +++ b/resolve/LinSolverDirectRocSolverRf.cpp @@ -82,11 +82,21 @@ namespace ReSolve { mem_.allocateArrayOnDevice(&d_P_, n); } + else + { + out::error() << "Trying to allocate permutation vector P in " << __func__ << " in LinSolverDirectRocSolverRf, but the permutation vector P is already allocated! " << std::endl; + return 1; + } if (d_Q_ == nullptr) { mem_.allocateArrayOnDevice(&d_Q_, n); } + else + { + out::error() << "Trying to allocate permutation vector Q in " << __func__ << " in LinSolverDirectRocSolverRf, but the permutation vector Q is already allocated! " << std::endl; + return 1; + } mem_.copyArrayHostToDevice(d_P_, P, n); mem_.copyArrayHostToDevice(d_Q_, Q, n); diff --git a/resolve/LinSolverIterative.cpp b/resolve/LinSolverIterative.cpp index c24af8f76..d94c5919b 100644 --- a/resolve/LinSolverIterative.cpp +++ b/resolve/LinSolverIterative.cpp @@ -26,6 +26,7 @@ namespace ReSolve { if (A == nullptr) { + out::error() << "LinSolverIterative::setup: input matrix A is nullptr" << std::endl; return 1; } this->A_ = A; diff --git a/resolve/MemoryUtils.hpp b/resolve/MemoryUtils.hpp index 702d80ea6..04f2ad92c 100644 --- a/resolve/MemoryUtils.hpp +++ b/resolve/MemoryUtils.hpp @@ -3,9 +3,12 @@ #include // <- declares `memcpy` #include +#include namespace ReSolve { + using out = ReSolve::io::Logger; + namespace memory { enum MemorySpace @@ -86,12 +89,22 @@ namespace ReSolve { std::size_t arraysize = static_cast(n) * sizeof(T); *v = new T[arraysize]; - return *v == nullptr ? 1 : 0; + if (*v == nullptr) + { + out::error() << "Memory allocation on host failed for size " << arraysize << " bytes." << std::endl; + return 1; + } + return 0; } template int deleteOnHost(T* v) { + if (v == nullptr) + { + out::error() << "Trying to delete nullptr on host!" << std::endl; + return 1; + } delete[] v; v = nullptr; return 0; diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp index 884c3fbd2..a714753cc 100644 --- a/resolve/SystemSolver.cpp +++ b/resolve/SystemSolver.cpp @@ -372,10 +372,6 @@ namespace ReSolve vectorHandler_, gs_); } - else - { - // do nothing - } return 0; } @@ -393,32 +389,46 @@ namespace ReSolve factorizationSolver_->setup(A_); return factorizationSolver_->analyze(); } + out::error() << "Analysis method " << factorizationSolver_ << " not recognized!\n" + << "Available options: klu.\n"; return 1; } int SystemSolver::factorize() { + if (A_ == nullptr) + { + out::error() << "System matrix not set!\n"; + return 1; + } if (factorizationMethod_ == "klu") { is_solve_on_device_ = false; return factorizationSolver_->factorize(); } + out::error() << "Factorization method << factorizationSolver_ << not recognized!\n" + << "Available options: klu.\n"; return 1; } int SystemSolver::refactorize() { + if (A_ == nullptr) + { + out::error() << "System matrix not set!\n"; + return 1; + } if (refactorizationMethod_ == "klu") { return factorizationSolver_->refactorize(); } - if (refactorizationMethod_ == "glu" || refactorizationMethod_ == "cusolverrf" || refactorizationMethod_ == "rocsolverrf") { is_solve_on_device_ = true; return refactorizationSolver_->refactorize(); } - + out::error() << "Refactorization method " << refactorizationSolver_ << " not recognized!\n" + << "Available options: klu, glu, cusolverrf, rocsolverrf.\n"; return 1; } diff --git a/resolve/matrix/io.cpp b/resolve/matrix/io.cpp index ae8502ae8..e6f9fab98 100644 --- a/resolve/matrix/io.cpp +++ b/resolve/matrix/io.cpp @@ -383,10 +383,14 @@ namespace ReSolve { if (!file) { - Logger::error() << "Empty input to updateArrayFromFile function ..." << std::endl; + Logger::error() << "Empty input to updateArrayFromFile function." << std::endl; + return; + } + if (p_rhs == nullptr) + { + Logger::error() << "Memory not allocated for updateArrayFromFile function." << std::endl; return; } - real_type* rhs = *p_rhs; std::stringstream ss; std::string line; @@ -399,11 +403,6 @@ namespace ReSolve } ss << line; ss >> n >> m; - - if (rhs == nullptr) - { - rhs = new real_type[n]; - } real_type a; index_type i = 0; while (file >> a) diff --git a/resolve/vector/Vector.cpp b/resolve/vector/Vector.cpp index 6d61d08ec..25412ad80 100644 --- a/resolve/vector/Vector.cpp +++ b/resolve/vector/Vector.cpp @@ -631,11 +631,11 @@ namespace ReSolve * * In case of multivectors, entire multivector is set to the constant. * - * @param[in] C - Constant (real number) - * @param[in] memspace - Memory space of the data to be set to 0 (HOST or DEVICE) + * @param[in] constant - Constant (real number) + * @param[in] memspace - Memory space of the data to be set to constant (HOST or DEVICE) * */ - int Vector::setToConst(real_type C, memory::MemorySpace memspace) + int Vector::setToConst(real_type constant, memory::MemorySpace memspace) { using namespace ReSolve::memory; switch (memspace) @@ -646,7 +646,7 @@ namespace ReSolve h_data_ = new real_type[n_capacity_ * k_]; owns_cpu_data_ = true; } - mem_.setArrayToConstOnHost(h_data_, C, n_size_ * k_); + mem_.setArrayToConstOnHost(h_data_, constant, n_size_ * k_); setHostUpdated(true); setDeviceUpdated(false); break; @@ -656,7 +656,7 @@ namespace ReSolve mem_.allocateArrayOnDevice(&d_data_, n_capacity_ * k_); owns_gpu_data_ = true; } - mem_.setArrayToConstOnDevice(d_data_, C, n_size_ * k_); + mem_.setArrayToConstOnDevice(d_data_, constant, n_size_ * k_); setHostUpdated(false); setDeviceUpdated(true); break; @@ -668,12 +668,12 @@ namespace ReSolve * @brief set the data of a single vector in a multivector to a given constant. * * @param[in] j - Index of a vector in a multivector - * @param[in] C - Constant (real number) + * @param[in] constant - Constant (real number) * @param[in] memspace - Memory space of the data to be set to 0 (HOST or DEVICE) * * @pre _j_ < _k_ i.e,, _j_ is smaller than the total number of vectors in multivector. */ - int Vector::setToConst(index_type j, real_type C, memory::MemorySpace memspace) + int Vector::setToConst(index_type j, real_type constant, memory::MemorySpace memspace) { using namespace ReSolve::memory; switch (memspace) @@ -681,20 +681,22 @@ namespace ReSolve case HOST: if (h_data_ == nullptr) { - h_data_ = new real_type[n_capacity_ * k_]; - owns_cpu_data_ = true; + out::error() << "In Vector setToConst: trying to set host data to constant, " + << "but host data not allocated!" << std::endl; + return 1; } - mem_.setArrayToConstOnHost(&h_data_[n_size_ * j], C, n_size_); + mem_.setArrayToConstOnHost(&h_data_[n_size_ * j], constant, n_size_); cpu_updated_[j] = true; gpu_updated_[j] = false; break; case DEVICE: if (d_data_ == nullptr) { - mem_.allocateArrayOnDevice(&d_data_, n_capacity_ * k_); - owns_gpu_data_ = true; + out::error() << "In Vector setToConst: trying to set device data to constant, " + << "but device data not allocated!" << std::endl; + return 1; } - mem_.setArrayToConstOnDevice(&d_data_[n_size_ * j], C, n_size_); + mem_.setArrayToConstOnDevice(&d_data_[n_size_ * j], constant, n_size_); cpu_updated_[j] = false; gpu_updated_[j] = true; break; diff --git a/tests/unit/matrix/MatrixHandlerTests.hpp b/tests/unit/matrix/MatrixHandlerTests.hpp index e6fb0eac9..e624ad6e7 100644 --- a/tests/unit/matrix/MatrixHandlerTests.hpp +++ b/tests/unit/matrix/MatrixHandlerTests.hpp @@ -8,12 +8,15 @@ #include #include #include +#include #include #include #include namespace ReSolve { + using out = io::Logger; + namespace tests { @@ -616,6 +619,7 @@ namespace ReSolve // Check if the matrix is valid if (A == nullptr) { + out::error() << "Matrix pointer is NULL in " << __func__ << std::endl; return false; } @@ -697,6 +701,7 @@ namespace ReSolve // Check if the matrix is valid if (A == nullptr) { + out::error() << "Matrix pointer is NULL in " << __func__ << std::endl; return false; } @@ -780,6 +785,7 @@ namespace ReSolve // Check if the vector is valid if (scaled_vec == nullptr) { + out::error() << "Vector pointer is NULL in " << __func__ << std::endl; return false; }