Expose ILU0 numeric boost configuration for CUDA and HIP

CHANGELOG.md

@@ -50,8 +50,21 @@
 - Fixed GPU refactorization allocation and synchronization errors, uninitialized GMRES example initial guesses, and CUDA ILU0 failures on stored numerical zero pivots.
 
 - Added cuDSS implementation of Cholesky solver (HyKKT) and refactorization
+
 - Optimized CGS2 coefficient accumulation in `GramSchmidt` by eliminating an unnecessary device-to-host synchronization.
 
+- Added examples/sysGmres.cpp to demonstrate how to use SystemSolver with GMRES.
+
+- Updated MatrixHandler::addConst to return integer error codes instead of void.
+
+- Added a preconditioner interface class so users can define their own preconditioners.
+
+- Added left preconditioning support for GMRES and a user-defined preconditioner class.
+
+- Added optional timing output to `gpuRefactor`, `kluRefactor`, `gluRefactor`, and `sysRefactor` plus benchmark utilities for parsing logs and generating timing/residual plots.
+
+- Added configurable zero-pivot and pivot-boost parameters for CUDA and HIP ILU0 and exposed the preconditioner solver through `SystemSolver` for configuring them.
+
 ## Changes to Re::Solve in release 0.99.2
 
 ### Major Features
@@ -105,13 +118,3 @@ It is seamless from the user perspective and fixed many bugs.
 13. Added LUSOL direct solver, which can factorize matrices and extract factors independently of KLU.
 
 14. Various Spack updates.
-
-15. Added examples/sysGmres.cpp to demonstrate how to use SystemSolver with GMRES. 
-
-16. Updated MatrixHandler::addConst to return integer error codes instead of void.
-
-17. Added a preconditioner interface class so users can define their own preconditioners.
-
-18. Added left preconditioning support for GMRES and a user-defined preconditioner class.
-
-19. Added optional timing output to `gpuRefactor`, `kluRefactor`, `gluRefactor`, and `sysRefactor` plus benchmark utilities for parsing logs and generating timing/residual plots.

examples/sysGmres.cpp

@@ -9,6 +9,7 @@
 #include <string>
 
 #include "ExampleHelper.hpp"
+#include <resolve/LinSolverDirect.hpp>
 #include <resolve/SystemSolver.hpp>
 #include <resolve/matrix/Csr.hpp>
 #include <resolve/matrix/io.hpp>
@@ -33,6 +34,10 @@ void printHelpInfo()
   std::cout << "\t-h \tPrints this message.\n";
   std::cout << "\t-i <iter method> \tIterative method: randgmres or fgmres (default 'randgmres').\n";
   std::cout << "\t-g <gs method> \tGram-Schmidt method: cgs1, cgs2, or mgs (default 'cgs2').\n";
+  std::cout << "\t-n <yes|no> \tEnable numeric boost on CUDA/HIP (default 'yes' on CUDA, 'no' on HIP).\n";
+  std::cout << "\t-t <boost tolerance> \tNumeric boost tolerance for CUDA/HIP (default '1e-6').\n";
+  std::cout << "\t-v <boost value> \tNumeric boost replacement value for CUDA/HIP (default '1e-6').\n";
+  std::cout << "\t-z <zero diagonal> \tZero diagonal replacement value for CPU (default '1e-6').\n";
   std::cout << "\t-s <sketching method> \tSketching method: count or fwht (default 'count')\n";
   std::cout << "\t-x <flexible> \tEnable flexible: yes or no (default 'yes')\n\n";
   std::cout << "\t-p <preconditioner side> \tPreconditioner side: left or right (default 'right')\n\n";
@@ -60,6 +65,11 @@ static void processInputs(std::string& method,
                           std::string& flexible,
                           std::string& side);
 
+/// Processes backend-specific ILU0 zero-pivot options
+static int processILU0Inputs(SystemSolver&      solver,
+                             const std::string& hw_backend,
+                             const CliOptions&  options);
+
 /// Main function selects example to be run
 int main(int argc, char* argv[])
 {
@@ -196,6 +206,12 @@ int sysGmres(int argc, char* argv[])
 
   solver.setGramSchmidtMethod(gs);
 
+  status = processILU0Inputs(solver, hw_backend, options);
+  if (status != 0)
+  {
+    return 1;
+  }
+
   // Read and open matrix and right-hand-side vector
   std::ifstream mat_file(matrix_pathname);
   if (!mat_file.is_open())
@@ -329,3 +345,69 @@ void processInputs(std::string& method, std::string& gs, std::string& sketch, st
     std::cout << "Setting preconditioning side to the default (right).\n\n";
   }
 }
+
+int processILU0Inputs(SystemSolver& solver, const std::string& hw_backend, const CliOptions& options)
+{
+  auto numeric_boost   = options.getParamFromKey("-n");
+  auto boost_tolerance = options.getParamFromKey("-t");
+  auto boost_value     = options.getParamFromKey("-v");
+  auto zero_diagonal   = options.getParamFromKey("-z");
+
+  if (hw_backend == "CPU")
+  {
+    if (numeric_boost || boost_tolerance || boost_value)
+    {
+      std::cout << "Options -n, -t, and -v are only supported by CUDA and HIP backends.\n"
+                << "For the CPU backend, use -z to set the zero diagonal value.\n";
+      return 1;
+    }
+
+    if (zero_diagonal)
+    {
+      if (zero_diagonal->second.empty())
+      {
+        std::cout << "Option -z requires a zero diagonal value.\n";
+        return 1;
+      }
+      return solver.getPreconditionerSolver().setCliParam("zero_diagonal", zero_diagonal->second);
+    }
+    return 0;
+  }
+
+  if (zero_diagonal)
+  {
+    std::cout << "Option -z is only supported by the CPU backend.\n"
+              << "For CUDA and HIP backends, use -n, -t, and -v.\n";
+    return 1;
+  }
+
+  int status = 0;
+  if (numeric_boost)
+  {
+    if ((numeric_boost->second != "yes") && (numeric_boost->second != "no"))
+    {
+      std::cout << "Numeric boost option must be 'yes' or 'no'.\n";
+      return 1;
+    }
+    status += solver.getPreconditionerSolver().setCliParam("numeric_boost", numeric_boost->second);
+  }
+  if (boost_tolerance)
+  {
+    if (boost_tolerance->second.empty())
+    {
+      std::cout << "Option -t requires a boost tolerance value.\n";
+      return 1;
+    }
+    status += solver.getPreconditionerSolver().setCliParam("boost_tolerance", boost_tolerance->second);
+  }
+  if (boost_value)
+  {
+    if (boost_value->second.empty())
+    {
+      std::cout << "Option -v requires a boost replacement value.\n";
+      return 1;
+    }
+    status += solver.getPreconditionerSolver().setCliParam("boost_value", boost_value->second);
+  }
+  return status;
+}

resolve/LinSolverDirectCpuILU0.cpp

@@ -21,6 +21,7 @@ namespace ReSolve
   LinSolverDirectCpuILU0::LinSolverDirectCpuILU0(LinAlgWorkspaceCpu* /* workspace */)
   // : workspace_(workspace)
   {
+    initParamList();
   }
 
   /**
@@ -381,23 +382,23 @@ namespace ReSolve
   }
 
   /**
-   * @brief Placeholder function for now.
+   * @brief Set Cli parameters for ILU0 solver.
    *
-   * The following switch (getParamId(Id)) cases always run the default and
-   * are currently redundant code (like an if (true)).
-   * In the future, they will be expanded to include more options.
+   * @param[in] id    - string ID for parameter to set
+   * @param[in] value - string value for parameter to set
    *
-   * @param id - string ID for parameter to set.
-   * @return int Error code.
+   * @return 0 if successful, 1 otherwise
    */
-  int LinSolverDirectCpuILU0::setCliParam(const std::string id, const std::string /* value */)
+  int LinSolverDirectCpuILU0::setCliParam(const std::string id, const std::string value)
   {
     switch (getParamId(id))
     {
+    case ZERO_DIAGONAL:
+      return setZeroDiagonal(atof(value.c_str()));
     default:
       std::cout << "Setting parameter failed!\n";
+      return 1;
     }
-    return 0;
   }
 
   /**
@@ -441,19 +442,18 @@ namespace ReSolve
   }
 
   /**
-   * @brief Placeholder function for now.
+   * @brief Get the real parameter for the ILU0 solver.
    *
-   * The following switch (getParamId(Id)) cases always run the default and
-   * are currently redundant code (like an if (true)).
-   * In the future, they will be expanded to include more options.
+   * @param[in] id - string ID for parameter to get
    *
-   * @param id - string ID for parameter to get.
-   * @return real_type Value of the real_type parameter to return.
+   * @return real_type - parameter value, or NaN if parameter is unknown
    */
   real_type LinSolverDirectCpuILU0::getCliParamReal(const std::string id) const
   {
     switch (getParamId(id))
     {
+    case ZERO_DIAGONAL:
+      return zero_diagonal_;
     default:
       out::error() << "Trying to get unknown real parameter " << id << "\n";
     }
@@ -480,15 +480,36 @@ namespace ReSolve
     return false;
   }
 
+  /**
+   * @brief Print the ILU0 Cli parameters.
+   *
+   * @param[in] id - string ID for parameter to print
+   *
+   * @return 0 if successful, 1 otherwise
+   */
   int LinSolverDirectCpuILU0::printCliParam(const std::string id) const
   {
     switch (getParamId(id))
     {
+    case ZERO_DIAGONAL:
+      std::cout << zero_diagonal_ << "\n";
+      break;
     default:
       out::error() << "Trying to print unknown parameter " << id << "\n";
       return 1;
     }
     return 0;
   }
 
+  /**
+   * @brief Initialize the parameter list for ILU0 solver.
+   *
+   * @post params_list_ is populated with the ILU0 solver parameters:
+   * - zero_diagonal
+   */
+  void LinSolverDirectCpuILU0::initParamList()
+  {
+    params_list_["zero_diagonal"] = ZERO_DIAGONAL;
+  }
+
 } // namespace ReSolve

resolve/LinSolverDirectCpuILU0.hpp

@@ -76,6 +76,13 @@ namespace ReSolve
     int         printCliParam(const std::string id) const override;
 
   private:
+    enum ParameterIDs
+    {
+      ZERO_DIAGONAL = 0
+    };
+
+    void initParamList();
+
     // MemoryHandler mem_; ///< Device memory manager object
     // LinAlgWorkspaceCpu* workspace_{nullptr};
 

resolve/LinSolverDirectCuDssRf.hpp

@@ -63,7 +63,7 @@ namespace ReSolve
     int  csc2csr(matrix::Csc* A_csc, matrix::Csr* A_csr);
 
   private:
-    enum ParamaterIDs
+    enum ParameterIDs
     {
       ZERO_PIVOT = 0,
       PIVOT_BOOST

resolve/LinSolverDirectCuSolverRf.hpp

@@ -59,7 +59,7 @@ namespace ReSolve
     int  csc2csr(matrix::Csc* A_csc, matrix::Csr* A_csr);
 
   private:
-    enum ParamaterIDs
+    enum ParameterIDs
     {
       ZERO_PIVOT = 0,
       PIVOT_BOOST