diff --git a/CMakeLists.txt b/CMakeLists.txt index c45bd56f..f4fb5a0a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,7 +11,7 @@ cmake_minimum_required(VERSION 3.22) # Adds version settings and set variable CMAKE_PROJECT_VERSION project(ReSolve VERSION "0.99.1") -set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD 14) set(PACKAGE_NAME "ReSolve") set(PACKAGE_TARNAME "resolve") diff --git a/tests/functionality/CMakeLists.txt b/tests/functionality/CMakeLists.txt index b4dc31bc..8ef4c630 100644 --- a/tests/functionality/CMakeLists.txt +++ b/tests/functionality/CMakeLists.txt @@ -6,6 +6,10 @@ ]] +# Create functionality test helper +add_library(FunctionalityTestHelper SHARED FunctionalityTestHelper.cpp) +target_link_libraries(FunctionalityTestHelper PRIVATE ReSolve) + # Build basic version test add_executable(version.exe testVersion.cpp) target_link_libraries(version.exe PRIVATE ReSolve) @@ -38,7 +42,7 @@ if(RESOLVE_USE_CUDA) # System solver test with cusolver rf and iterative refinement add_executable(sys_refactor_cuda_test.exe testSysRefactor.cpp) - target_link_libraries(sys_refactor_cuda_test.exe PRIVATE ReSolve) + target_link_libraries(sys_refactor_cuda_test.exe PRIVATE ReSolve FunctionalityTestHelper) # Build KLU+GLU test add_executable(klu_glu_test.exe testKLU_GLU.cpp) @@ -69,7 +73,7 @@ if(RESOLVE_USE_HIP) # System solver test with rocm rf and iterative refinement add_executable(sys_refactor_hip_test.exe testSysRefactor.cpp) - target_link_libraries(sys_refactor_hip_test.exe PRIVATE ReSolve) + target_link_libraries(sys_refactor_hip_test.exe PRIVATE ReSolve FunctionalityTestHelper ) endif(RESOLVE_USE_KLU) # Build randSolver test diff --git a/tests/functionality/FunctionalityTestHelper.cpp b/tests/functionality/FunctionalityTestHelper.cpp new file mode 100644 index 00000000..51b73c89 --- /dev/null +++ b/tests/functionality/FunctionalityTestHelper.cpp @@ -0,0 +1,433 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include "resolve/Common.hpp" +#include +#include +#include +#if defined RESOLVE_USE_KLU +#include +#endif +#if defined (RESOLVE_USE_CUDA) +#include + using workspace_type = ReSolve::LinAlgWorkspaceCUDA; + std::string memory_space("cuda"); +#elif defined (RESOLVE_USE_HIP) +#include + using workspace_type = ReSolve::LinAlgWorkspaceHIP; + std::string memory_space("hip"); +#else + using workspace_type = ReSolve::LinAlgWorkspaceCpu; + std::string memory_space("cpu"); +#endif +#include "FunctionalityTestHelper.hpp" + +namespace ReSolve +{ + + namespace tests + { + + AxEqualsRhsProblem::~AxEqualsRhsProblem() + { + delete A_; + delete vec_x_; + delete vec_rhs_; + } + + matrix::Csr* AxEqualsRhsProblem::getMatrix() const + { + return A_; + } + + vector::Vector* AxEqualsRhsProblem::getVector() const + { + return vec_x_; + } + + vector::Vector* AxEqualsRhsProblem::getRhs() const + { + return vec_rhs_; + } + + void AxEqualsRhsProblem::updateProblem(const std::string& matrix_filepath, + const std::string& rhs_filepath) + { + // Load the second matrix + std::ifstream mat2(matrix_filepath); + if(!mat2.is_open()) { + std::cout << "Failed to open file " << matrix_filepath << "\n"; + std::exit( 1 ); + } + + io::updateMatrixFromFile(mat2, A_); + + A_->syncData(memory::DEVICE); + + mat2.close(); + + // Load the second rhs vector + std::ifstream rhs2_file(rhs_filepath); + if(!rhs2_file.is_open()) { + std::cout << "Failed to open file " << rhs_filepath << "\n"; + std::exit( 1 ); + } + + // note: This intermediate allocation is clunky and should not be here + // in the original test, it was updateArrayFromFile: but we should not have to + // hold onto the intermediate rhs since it is just a temporary construction artifact + real_type* rhs = io::createArrayFromFile(rhs2_file); + //io::updateArrayFromFile(rhs2_file, &rhs); + + rhs2_file.close(); + + vec_rhs_->update(rhs, memory::HOST, memory::DEVICE); + + delete[] rhs; + } + + AxEqualsRhsProblem::AxEqualsRhsProblem(const std::string& matrix_filepath, + const std::string& rhs_filepath) + { + // Read first matrix + std::ifstream mat1(matrix_filepath); + if(!mat1.is_open()) { + std::cout << "Failed to open file " << matrix_filepath << "\n"; + std::exit(1); + } + + A_ = io::createCsrFromFile(mat1, true); + A_->syncData(memory::DEVICE); + mat1.close(); + + // Read first rhs vector + std::ifstream rhs1_file(rhs_filepath); + if(!rhs1_file.is_open()) { + std::cout << "Failed to open file " << rhs_filepath << "\n"; + std::exit(1); + } + real_type* rhs = io::createArrayFromFile(rhs1_file); + rhs1_file.close(); + + // setup/allocate testing workspace phase: + + // Create rhs, solution and residual vectors + vec_rhs_ = new vector::Vector(A_->getNumRows()); + vec_x_ = new vector::Vector(A_->getNumRows()); + + // Allocate solution vector + vec_x_->allocate(memory::HOST); //for KLU + vec_x_->allocate(memory::DEVICE); + + // Set RHS vector on CPU (update function allocates) + vec_rhs_->update(rhs, memory::HOST, memory::HOST); + + delete[] rhs; + } + + real_type FunctionalityTestHelper::calculateSolutionVectorNorm(vector::Vector& vec_x) + { + using namespace memory; + + // set member variable and also return in case this function is used outside of this class + norm_solution_vector_ = sqrt(vh_->dot(&vec_x, &vec_x, DEVICE)); + + return norm_solution_vector_; + } + + real_type FunctionalityTestHelper::calculateRhsVectorNorm(vector::Vector& vec_rhs) + { + using namespace memory; + + // set member variable and also return in case this function is used outside of this class + rhs_norm_ = sqrt(vh_->dot(&vec_rhs, &vec_rhs, DEVICE)); + + return rhs_norm_; + } + + real_type FunctionalityTestHelper::calculateResidualNorm(vector::Vector& vec_r) + { + using namespace memory; + + // set member variable and also return in case this function is used outside of this class + residual_norm_ = sqrt(vh_->dot(&vec_r, &vec_r, DEVICE)); + + return residual_norm_; + } + + real_type FunctionalityTestHelper::calculateDiffNorm(vector::Vector& vec_x) + { + using namespace memory; + using namespace ReSolve::constants; + + index_type n = vec_x.getSize(); + + vector::Vector vec_diff(n); + + vec_diff.setToConst(1.0, DEVICE); + + // why does this not return an error status if it fails? + vh_->axpy(&MINUSONE, &vec_x, &vec_diff, DEVICE); + + diff_norm_ = sqrt(vh_->dot(&vec_diff, &vec_diff, DEVICE)); + + return diff_norm_; + } + + real_type FunctionalityTestHelper::calculateTrueNorm(matrix::Csr& A, vector::Vector& vec_rhs) + { + using namespace memory; + using namespace ReSolve::constants; + + index_type n = A.getNumRows(); + + vector::Vector vec_test(n); + vector::Vector vec_tmp(n); + + //compute the residual using exact solution + vec_test.setToConst(1.0, DEVICE); + vec_tmp.update(vec_rhs.getData(HOST), HOST, DEVICE); + int status = mh_->matvec(&A, &vec_test, &vec_tmp, &ONE, &MINUSONE, DEVICE); + + if (status != 0) { + std::cout << "matvec failed" << std::endl; + std::exit( status ); + } + + true_residual_norm_ = sqrt(vh_->dot(&vec_tmp, &vec_tmp, DEVICE)); + + return true_residual_norm_; + } + + vector::Vector FunctionalityTestHelper::generateResidualVector(matrix::Csr& A, + vector::Vector& vec_x, + vector::Vector& vec_rhs) + { + using namespace memory; + using namespace ReSolve::constants; + + index_type n = A.getNumRows(); + + vector::Vector vec_r(n); + + vec_r.update(vec_rhs.getData(HOST), HOST, DEVICE); + + mh_->setValuesChanged(true, DEVICE); + + int status = mh_->matvec(&A, &vec_x, &vec_r, &ONE, &MINUSONE, DEVICE); + + if (status != 0) { + std::cout << "matvec from matrixhandler failed" << std::endl; + + std::exit(status); + } + + return vec_r; + } + + void FunctionalityTestHelper::printNorms(std::string &testname) + { + std::cout << "Results for " << testname << ":\n\n"; + + std::cout << std::scientific << std::setprecision(16); + + std::cout << "\t ||b-A*x||_2 : " + << residual_norm_ + << " (residual norm)\n"; + + std::cout << "\t ||b-A*x||_2/||b||_2 : " + << residual_norm_/rhs_norm_ + << " (relative residual norm)\n"; + + std::cout << "\t ||x-x_true||_2 : " + << diff_norm_ + << " (solution error)\n"; + + std::cout << "\t ||x-x_true||_2/||x_true||_2 : " + << diff_norm_/norm_solution_vector_ + << " (relative solution error)\n"; + + std::cout << "\t ||b-A*x_exact||_2 : " + << true_residual_norm_ + << " (control; residual norm with exact solution)\n"; + } + + int FunctionalityTestHelper::checkResidualNorm() + { + int error_sum = 0; + + if (!std::isfinite(residual_norm_/rhs_norm_)) { + std::cout << "Result is not a finite number!\n"; + error_sum++; + } + + if (residual_norm_/rhs_norm_ > tol_) { + std::cout << "Result inaccurate!\n"; + error_sum++; + } + + return error_sum; + } + + void FunctionalityTestHelper::calculateNorms( AxEqualsRhsProblem &problem ) + { + using namespace memory; + using namespace ReSolve::constants; + + matrix::Csr& A = *problem.getMatrix(); + vector::Vector& vec_x = *problem.getVector(); + vector::Vector& vec_rhs = *problem.getRhs(); + + int error_sum = 0; + + calculateSolutionVectorNorm(vec_x); + + // Compute residual norm + vector::Vector vec_r = generateResidualVector(A, vec_x, vec_rhs); + + calculateResidualNorm(vec_r); + + //for testing only - control + calculateRhsVectorNorm(vec_rhs); + + //compute ||x_diff|| = ||x - x_true|| norm + calculateDiffNorm(vec_x); + + calculateTrueNorm(A, vec_rhs); + } + + // note: this should be split into two functions for separate printing and checking + int FunctionalityTestHelper::checkResult(matrix::Csr& A, + vector::Vector& vec_rhs, + vector::Vector& vec_x, + SystemSolver& solver, + std::string testname) + { + using namespace memory; + using namespace ReSolve::constants; + + int error_sum = 0; + + error_sum += checkResidualNorm(); + + printNorms(testname); + + printIterativeSolverStats(solver); + + return error_sum; + } + + FunctionalityTestHelper::FunctionalityTestHelper(real_type tol_init, + workspace_type &workspace_init, + AxEqualsRhsProblem &problem) + : tol_(tol_init), + workspace_(workspace_init) + { + //mh_ = new MatrixHandler(&workspace_); + mh_ = std::make_unique(&workspace_); + + vh_ = std::make_unique(&workspace_); + + calculateNorms(problem); + } + + FunctionalityTestHelper::~FunctionalityTestHelper() + { + // no longer needed if mh_ is a unique_ptr, it will be deleted when out of scope + //delete mh_; + + //delete vh_; + } + + void FunctionalityTestHelper::printIterativeSolverStats(SystemSolver& solver) + { + // Get solver parameters + real_type tol = solver.getIterativeSolver().getTol(); + index_type restart = solver.getIterativeSolver().getRestart(); + index_type maxit = solver.getIterativeSolver().getMaxit(); + + // Get solver stats + index_type num_iter = solver.getIterativeSolver().getNumIter(); + real_type init_rnorm = solver.getIterativeSolver().getInitResidualNorm(); + real_type final_rnorm = solver.getIterativeSolver().getFinalResidualNorm(); + + std::cout << "\t IR iterations : " << num_iter << " (max " << maxit << ", restart " << restart << ")\n"; + std::cout << "\t IR starting res. norm : " << init_rnorm << "\n"; + std::cout << "\t IR final res. norm : " << final_rnorm << " (tol " << std::setprecision(2) << tol << ")\n\n"; + } + + int FunctionalityTestHelper::checkNormOfScaledResiduals(matrix::Csr& A, + vector::Vector& vec_rhs, + vector::Vector& vec_x, + SystemSolver& solver) + { + using namespace ReSolve::constants; + using namespace memory; + int error_sum = 0; + + // Compute residual norm for the second system + vector::Vector vec_r = generateResidualVector(A, vec_x, vec_rhs); + + // Compute norm of scaled residuals for the second system + real_type inf_norm_A = 0.0; + mh_->matrixInfNorm(&A, &inf_norm_A, DEVICE); + real_type inf_norm_x = vh_->infNorm(&vec_x, DEVICE); + real_type inf_norm_r = vh_->infNorm(&vec_r, DEVICE); + real_type nsr_norm = inf_norm_r / (inf_norm_A * inf_norm_x); + real_type nsr_system = solver.getNormOfScaledResiduals(&vec_rhs, &vec_x); + real_type error = std::abs(nsr_system - nsr_norm)/nsr_norm; + + if (error > 10.0*std::numeric_limits::epsilon()) { + std::cout << "Norm of scaled residuals computation failed:\n"; + std::cout << std::scientific << std::setprecision(16) + << "\tMatrix inf norm : " << inf_norm_A << "\n" + << "\tResidual inf norm : " << inf_norm_r << "\n" + << "\tSolution inf norm : " << inf_norm_x << "\n" + << "\tNorm of scaled residuals : " << nsr_norm << "\n" + << "\tNorm of scaled residuals (system): " << nsr_system << "\n\n"; + error_sum++; + } + + return error_sum; + } + + + // note: pass rel_residual_norm in as a function parameter rather than calculate here + int FunctionalityTestHelper::checkRelativeResidualNorm(vector::Vector& vec_rhs, + vector::Vector& vec_x, + SystemSolver& solver) + { + using namespace memory; + + int error_sum = 0; + + real_type rel_residual_norm = solver.getResidualNorm(&vec_rhs, &vec_x); + + real_type error = std::abs(rhs_norm_ * rel_residual_norm - residual_norm_)/residual_norm_; + + if (error > 10.0*std::numeric_limits::epsilon()) { + + std::cout << "Relative residual norm computation failed:\n"; + + std::cout << std::scientific << std::setprecision(16) + << "\tTest value : " << residual_norm_/rhs_norm_ << "\n" + << "\tSystemSolver computed : " << rel_residual_norm << "\n\n"; + + error_sum++; + } + + return error_sum; + } + + } //namespace tests +} //namespace ReSolve diff --git a/tests/functionality/FunctionalityTestHelper.hpp b/tests/functionality/FunctionalityTestHelper.hpp new file mode 100644 index 00000000..3ce8826d --- /dev/null +++ b/tests/functionality/FunctionalityTestHelper.hpp @@ -0,0 +1,94 @@ +namespace ReSolve +{ + namespace tests + { + /** + * @brief Class for compting and updating residual. + * + */ + class AxEqualsRhsProblem + { + public: + + AxEqualsRhsProblem(const std::string& matrix_filepath, const std::string& rhs_filepath); + ~AxEqualsRhsProblem(); + + ReSolve::matrix::Csr* getMatrix() const; + ReSolve::vector::Vector* getVector() const; + ReSolve::vector::Vector* getRhs() const; + + void updateProblem(const std::string& matrix_filepath, + const std::string& rhs_filepath); + + + private: + + ReSolve::matrix::Csr* A_; + ReSolve::vector::Vector* vec_x_; + ReSolve::vector::Vector* vec_rhs_; + }; + + /** + * @brief Class with utility functions for verifying test results. + * + */ + class FunctionalityTestHelper + { + public: + + FunctionalityTestHelper(ReSolve::real_type tol_init, + workspace_type &workspace_init, + AxEqualsRhsProblem &problem); + + ~FunctionalityTestHelper(); + + int checkResult(ReSolve::matrix::Csr& A, + ReSolve::vector::Vector& vec_rhs, + ReSolve::vector::Vector& vec_x, + ReSolve::SystemSolver& solver, + std::string testname); + + void printIterativeSolverStats(SystemSolver& solver); + + int checkNormOfScaledResiduals(ReSolve::matrix::Csr& A, + ReSolve::vector::Vector& vec_rhs, + ReSolve::vector::Vector& vec_x, + ReSolve::SystemSolver& solver); + + int checkRelativeResidualNorm(ReSolve::vector::Vector& vec_rhs, + ReSolve::vector::Vector& vec_x, + ReSolve::SystemSolver& solver); + + void calculateNorms(AxEqualsRhsProblem &problem); + + real_type calculateSolutionVectorNorm(ReSolve::vector::Vector& vec_x); + + real_type calculateRhsVectorNorm(ReSolve::vector::Vector& vec_x); + + real_type calculateResidualNorm(ReSolve::vector::Vector& vec_r); + + real_type calculateDiffNorm(ReSolve::vector::Vector& vec_x); + + int checkResidualNorm(); + + real_type calculateTrueNorm(ReSolve::matrix::Csr& A, ReSolve::vector::Vector& vec_rhs); + + void printNorms(std::string &testname); + + ReSolve::vector::Vector generateResidualVector(ReSolve::matrix::Csr& A, + ReSolve::vector::Vector& vec_x, + ReSolve::vector::Vector& vec_rhs); + + private: + workspace_type& workspace_; + std::unique_ptr mh_{nullptr}; + std::unique_ptr vh_{nullptr}; + ReSolve::real_type tol_{constants::DEFAULT_TOL}; + real_type residual_norm_{-1.0}; + real_type rhs_norm_{-1.0}; + real_type diff_norm_{-1.0}; + real_type norm_solution_vector_{-1.0}; + real_type true_residual_norm_{-1.0}; + }; + } // namespace tests +} // namespace ReSolve diff --git a/tests/functionality/testKLU.cpp b/tests/functionality/testKLU.cpp index 6b127be7..9860e93f 100644 --- a/tests/functionality/testKLU.cpp +++ b/tests/functionality/testKLU.cpp @@ -131,8 +131,9 @@ int main(int argc, char *argv[]) real_type normRmatrix1CPU = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::HOST)); std::cout<<"Results (first matrix): "<dot(vec_r, vec_r, ReSolve::memory::HOST)); std::cout<<"Results (second matrix): "<dot(vec_r, vec_r, ReSolve::memory::DEVICE)); std::cout<<"Results (first matrix): "<dot(vec_r, vec_r, ReSolve::memory::DEVICE)); std::cout<<"Results (second matrix): "< 10.0*std::numeric_limits::epsilon()) { - std::cout << "Relative residual norm computation failed:\n" << std::setprecision(16) + std::cout << "Relative residual norm computation failed:\n" + << std::scientific << std::setprecision(16) << "\tTest value : " << normRmatrix1/normB1 << "\n" << "\tSystemSolver computed : " << rel_residual_norm << "\n\n"; error_sum++; @@ -271,13 +272,15 @@ int main(int argc, char *argv[]) rel_residual_norm = solver.getResidualNorm(vec_rhs, vec_x); error = std::abs(normB2 * rel_residual_norm - normRmatrix2)/normRmatrix2; if (error > 10.0*std::numeric_limits::epsilon()) { - std::cout << "Relative residual norm computation failed:\n" << std::setprecision(16) + std::cout << "Relative residual norm computation failed:\n" + << std::scientific << std::setprecision(16) << "\tTest value : " << normRmatrix2/normB2 << "\n" << "\tSystemSolver computed : " << rel_residual_norm << "\n\n"; error_sum++; } std::cout << "Results (second matrix): " << std::endl << std::endl; + std::cout << std::scientific << std::setprecision(16); std::cout << "\t ||b-A*x||_2 : " << normRmatrix2 << " (residual norm)\n"; std::cout << "\t ||b-A*x||_2/||b||_2 : " << normRmatrix2/normB2 << " (relative residual norm)\n"; std::cout << "\t ||b-A*x||/(||A||*||x||) : " << nsr_norm << " (norm of scaled residuals)\n"; diff --git a/tests/functionality/testSysRefactor.cpp b/tests/functionality/testSysRefactor.cpp index 1ae563d7..89227072 100644 --- a/tests/functionality/testSysRefactor.cpp +++ b/tests/functionality/testSysRefactor.cpp @@ -36,309 +36,173 @@ std::string memory_space("cpu"); #endif +#include + + using namespace ReSolve::constants; +using namespace ReSolve::tests; using namespace ReSolve::colors; +using real_type = ReSolve::real_type; +using vector_type = ReSolve::vector::Vector; -int main(int argc, char *argv[]) +int verifyResult( AxEqualsRhsProblem &axb, + ReSolve::SystemSolver &solver, + workspace_type &workspace, + double tolerance ) { - // Use ReSolve data types. - using real_type = ReSolve::real_type; - using index_type = ReSolve::index_type; - using vector_type = ReSolve::vector::Vector; - - // Error sum needs to be 0 at the end for test to PASS. - // It is a FAIL otheriwse. int error_sum = 0; - int status = 0; + ReSolve::matrix::Csr* A = axb.getMatrix(); + vector_type* vec_x = axb.getVector(); + vector_type* vec_rhs = axb.getRhs(); - // Create workspace and initialize its handles. - workspace_type workspace; - workspace.initializeHandles(); + // larger tolerance than default 1e-17 because iterative refinement is not applied here + ReSolve::tests::FunctionalityTestHelper test_helper(tolerance, workspace, axb); - // Create linear algebra handlers - ReSolve::MatrixHandler matrix_handler(&workspace); - ReSolve::VectorHandler vector_handler(&workspace); + // Verify relative residual norm computation in SystemSolver + error_sum += test_helper.checkRelativeResidualNorm(*vec_rhs, *vec_x, solver); - // Create system solver - ReSolve::SystemSolver solver(&workspace); + // Compute norm of scaled residuals + error_sum += test_helper.checkNormOfScaledResiduals(*A, *vec_rhs, *vec_x, solver); - // Configure solver (CUDA-based solver needs slightly different - // settings than HIP-based one) - solver.setRefinementMethod("fgmres", "cgs2"); - solver.getIterativeSolver().setRestart(100); - if (memory_space == "hip") { - solver.getIterativeSolver().setMaxit(200); - } - if (memory_space == "cuda") { - solver.getIterativeSolver().setMaxit(400); - solver.getIterativeSolver().setTol(1e-17); - } - - // Input to this code is location of `data` directory where matrix files are stored - const std::string data_path = (argc == 2) ? argv[1] : "./"; + error_sum += + test_helper.checkResult(*A, *vec_rhs, *vec_x, solver, "first matrix"); + return error_sum; +} - std::string matrixFileName1 = data_path + "data/matrix_ACTIVSg2000_AC_00.mtx"; - std::string matrixFileName2 = data_path + "data/matrix_ACTIVSg2000_AC_02.mtx"; +void reportFinalResult(int const error_sum) +{ + std::cout << "Test KLU with Rf solver + IR "; + if (error_sum == 0) { + std::cout << GREEN << "PASSED" << CLEAR; + } else { + std::cout << RED << "FAILED" << CLEAR << ", error sum: "<< error_sum; + } + std::cout << std::endl << std::endl; +} - std::string rhsFileName1 = data_path + "data/rhs_ACTIVSg2000_AC_00.mtx.ones"; - std::string rhsFileName2 = data_path + "data/rhs_ACTIVSg2000_AC_02.mtx.ones"; +class FileInputs +{ + public: + + FileInputs(std::string const &data_path) + : + matrixFileName1_(data_path + "data/matrix_ACTIVSg2000_AC_00.mtx"), + matrixFileName2_(data_path + "data/matrix_ACTIVSg2000_AC_02.mtx"), + rhsFileName1_(data_path + "data/rhs_ACTIVSg2000_AC_00.mtx.ones"), + rhsFileName2_(data_path + "data/rhs_ACTIVSg2000_AC_02.mtx.ones") + { + } + std::string getMatrixFileName1() const + { + return matrixFileName1_; + } - // Read first matrix - std::ifstream mat1(matrixFileName1); - if(!mat1.is_open()) { - std::cout << "Failed to open file " << matrixFileName1 << "\n"; - return -1; + std::string getMatrixFileName2() const + { + return matrixFileName2_; } - ReSolve::matrix::Csr* A = ReSolve::io::createCsrFromFile(mat1, true); - A->syncData(ReSolve::memory::DEVICE); - mat1.close(); - - // Read first rhs vector - std::ifstream rhs1_file(rhsFileName1); - if(!rhs1_file.is_open()) { - std::cout << "Failed to open file " << rhsFileName1 << "\n"; - return -1; + + std::string getRhsFileName1() const + { + return rhsFileName1_; } - real_type* rhs = ReSolve::io::createArrayFromFile(rhs1_file); - rhs1_file.close(); - // Create rhs, solution and residual vectors - vector_type* vec_rhs = new vector_type(A->getNumRows()); - vector_type* vec_x = new vector_type(A->getNumRows()); - vector_type* vec_r = new vector_type(A->getNumRows()); + std::string getRhsFileName2() const + { + return rhsFileName2_; + } - // Allocate solution vector - vec_x->allocate(ReSolve::memory::HOST); //for KLU - vec_x->allocate(ReSolve::memory::DEVICE); + private: - // Set RHS vector on CPU (update function allocates) - vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); - vec_rhs->setDataUpdated(ReSolve::memory::HOST); + // Texas model + std::string const matrixFileName1_; + std::string const rhsFileName1_; - // Add system matrix to the solver - status = solver.setMatrix(A); - error_sum += status; + std::string const matrixFileName2_; + std::string const rhsFileName2_; +}; - // Solve the first system using KLU - status = solver.analyze(); - error_sum += status; +ReSolve::SystemSolver generateSolver(workspace_type *workspace, const AxEqualsRhsProblem &axb) +{ + ReSolve::SystemSolver solver(workspace); - status = solver.factorize(); - error_sum += status; + // Configure solver (CUDA-based solver needs slightly different + // settings than HIP-based one) + // cgs2 = classical Gram-Schmidt + solver.setRefinementMethod("fgmres", "cgs2"); - status = solver.solve(vec_rhs, vec_x); - error_sum += status; + solver.getIterativeSolver().setRestart(100); - // Evaluate the residual norm ||b-Ax|| on the device - vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - matrix_handler.setValuesChanged(true, ReSolve::memory::DEVICE); - status = matrix_handler.matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); - error_sum += status; - real_type normRmatrix1 = sqrt(vector_handler.dot(vec_r, vec_r, ReSolve::memory::DEVICE)); + if (memory_space == "hip") { + solver.getIterativeSolver().setMaxit(200); + } - // Compute norm of the rhs vector - real_type normB1 = sqrt(vector_handler.dot(vec_rhs, vec_rhs, ReSolve::memory::DEVICE)); + if (memory_space == "cuda") { + solver.getIterativeSolver().setMaxit(400); + solver.getIterativeSolver().setTol(1e-17); + } - // Compute solution vector norm - real_type normXtrue = sqrt(vector_handler.dot(vec_x, vec_x, ReSolve::memory::DEVICE)); + // next connect solver to problem: + int const status = solver.setMatrix(axb.getMatrix()); + + if(status != 0) { - // Compute norm of scaled residuals - real_type inf_norm_A = 0.0; - matrix_handler.matrixInfNorm(A, &inf_norm_A, ReSolve::memory::DEVICE); - real_type inf_norm_x = vector_handler.infNorm(vec_x, ReSolve::memory::DEVICE); - real_type inf_norm_r = vector_handler.infNorm(vec_r, ReSolve::memory::DEVICE); - real_type nsr_norm = inf_norm_r / (inf_norm_A * inf_norm_x); - real_type nsr_system = solver.getNormOfScaledResiduals(vec_rhs, vec_x); - real_type error = std::abs(nsr_system - nsr_norm)/nsr_norm; - - // Test norm of scaled residuals method in SystemSolver - if (error > 10.0*std::numeric_limits::epsilon()) { - std::cout << "Norm of scaled residuals computation failed:\n"; - std::cout << std::scientific << std::setprecision(16) - << "\tMatrix inf norm : " << inf_norm_A << "\n" - << "\tResidual inf norm : " << inf_norm_r << "\n" - << "\tSolution inf norm : " << inf_norm_x << "\n" - << "\tNorm of scaled residuals : " << nsr_norm << "\n" - << "\tNorm of scaled residuals (system): " << nsr_system << "\n\n"; - error_sum++; + std::cout << "solver.setMatrix(axb.getMatrix()) failed!" << std::endl; + std::exit(status); } - // Create reference vectors for testing purposes - vector_type* vec_test = new vector_type(A->getNumRows()); - vector_type* vec_diff = new vector_type(A->getNumRows()); + return solver; +} - // Set the reference solution vector (all ones) on both, CPU and GPU - real_type* x_data_ref = new real_type[A->getNumRows()]; - for (int i=0; igetNumRows(); ++i){ - x_data_ref[i] = 1.0; - } - vec_test->update(x_data_ref, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - vec_diff->update(x_data_ref, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - // compute ||x-x_true|| norm - vector_handler.axpy(&MINUSONE, vec_x, vec_diff, ReSolve::memory::DEVICE); - real_type normDiffMatrix1 = sqrt(vector_handler.dot(vec_diff, vec_diff, ReSolve::memory::DEVICE)); +int main(int argc, char *argv[]) +{ + // Input to this code is location of `data` directory where matrix files are stored + // build filenames for inputs + const std::string data_path = (argc == 2) ? argv[1] : "./"; - //compute the residual using exact solution - vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - status = matrix_handler.matvec(A, vec_test, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); - error_sum += status; - real_type exactSol_normRmatrix1 = sqrt(vector_handler.dot(vec_r, vec_r, ReSolve::memory::DEVICE)); + FileInputs const fileInputs( data_path ); - //evaluate the residual ON THE CPU using COMPUTED solution - vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); - status = matrix_handler.matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::HOST); - error_sum += status; - real_type normRmatrix1CPU = sqrt(vector_handler.dot(vec_r, vec_r, ReSolve::memory::HOST)); + // axb problem construction + AxEqualsRhsProblem axb(fileInputs.getMatrixFileName1(), + fileInputs.getRhsFileName1()); - // Verify relative residual norm computation in SystemSolver - real_type rel_residual_norm = solver.getResidualNorm(vec_rhs, vec_x); - error = std::abs(normB1 * rel_residual_norm - normRmatrix1)/normRmatrix1; - if (error > 10.0*std::numeric_limits::epsilon()) { - std::cout << "Relative residual norm computation failed:\n" << std::setprecision(16) - << "\tTest value : " << normRmatrix1/normB1 << "\n" - << "\tSystemSolver computed : " << rel_residual_norm << "\n\n"; - error_sum++; - } - - // Print out the result summary - std::cout << std::scientific << std::setprecision(16); - std::cout << "Results (first matrix): \n\n"; - std::cout << "\t ||b-A*x||_2 : " << normRmatrix1 << " (residual norm)" << std::endl; - std::cout << "\t ||b-A*x||_2 (CPU) : " << normRmatrix1CPU << " (residual norm)" << std::endl; - std::cout << "\t ||b-A*x||_2/||b||_2 : " << normRmatrix1/normB1 << " (relative residual norm)" << std::endl; - std::cout << "\t ||x-x_true||_2 : " << normDiffMatrix1 << " (solution error)" << std::endl; - std::cout << "\t ||x-x_true||_2/||x_true||_2 : " << normDiffMatrix1/normXtrue << " (relative solution error)" << std::endl; - std::cout << "\t ||b-A*x_exact||_2 : " << exactSol_normRmatrix1 << " (control; residual norm with exact solution)\n\n"; + // workspace construction + workspace_type workspace; + workspace.initializeHandles(); - // Now prepare the Rf solver - status = solver.refactorizationSetup(); - error_sum += status; - - // Load the second matrix - std::ifstream mat2(matrixFileName2); - if(!mat2.is_open()) { - std::cout << "Failed to open file " << matrixFileName2 << "\n"; - return -1; - } - ReSolve::io::updateMatrixFromFile(mat2, A); - A->syncData(ReSolve::memory::DEVICE); - mat2.close(); - - // Load the second rhs vector - std::ifstream rhs2_file(rhsFileName2); - if(!rhs2_file.is_open()) { - std::cout << "Failed to open file " << rhsFileName2 << "\n"; - return -1; - } - ReSolve::io::updateArrayFromFile(rhs2_file, &rhs); - rhs2_file.close(); + // solver construction + ReSolve::SystemSolver solver = generateSolver(&workspace, axb); - vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); + // Error sum must be 0 at the end for test to PASS. + int error_sum = 0; - status = solver.refactorize(); - error_sum += status; - - vec_x->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - status = solver.solve(vec_rhs, vec_x); - error_sum += status; - - // Compute residual norm for the second system - vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - matrix_handler.setValuesChanged(true, ReSolve::memory::DEVICE); - status = matrix_handler.matvec(A, vec_x, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); - error_sum += status; - real_type normRmatrix2 = sqrt(vector_handler.dot(vec_r, vec_r, ReSolve::memory::DEVICE)); - - //for testing only - control - real_type normB2 = sqrt(vector_handler.dot(vec_rhs, vec_rhs, ReSolve::memory::DEVICE)); - - // Compute norm of scaled residuals for the second system - inf_norm_A = 0.0; - matrix_handler.matrixInfNorm(A, &inf_norm_A, ReSolve::memory::DEVICE); - inf_norm_x = vector_handler.infNorm(vec_x, ReSolve::memory::DEVICE); - inf_norm_r = vector_handler.infNorm(vec_r, ReSolve::memory::DEVICE); - nsr_norm = inf_norm_r / (inf_norm_A * inf_norm_x); - nsr_system = solver.getNormOfScaledResiduals(vec_rhs, vec_x); - error = std::abs(nsr_system - nsr_norm)/nsr_norm; - - if (error > 10.0*std::numeric_limits::epsilon()) { - std::cout << "Norm of scaled residuals computation failed:\n"; - std::cout << std::scientific << std::setprecision(16) - << "\tMatrix inf norm : " << inf_norm_A << "\n" - << "\tResidual inf norm : " << inf_norm_r << "\n" - << "\tSolution inf norm : " << inf_norm_x << "\n" - << "\tNorm of scaled residuals : " << nsr_norm << "\n" - << "\tNorm of scaled residuals (system): " << nsr_system << "\n\n"; - error_sum++; - } + // Solve the first system + error_sum += solver.analyze(); - //compute ||x_diff|| = ||x - x_true|| norm - vec_diff->update(x_data_ref, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - vector_handler.axpy(&MINUSONE, vec_x, vec_diff, ReSolve::memory::DEVICE); - real_type normDiffMatrix2 = sqrt(vector_handler.dot(vec_diff, vec_diff, ReSolve::memory::DEVICE)); + error_sum += solver.factorize(); - //compute the residual using exact solution - vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE); - status = matrix_handler.matvec(A, vec_test, vec_r, &ONE, &MINUSONE, ReSolve::memory::DEVICE); - error_sum += status; - real_type exactSol_normRmatrix2 = sqrt(vector_handler.dot(vec_r, vec_r, ReSolve::memory::DEVICE)); + error_sum += solver.solve(axb.getRhs(), axb.getVector()); - // Verify relative residual norm computation in SystemSolver - rel_residual_norm = solver.getResidualNorm(vec_rhs, vec_x); - error = std::abs(normB2 * rel_residual_norm - normRmatrix2)/normRmatrix2; - if (error > 10.0*std::numeric_limits::epsilon()) { - std::cout << "Relative residual norm computation failed:\n" << std::setprecision(16) - << "\tTest value : " << normRmatrix2/normB2 << "\n" - << "\tSystemSolver computed : " << rel_residual_norm << "\n\n"; - error_sum++; - } - - // Get solver parameters - real_type tol = solver.getIterativeSolver().getTol(); - index_type restart = solver.getIterativeSolver().getRestart(); - index_type maxit = solver.getIterativeSolver().getMaxit(); - - // Get solver stats - index_type num_iter = solver.getIterativeSolver().getNumIter(); - real_type init_rnorm = solver.getIterativeSolver().getInitResidualNorm(); - real_type final_rnorm = solver.getIterativeSolver().getFinalResidualNorm(); + error_sum += verifyResult(axb, solver, workspace, 1e-12); + + // Now prepare the Rf solver + error_sum += solver.refactorizationSetup(); + + // update the Ax=b problem + axb.updateProblem(fileInputs.getMatrixFileName2(), fileInputs.getRhsFileName2()); + + error_sum += solver.refactorize(); + error_sum += solver.solve(axb.getRhs(), axb.getVector()); - std::cout << "Results (second matrix): " << std::endl << std::endl; - std::cout << "\t ||b-A*x||_2 : " << normRmatrix2 << " (residual norm)\n"; - std::cout << "\t ||b-A*x||_2/||b||_2 : " << normRmatrix2/normB2 << " (relative residual norm)\n"; - std::cout << "\t ||x-x_true||_2 : " << normDiffMatrix2 << " (solution error)\n"; - std::cout << "\t ||x-x_true||_2/||x_true||_2 : " << normDiffMatrix2/normXtrue << " (relative solution error)\n"; - std::cout << "\t ||b-A*x_exact||_2 : " << exactSol_normRmatrix2 << " (control; residual norm with exact solution)\n"; - std::cout << "\t IR iterations : " << num_iter << " (max " << maxit << ", restart " << restart << ")\n"; - std::cout << "\t IR starting res. norm : " << init_rnorm << "\n"; - std::cout << "\t IR final res. norm : " << final_rnorm << " (tol " << std::setprecision(2) << tol << ")\n\n"; - - if (!std::isfinite(normRmatrix1/normB1) || !std::isfinite(normRmatrix2/normB2)) { - std::cout << "Result is not a finite number!\n"; - error_sum++; - } - if ((normRmatrix1/normB1 > 1e-12 ) || (normRmatrix2/normB2 > 1e-15)) { - std::cout << "Result inaccurate!\n"; - error_sum++; - } - if (error_sum == 0) { - std::cout<<"Test KLU with Rf solver + IR " << GREEN << "PASSED" << CLEAR <