diff --git a/tests/functionality/FunctionalityTestHelper.hpp b/tests/functionality/FunctionalityTestHelper.hpp index cf60d37f..db60ff97 100644 --- a/tests/functionality/FunctionalityTestHelper.hpp +++ b/tests/functionality/FunctionalityTestHelper.hpp @@ -27,7 +27,9 @@ namespace ReSolve { class FunctionalityTestHelper { public: - FunctionalityTestHelper() + FunctionalityTestHelper( ReSolve::real_type tol_init = constants::DEFAULT_TOL ) + : + tol_( tol_init ) { workspace_.initializeHandles(); mh_ = new ReSolve::MatrixHandler(&workspace_); @@ -116,6 +118,8 @@ namespace ReSolve { index_type restart = solver.getIterativeSolver().getRestart(); index_type maxit = solver.getIterativeSolver().getMaxit(); + // note: these are the solver's tolerance, different from the testhelper's tolerance + // Get solver stats index_type num_iter = solver.getIterativeSolver().getNumIter(); real_type init_rnorm = solver.getIterativeSolver().getInitResidualNorm(); diff --git a/tests/functionality/testSysRefactor.cpp b/tests/functionality/testSysRefactor.cpp index 63bccfc7..76ddd952 100644 --- a/tests/functionality/testSysRefactor.cpp +++ b/tests/functionality/testSysRefactor.cpp @@ -48,7 +48,6 @@ int main(int argc, char *argv[]) { // 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. @@ -59,27 +58,20 @@ int main(int argc, char *argv[]) // Create workspace and initialize its handles. workspace_type workspace; - // note: nvidia/amd gpu init workspace.initializeHandles(); // Create linear algebra handlers - // note: numerical linalg ops ReSolve::MatrixHandler matrix_handler(&workspace); ReSolve::VectorHandler vector_handler(&workspace); // Create system solver - // :) - // user facing class ReSolve::SystemSolver solver(&workspace); - // solver multiplexes according to preprocessors above, which come from cmake // Configure solver (CUDA-based solver needs slightly different // settings than HIP-based one) // cgs2 = classical Gram-Schmidt solver.setRefinementMethod("fgmres", "cgs2"); - // builds Krylov subspace to project solution onto: - // cool! solver.getIterativeSolver().setRestart(100); if (memory_space == "hip") { solver.getIterativeSolver().setMaxit(200); @@ -102,8 +94,6 @@ int main(int argc, char *argv[]) std::string rhsFileName1 = data_path + "data/rhs_ACTIVSg2000_AC_00.mtx.ones"; std::string rhsFileName2 = data_path + "data/rhs_ACTIVSg2000_AC_02.mtx.ones"; - // note: don't forget - this is RHS, not ones!! - // Read first matrix std::ifstream mat1(matrixFileName1); @@ -129,7 +119,6 @@ int main(int argc, char *argv[]) // 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()); // Allocate solution vector vec_x->allocate(ReSolve::memory::HOST); //for KLU @@ -138,13 +127,6 @@ int main(int argc, char *argv[]) // Set RHS vector on CPU (update function allocates) vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST); - // note: below is deleted and everything still works! - // this is not needed in this case because we went through the API - // this is a manual update for if you secretly change the underlying memory value unknown - // to the vector class () - - //vec_rhs->setDataUpdated(ReSolve::memory::HOST); - // Add system matrix to the solver status = solver.setMatrix(A); error_sum += status; @@ -159,107 +141,16 @@ int main(int argc, char *argv[]) status = solver.solve(vec_rhs, vec_x); error_sum += status; - ReSolve::tests::FunctionalityTestHelper testhelper; - // error_sum += testhelper.checkRefactorizationResult(*A, *vec_rhs, *vec_x, solver, "first matrix"); - - // 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)); - - // Compute norm of the rhs vector - real_type normB1 = sqrt(vector_handler.dot(vec_rhs, vec_rhs, ReSolve::memory::DEVICE)); - - // Compute solution vector norm - real_type normXtrue = sqrt(vector_handler.dot(vec_x, vec_x, ReSolve::memory::DEVICE)); - - // 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; - - - // note: these norms above are testing that proper calculation happens in solver - // question: why is testing the norm functionality part of this test? Could it be a separate - // test focused only on norm calculation? <--- in the future look into this - - // 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++; - } - - // Create reference vectors for testing purposes - vector_type* vec_test = new vector_type(A->getNumRows()); - vector_type* vec_diff = new vector_type(A->getNumRows()); - - // 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)); - - //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)); + ReSolve::tests::FunctionalityTestHelper testhelper( 1e-12 ); - //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)); - - // 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::cout << std::scientific << std::setprecision(16) - << "\tTest value : " << normRmatrix1/normB1 << "\n" - << "\tSystemSolver computed : " << rel_residual_norm << "\n\n"; - error_sum++; - } - - // Print out the result summary - // track these - 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"; - - - // second section begins - needs the first section so that we can "refactorize" - // can't refactorize from scratch, so above section is necessary in this test + error_sum += + testhelper.checkRefactorizationResult(*A, *vec_rhs, *vec_x, solver, "first matrix"); // Now prepare the Rf solver status = solver.refactorizationSetup(); error_sum += status; - // note: this tests a different I/O setup + // note: this tests a different I/O setup than the first section above // Load the second matrix std::ifstream mat2(matrixFileName2); @@ -289,17 +180,9 @@ int main(int argc, char *argv[]) status = solver.solve(vec_rhs, vec_x); error_sum += status; - // does what is done manually in the first part of the code - error_sum += testhelper.checkRefactorizationResult(*A, *vec_rhs, *vec_x, solver, "second matrix"); + error_sum += + testhelper.checkRefactorizationResult(*A, *vec_rhs, *vec_x, solver, "second matrix"); - 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 <