Adding asymmetric tests (from the data folder) for KLU, CuSolver, and RocSolver.

resolve/LinSolverIterativeFGMRES.cpp

@@ -138,6 +138,13 @@ namespace ReSolve
     io::Logger::misc() << "it 0: norm of residual "
                        << std::scientific << std::setprecision(16)
                        << rnorm << " Norm of rhs: " << bnorm << "\n";
+    if (rnorm / bnorm < MACHINE_EPSILON)
+    {
+      std::cout << "Iterative solve skipped. Initial residual is accurate to machine precision.\n"
+                << std::scientific << std::setprecision(16) << rnorm
+                << " / " << bnorm << "\n";
+      return 0;
+    }
     initial_residual_norm_ = rnorm;
     while (outer_flag)
     {

resolve/LinSolverIterativeRandFGMRES.cpp

@@ -180,6 +180,13 @@ namespace ReSolve
     io::Logger::misc() << "it 0: norm of residual "
                        << std::scientific << std::setprecision(16)
                        << rnorm << " Norm of rhs: " << bnorm << "\n";
+    if (rnorm / bnorm < MACHINE_EPSILON)
+    {
+      std::cout << "Iterative solve skipped. Initial residual is accurate to machine precision.\n"
+                << std::scientific << std::setprecision(16) << rnorm
+                << " / " << bnorm << "\n";
+      return 0;
+    }
     initial_residual_norm_ = rnorm;
     while (outer_flag)
     {

resolve/SystemSolver.cpp

@@ -13,7 +13,6 @@
 #ifdef RESOLVE_USE_KLU
 #include <resolve/LinSolverDirectKLU.hpp>
 #endif
-
 #include <resolve/LinSolverIterativeRandFGMRES.hpp>
 
 #ifdef RESOLVE_USE_CUDA
@@ -563,7 +562,8 @@ namespace ReSolve
   {
     int status = 0;
 
-    status += iterativeSolver_->resetMatrix(A_);
+    // status += iterativeSolver_->resetMatrix(A_);
+    std::cout << "Refining solution with iterative refinement ...\n";
     status += iterativeSolver_->solve(rhs, x);
 
     return status;

resolve/lusol/README.md

@@ -0,0 +1,2 @@
+# LUSOL disclaimer
+LUSOL code is experimental and not actively maintained or supported.

tests/functionality/CMakeLists.txt

@@ -36,9 +36,6 @@ if(RESOLVE_USE_KLU)
     add_executable(sys_glu_test.exe testSysGLU.cpp)
     target_link_libraries(sys_glu_test.exe PRIVATE ReSolve)
 
-    add_executable(sys_glu_asym_test.exe testSysGLUAsym.cpp)
-    target_link_libraries(sys_glu_asym_test.exe PRIVATE ReSolve)
-
   endif(RESOLVE_USE_CUDA)
 endif(RESOLVE_USE_KLU)
 
@@ -62,7 +59,6 @@ if(RESOLVE_USE_KLU)
 
   if(RESOLVE_USE_CUDA)
     list(APPEND installable_tests sys_glu_test.exe)
-    list(APPEND installable_tests sys_glu_asym_test.exe)
   endif(RESOLVE_USE_CUDA)
 endif()
 
@@ -76,7 +72,12 @@ install(TARGETS ${installable_tests}
 # Install directory with data files
 install(DIRECTORY data DESTINATION bin/resolve/tests/functionality)
 
-set(test_data_dir ${CMAKE_SOURCE_DIR}/tests/functionality/)
+set(test_data_dir ${CMAKE_SOURCE_DIR}/tests/functionality)
+set(sym_matrix_path /data/matrix_ACTIVSg2000_AC_)
+set(sym_rhs_path /data/rhs_ACTIVSg2000_AC_ones)
+set(asym_matrix_path /data/ScaleMicrogrid_Jacobian_N1000_)
+set(asym_rhs_path /data/ScaleMicrogrid_Residual_N1000_ones_number)
+
 
 add_test(NAME version COMMAND $<TARGET_FILE:version.exe>)
 
@@ -117,17 +118,20 @@ add_test(NAME rand_gmres_test    COMMAND $<TARGET_FILE:rand_gmres_test.exe>)
 
 if(RESOLVE_USE_KLU)
   # Using KLU as refactorization solver
-  add_test(NAME klu_klu_test    COMMAND $<TARGET_FILE:klu_klu_test.exe> "-d" "${test_data_dir}")
-  add_test(NAME klu_klu_ir_test COMMAND $<TARGET_FILE:klu_klu_test.exe> "-d" "${test_data_dir}" "-i")
+  add_test(NAME klu_klu_test    COMMAND $<TARGET_FILE:klu_klu_test.exe> "-d" "${test_data_dir}" "-m" "${sym_matrix_path}" "-r" "${sym_rhs_path}") # expected fail, tolerance too strict
+  add_test(NAME klu_klu_ir_test COMMAND $<TARGET_FILE:klu_klu_test.exe> "-d" "${test_data_dir}" "-m" "${sym_matrix_path}" "-r" "${sym_rhs_path}" "-i")
+  add_test(NAME klu_klu_asym_test COMMAND $<TARGET_FILE:klu_klu_test.exe> "-d" "${test_data_dir}" "-m" "${asym_matrix_path}" "-r" "${asym_rhs_path}")
+  add_test(NAME klu_klu_asym_ir_test COMMAND $<TARGET_FILE:klu_klu_test.exe> "-d" "${test_data_dir}" "-m" "${asym_matrix_path}" "-r" "${asym_rhs_path}" "-i") # expected fail, IR should not be used when tolerance is reached
 
   # CUDA-SDK specific tests
   if(RESOLVE_USE_CUDA)
     add_test(NAME klu_rf_cuda_test    COMMAND $<TARGET_FILE:klu_rf_test.exe> "-d" "${test_data_dir}")
     add_test(NAME klu_rf_ir_cuda_test COMMAND $<TARGET_FILE:klu_rf_test.exe> "-d" "${test_data_dir}" "-i")
     add_test(NAME klu_glu_cuda_test   COMMAND $<TARGET_FILE:klu_rf_test.exe> "-d" "${test_data_dir}" "-s" "glu")
-    add_test(NAME sys_refactor_cuda_test COMMAND $<TARGET_FILE:sys_refactor_test.exe> "${test_data_dir}")
-    add_test(NAME sys_glu_test COMMAND $<TARGET_FILE:sys_glu_test.exe> "${test_data_dir}")
-    add_test(NAME sys_glu_asym_test COMMAND $<TARGET_FILE:sys_glu_asym_test.exe> "${test_data_dir}")
+    add_test(NAME sys_refactor_cuda_test COMMAND $<TARGET_FILE:sys_refactor_test.exe> "${test_data_dir}" "-d" "${test_data_dir}" "-m" "${sym_matrix_path}" "-r" "${sym_rhs_path}")
+    add_test(NAME sys_refactor_cuda_asym_test COMMAND $<TARGET_FILE:sys_refactor_test.exe> "${test_data_dir}" "-d" "${test_data_dir}" "-m" "${asym_matrix_path}" "-r" "${asym_rhs_path}") # expected fail, IR should not be used when tolerance is reached
+    add_test(NAME sys_glu_test COMMAND $<TARGET_FILE:sys_glu_test.exe> "${test_data_dir}" "-d" "${test_data_dir}" "-m" "${sym_matrix_path}" "-r" "${sym_rhs_path}")
+    add_test(NAME sys_glu_asym_test COMMAND $<TARGET_FILE:sys_glu_test.exe> "${test_data_dir}" "-d" "${test_data_dir}" "-m" "${asym_matrix_path}" "-r" "${asym_rhs_path}") 
   endif(RESOLVE_USE_CUDA)
 
   # ROCm specific tests
@@ -136,6 +140,7 @@ if(RESOLVE_USE_KLU)
     add_test(NAME rocsolver_rf_m1_test    COMMAND $<TARGET_FILE:klu_rf_test.exe> "-d" "${test_data_dir}" "-m" "rocsparse_trisolve")
     add_test(NAME rocsolver_rf_ir_test    COMMAND $<TARGET_FILE:klu_rf_test.exe> "-d" "${test_data_dir}" "-i")
     add_test(NAME rocsolver_rf_ir_m1_test COMMAND $<TARGET_FILE:klu_rf_test.exe> "-d" "${test_data_dir}" "-i" "-m" "rocsparse_trisolve")
-    add_test(NAME sys_refactor_hip_test  COMMAND $<TARGET_FILE:sys_refactor_test.exe>  "${test_data_dir}")
+    add_test(NAME sys_refactor_cuda_test COMMAND $<TARGET_FILE:sys_refactor_test.exe> "${test_data_dir}" "-d" "${test_data_dir}" "-m" "${sym_matrix_path}" "-r" "${sym_rhs_path}")
+    add_test(NAME sys_refactor_cuda_asym_test COMMAND $<TARGET_FILE:sys_refactor_test.exe> "${test_data_dir}" "-d" "${test_data_dir}" "-m" "${asym_matrix_path}" "-r" "${asym_rhs_path}") # expected fail, IR should not be used when tolerance is reached
   endif(RESOLVE_USE_HIP)
 endif(RESOLVE_USE_KLU)

tests/functionality/TestHelper.hpp

@@ -232,9 +232,6 @@ class TestHelper
       std::cout << "\t Residual norm on CPU    ||b-A*x|| (CPU)         : " << getNormResidualCpu() << "\n";
     }
     std::cout << "\t Relative residual norm  ||b-A*x||/||b||         : " << getNormResidualScaled() << "\n";
-    std::cout << "\t Error norm              ||x-x_true||            : " << getNormDiff() << "\n";
-    std::cout << "\t Relative error norm     ||x-x_true||/||x_true|| : " << getNormDiffScaled() << "\n";
-    std::cout << "\t Exact solution residual ||b-A*x_true||          : " << getNormResidualTrue() << "\n";
   }
 
   /// Summary of error norms for an iterative refinement test.
@@ -269,7 +266,6 @@ class TestHelper
   {
     int                error_sum = 0;
     ReSolve::real_type norm      = norm_res_ / norm_rhs_;
-
     if (!std::isfinite(norm))
     {
       std::cout << "Result is not a finite number!\n";

tests/functionality/testKlu.cpp

@@ -58,6 +58,30 @@ int runTest(int argc, char* argv[], std::string& solver_name)
   auto        opt       = options.getParamFromKey("-d");
   std::string data_path = opt ? (*opt).second : ".";
 
+  // Get matrix file name
+  opt = options.getParamFromKey("-m");
+  if (!opt)
+  {
+    std::cout << "Matrix file name not provided. Use -m <matrix_file_name>.\n";
+    return -1;
+  }
+  std::string matrix_temp = (*opt).second;
+
+  // Get rhs file name
+  opt = options.getParamFromKey("-r");
+  if (!opt)
+  {
+    std::cout << "RHS file name not provided. Use -r <rhs_file_name>.\n";
+    return -1;
+  }
+  std::string rhs_temp = (*opt).second;
+
+  // Construct matrix and rhs file names from inputs
+  std::string matrix_file_name_1 = data_path + matrix_temp + "01.mtx";
+  std::string matrix_file_name_2 = data_path + matrix_temp + "02.mtx";
+  std::string rhs_file_name_1    = data_path + rhs_temp + "01.mtx";
+  std::string rhs_file_name_2    = data_path + rhs_temp + "02.mtx";
+
   // Whether to use iterative refinement
   opt        = options.getParamFromKey("-i");
   bool is_ir = opt ? true : false;
@@ -78,13 +102,6 @@ int runTest(int argc, char* argv[], std::string& solver_name)
   ReSolve::GramSchmidt              GS(&vector_handler, ReSolve::GramSchmidt::CGS2);
   ReSolve::LinSolverIterativeFGMRES FGMRES(&matrix_handler, &vector_handler, &GS);
 
-  // Input data
-  std::string matrix_file_name_1 = data_path + "/data/matrix_ACTIVSg200_AC_10.mtx";
-  std::string matrix_file_name_2 = data_path + "/data/matrix_ACTIVSg200_AC_11.mtx";
-
-  std::string rhs_file_name_1 = data_path + "/data/rhs_ACTIVSg200_AC_10.mtx.ones";
-  std::string rhs_file_name_2 = data_path + "/data/rhs_ACTIVSg200_AC_11.mtx.ones";
-
   // Read first matrix
   std::ifstream mat1(matrix_file_name_1);
   if (!mat1.is_open())
@@ -125,8 +142,9 @@ int runTest(int argc, char* argv[], std::string& solver_name)
 
   status = KLU.solve(&vec_rhs, &vec_x);
   error_sum += status;
+  helper.setSystem(A, &vec_rhs, &vec_x);
 
-  if (is_ir)
+  if (is_ir && helper.checkResult(ReSolve::constants::MACHINE_EPSILON)) // only perform IR if you haven't reached machine accuracy
   {
     test_name += " + IR";
 
@@ -141,16 +159,16 @@ int runTest(int argc, char* argv[], std::string& solver_name)
   }
 
   // Compute error norms for the system
-  helper.setSystem(A, &vec_rhs, &vec_x);
+  helper.resetSystem(A, &vec_rhs, &vec_x);
 
   // Print result summary and check solution
   std::cout << "\nResults (first matrix): \n\n";
   helper.printSummary();
-  if (is_ir)
+  if (is_ir && helper.checkResult(ReSolve::constants::MACHINE_EPSILON)) // only perform IR if you haven't reached machine accuracy
   {
     helper.printIrSummary(&FGMRES);
   }
-  error_sum += helper.checkResult(ReSolve::constants::MACHINE_EPSILON);
+  error_sum += helper.checkResult(10 * ReSolve::constants::MACHINE_EPSILON); // tolerance increased to deal with system
 
   // Load the second matrix
   std::ifstream mat2(matrix_file_name_2);
@@ -177,8 +195,9 @@ int runTest(int argc, char* argv[], std::string& solver_name)
   error_sum += status;
   std::cout << "KLU refactorization status: " << status << std::endl;
 
-  if (is_ir)
+  if (is_ir && helper.checkResult(ReSolve::constants::MACHINE_EPSILON)) // only perform IR if you haven't reached machine accuracy
   {
+    std::cout << "Second matrix: Performing iterative refinement...\n";
     FGMRES.resetMatrix(A);
     status = FGMRES.setupPreconditioner("LU", &KLU);
     error_sum += status;
@@ -196,11 +215,11 @@ int runTest(int argc, char* argv[], std::string& solver_name)
 
   std::cout << "\nResults (second matrix): \n\n";
   helper.printSummary();
-  if (is_ir)
+  if (is_ir && helper.checkResult(ReSolve::constants::MACHINE_EPSILON)) // only perform IR if you haven't reached machine accuracy
   {
     helper.printIrSummary(&FGMRES);
   }
-  error_sum += helper.checkResult(ReSolve::constants::MACHINE_EPSILON);
+  error_sum += helper.checkResult(100 * ReSolve::constants::MACHINE_EPSILON); // tolerance increased to deal with system
 
   isTestPass(error_sum, test_name);
 

tests/functionality/testSysGLU.cpp

@@ -18,6 +18,7 @@
 #include <resolve/matrix/Csr.hpp>
 #include <resolve/matrix/MatrixHandler.hpp>
 #include <resolve/matrix/io.hpp>
+#include <resolve/utilities/params/CliOptions.hpp>
 #include <resolve/vector/Vector.hpp>
 #include <resolve/vector/VectorHandler.hpp>
 #include <resolve/workspace/LinAlgWorkspace.hpp>
@@ -42,38 +43,60 @@ int main(int argc, char* argv[])
   int error_sum = 0; ///< error sum zero means test passes, otherwise fails
   int status    = 0;
 
+  // Collect all command line options
+  ReSolve::CliOptions options(argc, argv);
+
+  // Get directory with input files
+  auto        opt       = options.getParamFromKey("-d");
+  std::string data_path = opt ? (*opt).second : ".";
+
+  // Get matrix file name
+  opt = options.getParamFromKey("-m");
+  if (!opt)
+  {
+    std::cout << "Matrix file name not provided. Use -m <matrix_file_name>.\n";
+    return -1;
+  }
+  std::string matrix_temp = (*opt).second;
+
+  // Get rhs file name
+  opt = options.getParamFromKey("-r");
+  if (!opt)
+  {
+    std::cout << "RHS file name not provided. Use -r <rhs_file_name>.\n";
+    return -1;
+  }
+  std::string rhs_temp = (*opt).second;
+
+  // Construct matrix and rhs file names from inputs
+  std::string matrix_file_name_1 = data_path + matrix_temp + "01.mtx";
+  std::string matrix_file_name_2 = data_path + matrix_temp + "02.mtx";
+  std::string rhs_file_name_1    = data_path + rhs_temp + "01.mtx";
+  std::string rhs_file_name_2    = data_path + rhs_temp + "02.mtx";
+
   workspace_type workspace;
   workspace.initializeHandles();
   ReSolve::MatrixHandler matrix_handler(&workspace);
   ReSolve::VectorHandler vector_handler(&workspace);
 
   ReSolve::SystemSolver solver(&workspace, "klu", "glu", "glu", "none", "none");
 
-  // Input to this code is location of `data` directory where matrix files are stored
-  const std::string data_path = (argc == 2) ? argv[1] : "./";
-
-  std::string matrixFileName1 = data_path + "data/matrix_ACTIVSg200_AC_10.mtx";
-  std::string matrixFileName2 = data_path + "data/matrix_ACTIVSg200_AC_11.mtx";
-
-  std::string rhsFileName1 = data_path + "data/rhs_ACTIVSg200_AC_10.mtx.ones";
-  std::string rhsFileName2 = data_path + "data/rhs_ACTIVSg200_AC_11.mtx.ones";
-
   // Read first matrix
-  std::ifstream mat1(matrixFileName1);
+  std::ifstream mat1(matrix_file_name_1);
   if (!mat1.is_open())
   {
-    std::cout << "Failed to open file " << matrixFileName1 << "\n";
+    std::cout << "Failed to open file " << matrix_file_name_1 << "\n";
     return -1;
   }
   ReSolve::matrix::Csr* A = ReSolve::io::createCsrFromFile(mat1);
   A->syncData(ReSolve::memory::DEVICE);
   mat1.close();
 
   // Read first rhs vector
-  std::ifstream rhs1_file(rhsFileName1);
+  std::ifstream rhs1_file(rhs_file_name_1);
   if (!rhs1_file.is_open())
   {
-    std::cout << "Failed to open file " << rhsFileName1 << "\n";
+    std::cout << "Failed to open file " << rhs_file_name_1 << "\n";
     return -1;
   }
   real_type*   rhs     = ReSolve::io::createArrayFromFile(rhs1_file);
@@ -202,21 +225,21 @@ int main(int argc, char* argv[])
   std::cout << "\t ||b-A*x_true||_2  (control) : " << exactSol_normRmatrix1 << " (residual norm with exact solution)\n\n";
 
   // Load the second matrix
-  std::ifstream mat2(matrixFileName2);
+  std::ifstream mat2(matrix_file_name_2);
   if (!mat2.is_open())
   {
-    std::cout << "Failed to open file " << matrixFileName2 << "\n";
+    std::cout << "Failed to open file " << matrix_file_name_2 << "\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);
+  std::ifstream rhs2_file(rhs_file_name_2);
   if (!rhs2_file.is_open())
   {
-    std::cout << "Failed to open file " << rhsFileName2 << "\n";
+    std::cout << "Failed to open file " << rhs_file_name_2 << "\n";
     return -1;
   }
   ReSolve::io::updateArrayFromFile(rhs2_file, &rhs);