Functionality tests utils dev

tests/functionality/CMakeLists.txt

@@ -6,6 +6,11 @@
 
 ]]
 
+# Captain!
+# 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)
@@ -36,9 +41,10 @@ if(RESOLVE_USE_CUDA)
     add_executable(klu_rf_fgmres_test.exe testKLU_Rf_FGMRES.cpp)
     target_link_libraries(klu_rf_fgmres_test.exe PRIVATE ReSolve)
 
+    # Captain!
     # 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 +75,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

tests/functionality/FunctionalityTestHelper.cpp

@@ -0,0 +1,197 @@
+#include <string>
+#include <iostream>
+#include <iomanip>
+
+#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>
+
+#if defined (RESOLVE_USE_CUDA)
+#include <resolve/LinSolverDirectCuSolverRf.hpp>
+  using workspace_type = ReSolve::LinAlgWorkspaceCUDA;
+  std::string memory_space("cuda");
+#elif defined (RESOLVE_USE_HIP)
+#include <resolve/LinSolverDirectRocSolverRf.hpp>
+  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{
+
+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;
+}
+
+FunctionalityTestHelper::FunctionalityTestHelper( 
+  ReSolve::real_type tol_init )
+  :
+  tol_( tol_init )
+{
+  workspace_.initializeHandles();
+  mh_ = new ReSolve::MatrixHandler(&workspace_);
+  vh_ = new ReSolve::VectorHandler(&workspace_);
+}
+
+FunctionalityTestHelper::~FunctionalityTestHelper()
+{
+  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();
+
+  // 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;
+}
+
+}
+}

tests/functionality/FunctionalityTestHelper.hpp

@@ -0,0 +1,45 @@
+namespace ReSolve 
+{
+  namespace tests
+  {
+    class FunctionalityTestHelper
+    {
+      public:
+
+        FunctionalityTestHelper( ReSolve::real_type tol_init = constants::DEFAULT_TOL );
+
+        ~FunctionalityTestHelper();
+
+        int checkRefactorizationResult(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::vector::Vector& vec_r,
+                                      ReSolve::SystemSolver& solver);
+
+        int checkRelativeResidualNorm(ReSolve::vector::Vector& vec_rhs,
+                                      ReSolve::vector::Vector& vec_x,
+                                      const real_type residual_norm,
+                                      const real_type rhs_norm,
+                                      ReSolve::SystemSolver& solver);
+
+      private:
+        workspace_type workspace_;
+        ReSolve::MatrixHandler* mh_{nullptr};
+        ReSolve::VectorHandler* 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 true_norm_{-1.0};
+        real_type true_residual_norm_{-1.0};
+    };
+  } // namespace tests
+} // namespace ReSolve

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): "<<std::endl<<std::endl;
-  std::cout<<"\t ||b-A*x||_2                 : " << std::setprecision(16) << normRmatrix1    << " (residual norm)" << std::endl;
-  std::cout<<"\t ||b-A*x||_2  (CPU)          : " << std::setprecision(16) << normRmatrix1CPU << " (residual norm)" << std::endl;
+  std::cout << std::scientific << std::setprecision(16);
+  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   << " (scaled 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 << " (scaled solution error)"        << std::endl;
@@ -190,6 +191,7 @@ int main(int argc, char *argv[])
   real_type exactSol_normRmatrix2 = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::HOST));
 
   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)"<<std::endl;
   std::cout<<"\t ||b-A*x||_2/||b||_2         : "<<normRmatrix2/normB2<<" (scaled residual norm)"<<std::endl;
   std::cout<<"\t ||x-x_true||_2              : "<<normDiffMatrix2<<" (solution error)"<<std::endl;

tests/functionality/testKLU_RocSolver.cpp

@@ -158,8 +158,9 @@ int main(int argc, char *argv[])
   real_type normRmatrix1CPU = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE));
 
   std::cout<<"Results (first matrix): "<<std::endl<<std::endl;
-  std::cout<<"\t ||b-A*x||_2                 : " << std::setprecision(16) << normRmatrix1    << " (residual norm)" << std::endl;
-  std::cout<<"\t ||b-A*x||_2  (CPU)          : " << std::setprecision(16) << normRmatrix1CPU << " (residual norm)" << std::endl;
+  std::cout << std::scientific << std::setprecision(16);
+  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   << " (scaled 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 << " (scaled solution error)"        << std::endl;
@@ -220,6 +221,7 @@ int main(int argc, char *argv[])
   real_type exactSol_normRmatrix2 = sqrt(vector_handler->dot(vec_r, vec_r, ReSolve::memory::DEVICE));
 
   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)"<<std::endl;
   std::cout<<"\t ||b-A*x||_2/||b||_2         : "<<normRmatrix2/normB2<<" (scaled residual norm)"<<std::endl;
   std::cout<<"\t ||x-x_true||_2              : "<<normDiffMatrix2<<" (solution error)"<<std::endl;
@@ -235,9 +237,9 @@ int main(int argc, char *argv[])
     error_sum++;
   }
   if (error_sum == 0) {
-    std::cout << "Test KLU with cuSolverGLU refactorization " << GREEN << "PASSED" << CLEAR << std::endl;
+    std::cout << "Test KLU with rocsolverRf refactorization " << GREEN << "PASSED" << CLEAR << std::endl;
   } else {
-    std::cout << "Test KLU with cuSolverGLU refactorization " << RED << "FAILED" << CLEAR
+    std::cout << "Test KLU with rocsolverRf refactorization " << RED << "FAILED" << CLEAR
               << ", error sum: " << error_sum << std::endl;
   }
 

tests/functionality/testLUSOL.cpp

@@ -142,7 +142,8 @@ int main(int argc, char* argv[])
   real_type exactSol_normRmatrix = sqrt(vector_handler.dot(&vec_r, &vec_r, ReSolve::memory::HOST));
 
   std::cout << "Results: \n";
-  std::cout << "\t ||b-A*x||_2                 : " << std::setprecision(16) << normRmatrix << " (residual norm)\n";
+  std::cout << std::scientific << std::setprecision(16);
+  std::cout << "\t ||b-A*x||_2                 : " << normRmatrix << " (residual norm)\n";
   std::cout << "\t ||b-A*x||_2/||b||_2         : " << normRmatrix / normB << " (scaled residual norm)\n";
   std::cout << "\t ||x-x_true||_2              : " << normDiffMatrix << " (solution error)\n";
   std::cout << "\t ||x-x_true||_2/||x_true||_2 : " << normDiffMatrix / normXtrue << " (scaled solution error)\n";
@@ -243,7 +244,8 @@ int main(int argc, char* argv[])
   exactSol_normRmatrix = sqrt(vector_handler.dot(&vec_r, &vec_r, ReSolve::memory::HOST));
 
   std::cout << "Results: \n";
-  std::cout << "\t ||b-A*x||_2                 : " << std::setprecision(16) << normRmatrix << " (residual norm)\n";
+  std::cout << std::scientific << std::setprecision(16);
+  std::cout << "\t ||b-A*x||_2                 : " << normRmatrix << " (residual norm)\n";
   std::cout << "\t ||b-A*x||_2/||b||_2         : " << normRmatrix / normB << " (scaled residual norm)\n";
   std::cout << "\t ||x-x_true||_2              : " << normDiffMatrix << " (solution error)\n";
   std::cout << "\t ||x-x_true||_2/||x_true||_2 : " << normDiffMatrix / normXtrue << " (scaled solution error)\n";

tests/functionality/testSysGLU.cpp

@@ -180,7 +180,8 @@ int main(int argc, char *argv[])
   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<real_type>::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<real_type>::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";