Support left preconditioning in GMRES

CHANGELOG.md

@@ -80,4 +80,6 @@ It is seamless from the user perspective and fixed many bugs.
 
 16. Updated MatrixHandler::addConst to return integer error codes instead of void.
 
-17. Added a preconditioner interface class so users can define thier own preconditioners.
+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.

examples/ExampleHelper.hpp

@@ -141,12 +141,12 @@ namespace ReSolve
         }
 
         res_->copyFromExternal(r_, memspace_, memspace_);
-        real_type norm  = computeResidualNorm(*A_, *x_, *res_, memspace_);
-        real_type rnorm = norm2(*r_, memspace_);
+        real_type norm     = computeResidualNorm(*A_, *x_, *res_, memspace_);
+        real_type rhs_norm = norm2(*r_, memspace_);
 
         std::cout << "\t2-Norm of the residual: "
                   << std::scientific << std::setprecision(16)
-                  << norm / rnorm << "\n";
+                  << norm / rhs_norm << "\n";
       }
 
       /// Summary of direct solve

examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp

@@ -186,8 +186,8 @@ int main(int argc, char* argv[])
                   << status << std::endl;
         vec_rhs->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
         status = Rf->solve(vec_rhs, vec_x);
-        ReSolve::PreconditionerLU precond_lu(Rf);
-        FGMRES->setPreconditioner(&precond_lu);
+        ReSolve::PreconditionerLU preconditioner(Rf);
+        FGMRES->setPreconditioner(&preconditioner);
       }
       // if (i%2!=0)  vec_x->setToZero(ReSolve::memory::DEVICE);
       real_type norm_x = vector_handler->dot(vec_x, vec_x, ReSolve::memory::DEVICE);

examples/gpuRefactor.cpp

@@ -297,8 +297,8 @@ int gpuRefactor(int argc, char* argv[])
       {
         // Setup iterative refinement
         FGMRES.resetMatrix(A);
-        ReSolve::PreconditionerLU precond_lu(&Rf);
-        FGMRES.setPreconditioner(&precond_lu);
+        ReSolve::PreconditionerLU preconditioner(&Rf);
+        FGMRES.setPreconditioner(&preconditioner);
 
         // If refactorization produced finite solution do iterative refinement
         if (std::isfinite(helper.getNormRelativeResidual()))

examples/kluFactor.cpp

@@ -191,8 +191,8 @@ int main(int argc, char* argv[])
     {
       // Setup iterative refinement
       FGMRES.setup(A);
-      ReSolve::PreconditionerLU precond_lu(&KLU);
-      FGMRES.setPreconditioner(&precond_lu);
+      ReSolve::PreconditionerLU preconditioner(&KLU);
+      FGMRES.setPreconditioner(&preconditioner);
 
       // If refactorization produced finite solution do iterative refinement
       if (std::isfinite(helper.getNormRelativeResidual()))

examples/kluRefactor.cpp

@@ -194,8 +194,8 @@ int main(int argc, char* argv[])
     {
       // Setup iterative refinement
       FGMRES.setup(A);
-      ReSolve::PreconditionerLU precond_lu(KLU);
-      FGMRES.setPreconditioner(&precond_lu);
+      ReSolve::PreconditionerLU preconditioner(KLU);
+      FGMRES.setPreconditioner(&preconditioner);
 
       // If refactorization produced finite solution do iterative refinement
       if (std::isfinite(helper.getNormRelativeResidual()))

examples/randGmres.cpp

@@ -116,7 +116,7 @@ int runGmresExample(int argc, char* argv[])
 
   GramSchmidt GS(&vector_handler, GramSchmidt::CGS2);
 
-  precon_type                  Precond(&workspace);
+  precon_type                  precondition_solver(&workspace);
   LinSolverIterativeRandFGMRES FGMRES(&matrix_handler,
                                       &vector_handler,
                                       LinSolverIterativeRandFGMRES::cs,
@@ -170,15 +170,16 @@ int runGmresExample(int argc, char* argv[])
 
   matrix_handler.setValuesChanged(true, memspace);
 
-  Precond.setup(A);
   FGMRES.setRestart(150);
   FGMRES.setMaxit(2500);
   FGMRES.setTol(1e-12);
   FGMRES.setup(A);
-
   FGMRES.resetMatrix(A);
-  ReSolve::PreconditionerLU precond_lu(&Precond);
-  FGMRES.setPreconditioner(&precond_lu);
+
+  ReSolve::PreconditionerLU preconditioner(&precondition_solver);
+  preconditioner.setup(A);
+
+  FGMRES.setPreconditioner(&preconditioner);
   FGMRES.setFlexible(1);
 
   FGMRES.solve(vec_rhs, vec_x);

examples/sysGmres.cpp

@@ -35,6 +35,7 @@ void printHelpInfo()
   std::cout << "\t-g <gs method> \tGram-Schmidt method: cgs1, cgs2, or mgs (default 'cgs2').\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";
 }
 
 //
@@ -56,7 +57,8 @@ static int sysGmres(int argc, char* argv[]);
 static void processInputs(std::string& method,
                           std::string& gs,
                           std::string& sketch,
-                          std::string& flexible);
+                          std::string& flexible,
+                          std::string& side);
 
 /// Main function selects example to be run
 int main(int argc, char* argv[])
@@ -164,7 +166,10 @@ int sysGmres(int argc, char* argv[])
   opt                  = options.getParamFromKey("-x");
   std::string flexible = opt ? (*opt).second : "yes";
 
-  processInputs(method, gs, sketch, flexible);
+  opt              = options.getParamFromKey("-p");
+  std::string side = opt ? (*opt).second : "right";
+
+  processInputs(method, gs, sketch, flexible, side);
 
   std::cout << "Matrix file: " << matrix_pathname << "\n"
             << "RHS file: " << rhs_pathname << "\n";
@@ -247,7 +252,7 @@ int sysGmres(int argc, char* argv[])
   // Set up the preconditioner
   if (return_code == 0)
   {
-    status = solver.preconditionerSetup();
+    status = solver.preconditionerSetup(side);
     std::cout << "solver.preconditionerSetup returned status: " << status << "\n";
     if (status != 0)
     {
@@ -281,7 +286,7 @@ int sysGmres(int argc, char* argv[])
   return return_code;
 }
 
-void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& flexible)
+void processInputs(std::string& method, std::string& gs, std::string& sketch, std::string& flexible, std::string& side)
 {
   if (method == "randgmres")
   {
@@ -314,4 +319,10 @@ void processInputs(std::string& method, std::string& gs, std::string& sketch, st
     std::cout << "Setting flexible to the default (yes).\n\n";
     flexible = "yes";
   }
+
+  if ((side != "left") && (side != "right"))
+  {
+    std::cout << "Preconditioning side " << side << " not recognized.\n";
+    std::cout << "Setting preconditioning side to the default (right).\n\n";
+  }
 }

resolve/CMakeLists.txt

@@ -21,6 +21,7 @@ set(ReSolve_SRC
     SystemSolver.cpp
     Preconditioner.cpp
     PreconditionerLU.cpp
+    PreconditionerUserMatrix.cpp
 )
 
 set(ReSolve_KLU_SRC LinSolverDirectKLU.cpp)
@@ -52,6 +53,7 @@ set(ReSolve_HEADER_INSTALL
     MemoryUtils.hpp
     Preconditioner.hpp
     PreconditionerLU.hpp
+    PreconditionerUserMatrix.hpp
 )
 
 set(ReSolve_KLU_HEADER_INSTALL LinSolverDirectKLU.hpp)