Refactor and stand alone examples for KLU class

examples/CMakeLists.txt

@@ -11,6 +11,11 @@ add_executable(rand_gmres.exe randGmres.cpp)
 target_link_libraries(rand_gmres.exe PRIVATE ReSolve)
 
 if(RESOLVE_USE_KLU)
+
+  # Build example with KLU factorization and KLU refactorization
+  add_executable(kluRefactor.exe kluRefactor.cpp)
+  target_link_libraries(kluRefactor.exe PRIVATE ReSolve)
+
   # Build example with KLU factorization and KLU refactorization
   add_executable(klu_klu.exe r_KLU_KLU.cpp)
   target_link_libraries(klu_klu.exe PRIVATE ReSolve)
@@ -19,6 +24,10 @@ if(RESOLVE_USE_KLU)
   add_executable(klu_klu_standalone.exe r_KLU_KLU_standalone.cpp)
   target_link_libraries(klu_klu_standalone.exe PRIVATE ReSolve)
 
+  # Read ONE matrix,solve with KLU
+  add_executable(kluStandAlone.exe kluStandAlone.cpp)
+  target_link_libraries(kluStandAlone.exe PRIVATE ReSolve)
+
   # Build a CPU example with a configurable system solver
   add_executable(system.exe r_SysSolver.cpp)
   target_link_libraries(system.exe PRIVATE ReSolve)
@@ -51,14 +60,14 @@ if(RESOLVE_USE_KLU)
     add_executable(klu_rf_fgmres_reuse_refactorization.exe r_KLU_rf_FGMRES_reuse_factorization.cpp)
     target_link_libraries(klu_rf_fgmres_reuse_refactorization.exe PRIVATE ReSolve)
 
-    # Build example with a configurable system solver
-    add_executable(system_cuda.exe r_SysSolverCuda.cpp)
-    target_link_libraries(system_cuda.exe PRIVATE ReSolve)
-
     # Build example where matrix data is updated
     add_executable(klu_glu_values_update.exe r_KLU_GLU_matrix_values_update.cpp)
     target_link_libraries(klu_glu_values_update.exe PRIVATE ReSolve)
 
+    # Build example with a configurable system solver
+    add_executable(system_cuda.exe r_SysSolverCuda.cpp)
+    target_link_libraries(system_cuda.exe PRIVATE ReSolve)
+
     # Example in which factorization is redone if solution is bad
     add_executable(klu_cusolverrf_check_redo.exe r_KLU_cusolverrf_redo_factorization.cpp)
     target_link_libraries(klu_cusolverrf_check_redo.exe PRIVATE ReSolve)
@@ -90,9 +99,12 @@ set(installable_executables "")
 list(APPEND installable_executables  rand_gmres.exe)
 
 if(RESOLVE_USE_KLU)
-  list(APPEND installable_executables klu_klu.exe
+  list(APPEND installable_executables kluRefactor.exe
+                                      klu_klu.exe
                                       klu_klu_standalone.exe
-                                      sysRefactor.exe)
+                                      sysRefactor.exe
+                                      kluStandAlone.exe
+                                      system.exe)
 
   if(RESOLVE_USE_GPU)
     list(APPEND installable_executables gpuRefactor.exe)

examples/kluRefactor.cpp

@@ -0,0 +1,181 @@
+#include <string>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <sstream>
+#include <resolve/matrix/Coo.hpp>
+#include <resolve/matrix/Csr.hpp>
+#include <resolve/matrix/Csc.hpp>
+#include <resolve/vector/Vector.hpp>
+#include <resolve/matrix/io.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 "ExampleHelper.hpp"
+#include <resolve/utilities/params/CliOptions.hpp>
+
+using namespace ReSolve::constants;
+
+/// Prints help message describing system usage.
+void printHelpInfo()
+{
+  std::cout << "\nLoads from files and solves a series of linear systems.\n\n";
+  std::cout << "System matrices are in files with names <pathname>XX.mtx, where XX are\n";
+  std::cout << "consecutive integer numbers 00, 01, 02, ...\n\n";
+  std::cout << "System right hand side vectors are stored in files with matching numbering.\n";
+  std::cout << "and file extension.\n\n";
+  std::cout << "Usage:\n\t./";
+  std::cout << "kluRefactor.exe -m <matrix pathname> -r <rhs pathname> -n <number of systems>\n\n";
+  std::cout << "Optional features:\n\t-h\tPrints this message.\n";
+  std::cout << "\t-i\tEnables iterative refinement.\n\n";
+}
+
+int main(int argc, char *argv[])
+{
+  // Use the same data types as those you specified in ReSolve build.
+  using index_type = ReSolve::index_type;
+  using real_type  = ReSolve::real_type;
+  using vector_type = ReSolve::vector::Vector;
+  using namespace ReSolve::examples;
+  using namespace ReSolve;
+  CliOptions options(argc, argv);
+  CliOptions::Option* opt = nullptr;
+
+  bool is_help = options.hasKey("-h");
+  if (is_help) {
+    printHelpInfo();
+    return 0;
+  }
+
+  bool is_iterative_refinement = options.hasKey("-i");
+
+  index_type num_systems = 0;
+  opt = options.getParamFromKey("-n");
+  if (opt) {
+    num_systems = std::stoi((opt->second).c_str());
+  } else {
+    std::cout << "Incorrect input!\n";
+    printHelpInfo();
+    return 1;
+  }
+
+  std::string matrix_path_name("");
+  opt = options.getParamFromKey("-m");
+  if (opt) {
+    matrix_path_name = opt->second;
+  } else {
+    std::cout << "Incorrect input!\n";
+    printHelpInfo();
+    return 1;
+  }
+
+  std::string rhs_path_name("");
+  opt = options.getParamFromKey("-r");
+  if (opt) {
+    rhs_path_name = opt->second;
+  } else {
+    std::cout << "Incorrect input!\n";
+    printHelpInfo();
+    return 1;
+  }
+
+  std::string fileId;
+  std::string rhsId;
+  std::string matrix_file_name_full;
+  std::string rhs_file_name_full;
+
+  matrix::Csr* A = nullptr;
+  LinAlgWorkspaceCpu workspace;
+  ExampleHelper<LinAlgWorkspaceCpu> helper(workspace);
+  MatrixHandler matrix_handler(&workspace);
+  VectorHandler vector_handler(&workspace);
+
+  vector_type* vec_rhs = nullptr;
+  vector_type* vec_x   = nullptr;
+
+  LinSolverDirectKLU* KLU = new LinSolverDirectKLU;
+  GramSchmidt GS(&vector_handler, GramSchmidt::CGS2);
+  LinSolverIterativeFGMRES FGMRES(&matrix_handler, &vector_handler, &GS);
+  for (int i = 0; i < num_systems; ++i)
+  {
+    std::cout << "System " << i << ":\n";
+
+    std::ostringstream matname;
+    std::ostringstream rhsname;
+    matname << matrix_path_name << std::setfill('0') << std::setw(2) << i << ".mtx";
+    rhsname << rhs_path_name    << std::setfill('0') << std::setw(2) << i << ".mtx";
+    matrix_file_name_full = matname.str();
+    rhs_file_name_full = rhsname.str();
+    std::ifstream mat_file(matrix_file_name_full);
+    if(!mat_file.is_open())
+    {
+      std::cout << "Failed to open file " << matrix_file_name_full << "\n";
+      return 1;
+    }
+    std::ifstream rhs_file(rhs_file_name_full);
+    if(!rhs_file.is_open())
+    {
+      std::cout << "Failed to open file " << rhs_file_name_full << "\n";
+      return 1;
+    }
+    bool is_expand_symmetric = true;
+    if (i == 0) {
+      A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric);
+
+      vec_rhs = ReSolve::io::createVectorFromFile(rhs_file);
+      vec_x = new vector_type(A->getNumRows());
+    } else {
+      ReSolve::io::updateMatrixFromFile(mat_file, A);
+      ReSolve::io::updateVectorFromFile(rhs_file, vec_rhs);
+    }
+    printSystemInfo(matrix_file_name_full, A);
+    mat_file.close();
+    rhs_file.close();
+
+    // Update data.
+    if (i < 2) {
+      vec_rhs->setDataUpdated(ReSolve::memory::HOST);
+    }
+    std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnz()<<std::endl;
+    //Now call direct solver
+    int status;
+
+    if (i < 2){
+      KLU->setup(A);
+      status = KLU->analyze();
+      std::cout<<"KLU analysis status: "<<status<<std::endl;
+      status = KLU->factorize();
+      std::cout<<"KLU factorization status: "<<status<<std::endl;
+    } else {
+      status =  KLU->refactorize();
+      std::cout<<"KLU re-factorization status: "<<status<<std::endl;
+    }
+    status = KLU->solve(vec_rhs, vec_x);
+    std::cout<<"KLU solve status: "<<status<<std::endl;
+
+    helper.resetSystem(A, vec_rhs, vec_x);
+    helper.printShortSummary();
+    if (is_iterative_refinement) {
+      // Setup iterative refinement
+      FGMRES.setup(A);
+      FGMRES.setupPreconditioner("LU", KLU);
+
+      // If refactorization produced finite solution do iterative refinement
+      if (std::isfinite(helper.getNormRelativeResidual())) {
+        FGMRES.solve(vec_rhs, vec_x);
+
+        // Print summary
+        helper.printIrSummary(&FGMRES);
+      }
+    }
+  }
+
+  //now DELETE
+  delete A;
+  delete KLU;
+  delete vec_rhs;
+  delete vec_x;
+  return 0;
+}

examples/kluStandAlone.cpp

@@ -0,0 +1,138 @@
+#include <string>
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+#include <resolve/matrix/Csr.hpp>
+#include <resolve/vector/Vector.hpp>
+#include <resolve/matrix/io.hpp>
+#include <resolve/matrix/MatrixHandler.hpp>
+#include <resolve/vector/VectorHandler.hpp>
+#include <resolve/LinSolverDirectKLU.hpp>
+#include <resolve/workspace/LinAlgWorkspace.hpp>
+#include <resolve/LinSolverIterativeFGMRES.hpp>
+#include "ExampleHelper.hpp"
+#include <resolve/utilities/params/CliOptions.hpp>
+
+using namespace ReSolve::constants;
+
+/// Prints help message describing system usage.
+void printHelpInfo()
+{
+  std::cout << "\nLoads from files and solves a linear system.\n\n";
+  std::cout << "Usage:\n\t./kluStandAlone.exe -m <matrix pathname> -r <rhs pathname>\n\n";
+  std::cout << "Optional features:\n\t-h\tPrints this message.\n";
+  std::cout << "\t-i\tEnables iterative refinement.\n\n";
+}
+
+int main(int argc, char *argv[])
+{
+  using index_type = ReSolve::index_type;
+  using real_type  = ReSolve::real_type;
+  using vector_type = ReSolve::vector::Vector;
+  using namespace ReSolve;
+  using namespace ReSolve::examples;
+
+  CliOptions options(argc, argv);
+  CliOptions::Option* opt = nullptr;
+
+  bool is_help = options.hasKey("-h");
+  if (is_help) {
+    printHelpInfo();
+    return 0;
+  }
+
+  bool is_iterative_refinement = options.hasKey("-i");
+
+  std::string matrix_path_name("");
+  opt = options.getParamFromKey("-m");
+  if (opt) {
+    matrix_path_name = opt->second;
+  } else {
+    std::cout << "Incorrect input!\n";
+    printHelpInfo();
+    return 1;
+  }
+
+  std::string rhs_path_name("");
+  opt = options.getParamFromKey("-r");
+  if (opt) {
+    rhs_path_name = opt->second;
+  } else {
+    std::cout << "Incorrect input!\n";
+    printHelpInfo();
+    return 1;
+  }
+
+  matrix::Csr* A = nullptr;
+  LinAlgWorkspaceCpu workspace;
+  ExampleHelper<LinAlgWorkspaceCpu> helper(workspace);
+  MatrixHandler matrix_handler(&workspace);
+  VectorHandler vector_handler(&workspace);
+  vector_type* vec_rhs = nullptr;
+  vector_type* vec_x = nullptr;
+
+  LinSolverDirectKLU* KLU = new LinSolverDirectKLU;
+
+  // Load the system
+  std::cout << "Solving the system:\n";
+
+  std::ifstream mat_file(matrix_path_name);
+  if(!mat_file.is_open())
+  {
+    std::cout << "Failed to open file " << matrix_path_name << "\n";
+    return 1;
+  }
+
+  std::ifstream rhs_file(rhs_path_name);
+  if(!rhs_file.is_open())
+  {
+    std::cout << "Failed to open file " << rhs_path_name << "\n";
+    return 1;
+  }
+
+  bool is_expand_symmetric = true;
+  A = ReSolve::io::createCsrFromFile(mat_file, is_expand_symmetric);
+  vec_rhs = ReSolve::io::createVectorFromFile(rhs_file);
+  vec_x = new vector_type(A->getNumRows());
+
+  printSystemInfo(matrix_path_name, A);
+  mat_file.close();
+  rhs_file.close();
+
+  vec_rhs->setDataUpdated(ReSolve::memory::HOST);
+  std::cout << "COO to CSR completed. Expanded NNZ: " << A->getNnz() << std::endl;
+
+  // Direct solver
+  int status;
+  KLU->setup(A);
+  status = KLU->analyze();
+  std::cout << "KLU analysis status: " << status << std::endl;
+  status = KLU->factorize();
+  std::cout << "KLU factorization status: " << status << std::endl;
+
+  // Solve the system
+  status = KLU->solve(vec_rhs, vec_x);
+  std::cout << "KLU solve status: " << status << std::endl;
+
+  helper.resetSystem(A, vec_rhs, vec_x);
+  helper.printShortSummary();
+  if (is_iterative_refinement) {
+    GramSchmidt GS(&vector_handler, GramSchmidt::CGS2);
+    LinSolverIterativeFGMRES FGMRES(&matrix_handler, &vector_handler, &GS);
+    // Setup iterative refinement
+    FGMRES.setup(A);
+    FGMRES.setupPreconditioner("LU", KLU);
+    // If refactorization produced finite solution do iterative refinement
+    if (std::isfinite(helper.getNormRelativeResidual())) {
+      FGMRES.solve(vec_rhs, vec_x);
+      helper.printIrSummary(&FGMRES);
+    }
+  }
+
+  // Cleanup
+  delete A;
+  delete KLU;
+  delete vec_rhs;
+  delete vec_x;
+  return 0;
+}