-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Harry Hughes
committed
Oct 16, 2024
1 parent
dfb4d39
commit c63d9ef
Showing
4 changed files
with
256 additions
and
33 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,196 @@ | ||
| #include "FunctionalityTestHelper.hpp" | ||
| #include <resolve/vector/Vector.hpp> | ||
| #include <resolve/matrix/io.hpp> | ||
| #include <resolve/matrix/Coo.hpp> | ||
| #include <resolve/matrix/Csr.hpp> | ||
| #include <resolve/matrix/Csc.hpp> | ||
| #include <resolve/matrix/MatrixHandler.hpp> | ||
| #include <resolve/vector/VectorHandler.hpp> | ||
| #include <resolve/LinSolverDirectKLU.hpp> | ||
| #include <resolve/LinSolverIterativeFGMRES.hpp> | ||
| #include <resolve/workspace/LinAlgWorkspace.hpp> | ||
| #include <resolve/SystemSolver.hpp> | ||
| #include <iostream> | ||
| #include <iomanip> | ||
|
|
||
| namespace ReSolve | ||
| { | ||
| namespace tests | ||
| { | ||
|
|
||
|
|
||
| // move definitions into here first, setup a PR, merge it in | ||
| FunctionalityTestHelper::FunctionalityTestHelper( ReSolve::real_type tol_init ) | ||
| : | ||
| tol_( tol_init ) | ||
| { | ||
| workspace_.initializeHandles(); | ||
| mh_ = new ReSolve::MatrixHandler(&workspace_); | ||
| vh_ = new ReSolve::VectorHandler(&workspace_); | ||
| } | ||
|
|
||
| int FunctionalityTestHelper::checkRefactorizationResult( | ||
| ReSolve::matrix::Csr& A, | ||
| ReSolve::vector::Vector& vec_rhs, | ||
| ReSolve::vector::Vector& vec_x, | ||
| ReSolve::SystemSolver& solver, | ||
| std::string testname) | ||
| { | ||
| using namespace memory; | ||
| using namespace ReSolve::constants; | ||
|
|
||
| int status = 0; | ||
| int error_sum = 0; | ||
|
|
||
| index_type n = A.getNumRows(); | ||
|
|
||
| true_norm_ = sqrt(vh_->dot(&vec_x, &vec_x, DEVICE)); | ||
|
|
||
| // Allocate vectors | ||
| ReSolve::vector::Vector vec_r(n); | ||
| ReSolve::vector::Vector vec_diff(n); | ||
| ReSolve::vector::Vector vec_test(n); | ||
|
|
||
| // Compute residual norm for the second system | ||
| vec_r.update(vec_rhs.getData(HOST), HOST, DEVICE); | ||
| mh_->setValuesChanged(true, DEVICE); | ||
| status = mh_->matvec(&A, &vec_x, &vec_r, &ONE, &MINUSONE, DEVICE); | ||
| error_sum += status; | ||
| residual_norm_ = sqrt(vh_->dot(&vec_r, &vec_r, DEVICE)); | ||
|
|
||
| //for testing only - control | ||
| rhs_norm_ = sqrt(vh_->dot(&vec_rhs, &vec_rhs, DEVICE)); | ||
|
|
||
| // Compute norm of scaled residuals: | ||
| // NSR = ||r||_inf / (||A||_inf * ||x||_inf) | ||
| error_sum += checkNormOfScaledResiduals(A, vec_rhs, vec_x, vec_r, solver); | ||
|
|
||
| //compute ||x_diff|| = ||x - x_true|| norm | ||
| vec_diff.setToConst(1.0, DEVICE); | ||
| vh_->axpy(&MINUSONE, &vec_x, &vec_diff, DEVICE); | ||
| diff_norm_ = sqrt(vh_->dot(&vec_diff, &vec_diff, DEVICE)); | ||
|
|
||
| //compute the residual using exact solution | ||
| vec_test.setToConst(1.0, DEVICE); | ||
| vec_r.update(vec_rhs.getData(HOST), HOST, DEVICE); | ||
| status = mh_->matvec(&A, &vec_test, &vec_r, &ONE, &MINUSONE, DEVICE); | ||
| error_sum += status; | ||
| true_residual_norm_ = sqrt(vh_->dot(&vec_r, &vec_r, DEVICE)); | ||
|
|
||
| // Verify relative residual norm computation in SystemSolver | ||
| error_sum += checkRelativeResidualNorm(vec_rhs, vec_x, residual_norm_, rhs_norm_, solver); | ||
|
|
||
| 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_/true_norm_ << " (relative solution error)\n"; | ||
| std::cout << "\t ||b-A*x_exact||_2 : " << true_residual_norm_ << " (control; residual norm with exact solution)\n"; | ||
| printIterativeSolverStats(solver); | ||
|
|
||
| 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::printIterativeSolverStats(SystemSolver& solver) | ||
| { | ||
| // Get solver parameters | ||
| real_type tol = solver.getIterativeSolver().getTol(); | ||
| 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(); | ||
| 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(ReSolve::matrix::Csr& A, | ||
| ReSolve::vector::Vector& vec_rhs, | ||
| ReSolve::vector::Vector& vec_x, | ||
| ReSolve::vector::Vector& vec_r, | ||
| ReSolve::SystemSolver& solver) | ||
| { | ||
| using namespace ReSolve::constants; | ||
| using namespace memory; | ||
| int error_sum = 0; | ||
|
|
||
| // 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<real_type>::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; | ||
| } | ||
|
|
||
| int FunctionalityTestHelper::checkRelativeResidualNorm(ReSolve::vector::Vector& vec_rhs, | ||
| ReSolve::vector::Vector& vec_x, | ||
| const real_type residual_norm, | ||
| const real_type rhs_norm, | ||
| ReSolve::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<real_type>::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; | ||
| } | ||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
||
| } | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters