Support unpreconditioned GMRES

CHANGELOG.md

@@ -4,6 +4,7 @@
 
 - Refactored system solver: Added helper functions to check and fix conflicts in configuration options and to reduce code duplication.
 - Code returns errors instead of warnings when encountering ambiguous erronneus inputs from the user.
+- Added an identity preconditioner for unpreconditioned GMRES.
 
 ## HyKKT Release changes
 

examples/randGmres.cpp

@@ -116,7 +116,7 @@ int runGmresExample(int argc, char* argv[])
 
   GramSchmidt GS(&vector_handler, GramSchmidt::CGS2);
 
-  precon_type                  precondition_solver(&workspace);
+  precon_type                  preconditioner_solver(&workspace);
   LinSolverIterativeRandFGMRES FGMRES(&matrix_handler,
                                       &vector_handler,
                                       LinSolverIterativeRandFGMRES::cs,
@@ -179,7 +179,7 @@ int runGmresExample(int argc, char* argv[])
   FGMRES.setup(A);
   FGMRES.resetMatrix(A);
 
-  ReSolve::PreconditionerLU preconditioner(&precondition_solver);
+  ReSolve::PreconditionerLU preconditioner(&preconditioner_solver);
   preconditioner.setup(A);
 
   FGMRES.setPreconditioner(&preconditioner);

examples/sysGmres.cpp

@@ -31,6 +31,7 @@ void printHelpInfo()
   std::cout << "sysGmres.exe -m <matrix file> -r <rhs file>\n\n";
   std::cout << "Optional features:\n";
   std::cout << "\t-b <cpu|cuda|hip> \tSelects hardware backend.\n";
+  std::cout << "\t-c <ilu0|none> \tSelects preconditioner (default 'ilu0').\n";
   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";
@@ -63,7 +64,8 @@ static void processInputs(std::string& method,
                           std::string& gs,
                           std::string& sketch,
                           std::string& flexible,
-                          std::string& side);
+                          std::string& side,
+                          std::string& preconditioner);
 
 /// Processes backend-specific ILU0 zero-pivot options
 static int processILU0Inputs(SystemSolver&      solver,
@@ -179,10 +181,14 @@ int sysGmres(int argc, char* argv[])
   opt              = options.getParamFromKey("-p");
   std::string side = opt ? (*opt).second : "right";
 
-  processInputs(method, gs, sketch, flexible, side);
+  opt                        = options.getParamFromKey("-c");
+  std::string preconditioner = opt ? (*opt).second : "ilu0";
+
+  processInputs(method, gs, sketch, flexible, side, preconditioner);
 
   std::cout << "Matrix file: " << matrix_pathname << "\n"
-            << "RHS file: " << rhs_pathname << "\n";
+            << "RHS file: " << rhs_pathname << "\n"
+            << "Preconditioner: " << (preconditioner == "none" ? "none (identity)" : "ILU0") << "\n";
 
   // Create workspace
   workspace_type workspace;
@@ -201,15 +207,18 @@ int sysGmres(int argc, char* argv[])
                       "none",
                       "none",
                       method,
-                      "ilu0",
+                      preconditioner,
                       "none");
 
   solver.setGramSchmidtMethod(gs);
 
-  status = processILU0Inputs(solver, hw_backend, options);
-  if (status != 0)
+  if (preconditioner == "ilu0")
   {
-    return 1;
+    status = processILU0Inputs(solver, hw_backend, options);
+    if (status != 0)
+    {
+      return 1;
+    }
   }
 
   // Read and open matrix and right-hand-side vector
@@ -268,8 +277,8 @@ int sysGmres(int argc, char* argv[])
   solver.getIterativeSolver().setCliParam("flexible", flexible);
   solver.getIterativeSolver().setCliParam("restart", "200");
 
-  // Set up the preconditioner
-  if (return_code == 0)
+  // Set up the requested non-identity preconditioner
+  if (return_code == 0 && preconditioner == "ilu0")
   {
     status = solver.preconditionerSetup(side);
     std::cout << "solver.preconditionerSetup returned status: " << status << "\n";
@@ -305,7 +314,12 @@ int sysGmres(int argc, char* argv[])
   return return_code;
 }
 
-void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& flexible, std::string& side)
+void processInputs(std::string& method,
+                   std::string& gs,
+                   std::string& sketch,
+                   std::string& flexible,
+                   std::string& side,
+                   std::string& preconditioner)
 {
   if (method == "randgmres")
   {
@@ -344,6 +358,13 @@ void processInputs(std::string& method, std::string& gs, std::string& sketch, st
     std::cout << "Preconditioning side " << side << " not recognized.\n";
     std::cout << "Setting preconditioning side to the default (right).\n\n";
   }
+
+  if (preconditioner != "ilu0" && preconditioner != "none")
+  {
+    std::cout << "Preconditioner " << preconditioner << " not recognized.\n";
+    std::cout << "Setting preconditioner to the default (ILU0).\n\n";
+    preconditioner = "ilu0";
+  }
 }
 
 int processILU0Inputs(SystemSolver& solver, const std::string& hw_backend, const CliOptions& options)

resolve/CMakeLists.txt

@@ -20,6 +20,7 @@ set(ReSolve_SRC
     LinSolverIterativeRandFGMRES.cpp
     SystemSolver.cpp
     Preconditioner.cpp
+    PreconditionerIdentity.cpp
     PreconditionerLU.cpp
     PreconditionerUserMatrix.cpp
 )
@@ -55,6 +56,7 @@ set(ReSolve_HEADER_INSTALL
     GramSchmidt.hpp
     MemoryUtils.hpp
     Preconditioner.hpp
+    PreconditionerIdentity.hpp
     PreconditionerLU.hpp
     PreconditionerUserMatrix.hpp
 )

resolve/PreconditionerIdentity.cpp

@@ -0,0 +1,75 @@
+/**
+ * @file   PreconditionerIdentity.cpp
+ * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp)
+ * @brief  Implementation of identity preconditioner class.
+ *
+ */
+
+#include "PreconditionerIdentity.hpp"
+
+#include <cassert>
+
+#include <resolve/MemoryUtils.hpp>
+#include <resolve/vector/Vector.hpp>
+
+namespace ReSolve
+{
+  /**
+   * @brief Constructs an identity preconditioner.
+   *
+   * @param[in] memspace Memory space used by input and output vectors.
+   */
+  PreconditionerIdentity::PreconditionerIdentity(memory::MemorySpace memspace)
+    : memspace_(memspace)
+  {
+  }
+
+  /**
+   * @brief Sets up the identity preconditioner.
+   *
+   * The identity preconditioner has no matrix-dependent state, so setup is a
+   * no-op and accepts a null matrix pointer.
+   *
+   * @return Always returns zero.
+   */
+  int PreconditionerIdentity::setup(matrix_type* /* A */)
+  {
+    return 0;
+  }
+
+  /**
+   * @brief Resets the identity preconditioner.
+   *
+   * The identity preconditioner has no matrix-dependent state, so reset is a
+   * no-op and accepts a null matrix pointer.
+   *
+   * @return Always returns zero.
+   */
+  int PreconditionerIdentity::reset(matrix_type* /* A */)
+  {
+    return 0;
+  }
+
+  /**
+   * @brief Applies the identity preconditioner.
+   *
+   * Copies rhs to x in the memory space specified at construction.
+   *
+   * @param[in]  rhs Input vector
+   * @param[out] x   Output vector
+   *
+   * @return Zero on success, one if the vectors are invalid or cannot be copied.
+   */
+  int PreconditionerIdentity::apply(vector_type* rhs, vector_type* x)
+  {
+    if (rhs == nullptr || x == nullptr)
+    {
+      return 1;
+    }
+
+    assert(rhs->getSize() == x->getSize() && "Vector sizes do not match!");
+    assert(rhs->getNumVectors() == x->getNumVectors() && "Vector counts do not match!");
+
+    return x->copyFromExternal(rhs, memspace_, memspace_) == 0 ? 0 : 1;
+  }
+} // namespace ReSolve

resolve/PreconditionerIdentity.hpp

@@ -0,0 +1,36 @@
+/**
+ * @file   PreconditionerIdentity.hpp
+ * @author Kakeru Ueda (k.ueda.2290@m.isct.ac.jp)
+ * @brief  Declaration of identity preconditioner class.
+ *
+ */
+
+#pragma once
+
+#include <resolve/MemoryUtils.hpp>
+#include <resolve/Preconditioner.hpp>
+
+namespace ReSolve
+{
+  /**
+   * @brief Identity preconditioner.
+   *
+   * This preconditioner uses the identity matrix,
+   * so applying it copies rhs to x.
+   */
+  class PreconditionerIdentity : public Preconditioner
+  {
+  public:
+    using vector_type = vector::Vector;
+    using matrix_type = matrix::Sparse;
+
+    explicit PreconditionerIdentity(memory::MemorySpace memspace);
+
+    int setup(matrix_type* A) override;
+    int reset(matrix_type* A) override;
+    int apply(vector_type* rhs, vector_type* x) override;
+
+  private:
+    memory::MemorySpace memspace_;
+  };
+} // namespace ReSolve