Support left preconditioning in GMRES by kakeueda · Pull Request #417 · ORNL/ReSolve

kakeueda · 2026-03-04T21:32:36Z

Description

This change allows users to solve a left-preconditioned system $M^{-1}Ax = M^{-1}b$ using GMRES. It also adds a PreconditionerUserMatrix class.

@JeffZ594 My understanding is that we eventually want to add regularization to the GMRES solvers. We can follow up with another PR that adds regularization support with testTripsProblem.

Proposed changes

Updated GMRES to support left preconditioning.
Added a PreconditionerUserMatrix class so users can define $M^{-1}$ explicitly.
Added unit tests for lu and usermatrix preconditioners.
Added a new command line option -p left (or right) for the GMRES solver (default: right).

Checklist

Further comments

Non-flexible GMRES only. If FGMRES is selected with -p left, the solver reports an error.

Copilot

Pull request overview

This PR extends ReSolve’s GMRES implementations to support left preconditioning (M^{-1}Ax = M^{-1}b), adds a user-supplied explicit preconditioner matrix via PreconditionerUserMatrix, and wires the feature through examples/tests and CLI (-p left|right).

Changes:

Add Preconditioner::getSide()/setSide() and plumb side selection through SystemSolver::preconditionerSetup(...).
Implement left-preconditioned variants of the GMRES/Randomized GMRES iteration logic (with fallback to right preconditioning for flexible GMRES).
Add unit/functionality tests and example updates for left preconditioning and the new user-matrix preconditioner.

Reviewed changes

Copilot reviewed 24 out of 24 changed files in this pull request and generated 11 comments.

Show a summary per file

File	Description
`resolve/Preconditioner.hpp`	Adds preconditioning side API (`left`/`right`).
`resolve/Preconditioner.cpp`	Implements `getSide()` / `setSide()`.
`resolve/SystemSolver.hpp`	Extends preconditioner setup signature and exposes `getPreconditioner()`.
`resolve/SystemSolver.cpp`	Applies preconditioner setup with chosen side; exposes preconditioner getter.
`resolve/PreconditionerUserMatrix.hpp`	Declares user-provided explicit preconditioner matrix wrapper.
`resolve/PreconditionerUserMatrix.cpp`	Implements user-matrix preconditioner via `MatrixHandler::matvec`.
`resolve/LinSolverIterativeFGMRES.hpp/.cpp`	Implements left-preconditioning behavior (non-flexible GMRES) and warns/falls back for flexible GMRES.
`resolve/LinSolverIterativeRandFGMRES.hpp/.cpp`	Implements left-preconditioning behavior in randomized GMRES variant.
`resolve/PreconditionerLU.hpp`	Updates file header comment.
`resolve/PreconditionerLU.cpp`	Updates documentation strings/comments for LU preconditioner.
`resolve/CMakeLists.txt`	Adds new preconditioner source/header to build/install.
`tests/unit/preconditioner/*`	Adds LU and user-matrix preconditioner unit tests and CMake integration.
`tests/unit/CMakeLists.txt`	Adds preconditioner unit test subdirectory.
`tests/functionality/testSysGmres.cpp`	Adds `-p` parsing and conditional residual checking for left preconditioning.
`tests/functionality/CMakeLists.txt`	Adds functionality tests covering left-preconditioned runs.
`examples/sysGmres.cpp`	Adds `-p` parsing and passes side into `preconditionerSetup`.
`examples/randGmres.cpp`	Refactors preconditioner setup to use `PreconditionerLU` consistently.
`CHANGELOG.md`	Adds entry documenting left preconditioning + user-matrix preconditioner.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2026-03-04T22:08:10Z

+    matrix_handler_->matvec(B_, rhs, x, &ONE, &ZERO, memspace_);
+    return 0;
+  }
+
+  void PreconditionerUserMatrix::setMemorySpace()
+  {
+    bool is_matrix_handler_cuda = matrix_handler_->getIsCudaEnabled();
+    bool is_matrix_handler_hip  = matrix_handler_->getIsHipEnabled();
+


PreconditionerUserMatrix::apply() unconditionally calls MatrixHandler::matvec() and then returns 0. If B_ or matrix_handler_ is null this will crash, and any nonzero status from matvec is currently ignored. Please add null checks (e.g., require setPrecMatrix() to have been called) and propagate/return the matvec status so callers can detect failures.

Suggested change

matrix_handler_->matvec(B_, rhs, x, &ONE, &ZERO, memspace_);

return 0;

}

void PreconditionerUserMatrix::setMemorySpace()

{

bool is_matrix_handler_cuda = matrix_handler_->getIsCudaEnabled();

bool is_matrix_handler_hip = matrix_handler_->getIsHipEnabled();

// Ensure that the matrix handler and preconditioner matrix have been set

if (matrix_handler_ == nullptr || B_ == nullptr)

{

// Nonzero return value indicates failure (e.g., missing setup)

return 1;

}

int status = matrix_handler_->matvec(B_, rhs, x, &ONE, &ZERO, memspace_);

return status;

}

void PreconditionerUserMatrix::setMemorySpace()

{

bool is_matrix_handler_cuda = false;

bool is_matrix_handler_hip = false;

if (matrix_handler_ != nullptr)

{

is_matrix_handler_cuda = matrix_handler_->getIsCudaEnabled();

is_matrix_handler_hip = matrix_handler_->getIsHipEnabled();

}

JeffZ594 · 2026-03-09T17:29:47Z

@kakeueda

Thank you for making the unit tests for the preconditioners! I have code that implements regularization in Eigen that needs to be translated that would be better for testTripsProblem. I think these tests work better for testing just the preconditioners.

kakeueda · 2026-03-10T20:24:13Z

@kakeueda

Thank you for making the unit tests for the preconditioners! I have code that implements regularization in Eigen that needs to be translated that would be better for testTripsProblem. I think these tests work better for testing just the preconditioners.

Sure, let me first finish this PR and I'll take a look at your implementation .
@pelesh Are you okay with adding Eigen as a new dependency?

JeffZ594 · 2026-03-10T20:27:43Z

@kakeueda I believe Dr Peles does not want to use Eigen as a dependency. Right now I only use Eigen to calculate the SVD of the Hessenberg Matrix for GMRES but I believe there are alternatives like LAPACK that can also do this.

pelesh · 2026-03-18T22:00:57Z

@kakeueda I believe Dr Peles does not want to use Eigen as a dependency. Right now I only use Eigen to calculate the SVD of the Hessenberg Matrix for GMRES but I believe there are alternatives like LAPACK that can also do this.

Hessenberg matrix is dense, so SVD from LAPACK should work in this case. I wouldn't add Eigen as a dependency unless there is compelling enough reason to do so.

Co-authored-by: Kakeru Ueda <kakecoco2580@icloud.com>

Fixed some small bugs

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

kakeueda · 2026-03-19T13:24:02Z

For now the default behavior for left preconditioning is to use the preconditioned residual $||M^{-1}(b - Ax)||$ for convergence check and the true residual $||(b - Ax)||$ for monitoring (final_residual_norm_ stores this).

kakeueda · 2026-03-19T13:30:37Z

This is outside the scope of this PR but 'solve()' in GMRES is starting to get a bit cluttered. It might be worth considering adding a base GMRES class and having the current RandFGMRES and regularized GMRES derive from it. We could also clean up the convergence check for GMRES which I believe this was mentioned before.

shakedregev · 2026-03-19T14:35:45Z

+        true_rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_);
+        true_rnorm = std::sqrt(true_rnorm);


maybe we should change this to true_res_norm? r can be anything, it's not clear from context (and there's no separator from norm)

Agree, also changed bnorm to rhs_norm

shakedregev · 2026-03-19T14:43:34Z

+      /**
+       * @brief Test apply() performs correct matvec with B.
+       *
+       * B = diag(2, 3, 4), x = [1, 1, 1]^T => y = [2, 3, 4]^T.


Why tranpose?

It was meant to show a column vector, but maybe unnecessary

Worse, it's misleading. Vectors are by default column, unless you put a transpose. That's what didn't make sense to me. I knew it had to be a column, but this is the convention for row.

Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com>

shakedregev

Some minor comments. I need to test and I will change this to approve.

shakedregev

Great work.

kakeueda · 2026-03-19T17:20:11Z

returned 503: Egress is over the account limit.

I got this error in the test. Any idea what might be causing this? @pelesh

nkoukpaizan · 2026-03-19T17:42:43Z

returned 503: Egress is over the account limit.

I got this error in the test. Any idea what might be causing this? @pelesh

Looks like a temporary failure in GH-Actions. It's working now!

pelesh

Very nice work and good testing coverage of the new preconditioner features. A few things that need to be done before we merge this:

The GMRES code needs to be better commented to explain the workflow.
Use enums instead of strings to select preconditioner type.
There should be a warning message if user tries to use e.g. FGMRES with unsupported preconditioner and the preconditioner choice should fall back to default.

Also see some more specific comments below.

pelesh · 2026-03-26T20:10:14Z

+    // Left preconditioning uses preconditioned norms for convergence
+    if (!flexible_ && preconditioner_->getSide() == "left")
+    {
+      // Left-preconditioned residual norm ||M^{-1}*(b - A*x0)||
+      preconditioner_->apply(&vec_v, &vec_z);
+      vec_v.copyFromExternal(&vec_z, memspace_, memspace_);
+      res_norm = vector_handler_->dot(vec_V_, vec_V_, memspace_);
+      res_norm = std::sqrt(res_norm);
+
+      // Left-preconditioned right-hand side norm ||M^{-1}*b||
+      vec_v.setData(rhs->getData(memspace_), memspace_);
+      preconditioner_->apply(&vec_v, &vec_z);
+      rhs_norm = vector_handler_->dot(&vec_z, &vec_z, memspace_);
+      rhs_norm = std::sqrt(rhs_norm);
+    }
+    else
+    {
+      res_norm = true_res_norm;
+      rhs_norm = true_rhs_norm;
+    }


Why is left preconditioner not applied to flexible GMRES?

FGMRES is generally defined as a right-preconditioned method. Please see Section 9.4 (https://www.mat.uc.pt/~alma/aulas/NewTrends/Books/Saad.pd). One difficulty with left preconditioning in the flexible setting is that the residual $M^{-1}(Ax-b)$ can change at each iteration as the preconditioner changes. petsc also supports only right preconditioning for FGMRES.

pelesh · 2026-03-26T20:15:18Z

+      vec_v.setData(rhs->getData(memspace_), memspace_);
+      preconditioner_->apply(&vec_v, &vec_z);


Why use setData function here? This method will be removed and should not be used in new contributions.

Is there currently an equivalent to setData that we should use instead?
I thought we would still need setData because vec_v should be a view into the multivector vec_V_, not a copy.

pelesh · 2026-03-26T20:20:11Z

+        if (preconditioner_->getSide() == "right")
+        {


It is safer to use here enums instead of strings.

pelesh · 2026-03-26T20:22:11Z

+    real_type   true_res_norm = 0.0;
+    real_type   true_rhs_norm = 0.0;


What is "true" residual norm? This needs to be documented in the code.

Added documentation explaining what the true residual is. This separates the residual used for convergence checking from the residual reported to the user. For right preconditioning the true residual $Ax - b$ is used for both, but for left preconditioning the preconditioned residual $M^{-1}(Ax - b)$ is used for convergence checking, while $Ax - b$ is still reported so users can monitor the residual of the original system $Ax = b$. We could also add an option to report the preconditioned residual but not the scope of this PR because TestHelper would also need to be updated to compute the preconditioned residual.

pelesh · 2026-03-26T20:25:14Z

+          if (preconditioner_->getSide() == "right")
+          {


I suggest using sitch statement and enums instead of strings and if-else statements here.

pelesh · 2026-03-26T20:26:25Z

+          // The correction is M^{-1}*vec_z for right preconditioning
+          preconditioner_->apply(&vec_z, &vec_v);
+          // Add the correction to x
+          vector_handler_->axpy(ONE, &vec_v, x, memspace_);
+        }
+        else
+        {
+          // Add the correction to x directly for left preconditioning.
+          vector_handler_->axpy(ONE, &vec_z, x, memspace_);


The comments in the code need to be more specific explaining steps clearly.

pelesh · 2026-03-26T20:27:16Z

+        // Report the true relative residual norm
+        final_residual_norm_ = true_res_norm / true_rhs_norm;


Really need to explain what is "true" residual norm.

pelesh · 2026-03-26T20:31:14Z

+add_test(NAME sys_rand_count_gmres_left_cgs2_test
+         COMMAND $<TARGET_FILE:sys_rand_gmres_test.exe> "-x" "no" "-i"
+                 "randgmres" "-g" "cgs2" "-s" "count" "-p" "left"
+)
+add_test(NAME sys_rand_count_gmres_left_mgs_test
+         COMMAND $<TARGET_FILE:sys_rand_gmres_test.exe> "-x" "no" "-i"
+                 "randgmres" "-g" "mgs" "-s" "count" "-p" "left"
+)
+add_test(NAME sys_gmres_left_cgs2_test
+         COMMAND $<TARGET_FILE:sys_rand_gmres_test.exe> "-x" "no" "-i" "fgmres"
+                 "-g" "cgs2" "-p" "left"
+)
+add_test(NAME sys_gmres_left_mgs_test
+         COMMAND $<TARGET_FILE:sys_rand_gmres_test.exe> "-x" "no" "-i" "fgmres"
+                 "-g" "mgs" "-p" "left"
+)


We should consider adding more combinations of different preconditioners and other solver options.

kakeueda · 2026-04-02T13:16:44Z

@pelesh Thank you for your suggestions.

Replaced std::string with an enum for the preconditioning side
Improved the comments explaining the workflow and residuals
More tests are added for left preconditioning
Changed the behavior to fall back to the default option (right) instead of throwing an error when FGMRES and left preconditioning are selected.

Copilot AI review requested due to automatic review settings March 4, 2026 21:32

kakeueda self-assigned this Mar 4, 2026

Copilot started reviewing on behalf of kakeueda March 4, 2026 21:33 View session

kakeueda requested review from Copilot and removed request for Copilot March 4, 2026 21:35

Copilot started reviewing on behalf of kakeueda March 4, 2026 22:00 View session

Copilot AI reviewed Mar 4, 2026

View reviewed changes

kakeueda marked this pull request as draft March 4, 2026 23:52

JeffZ594 and others added 18 commits March 19, 2026 15:06

ABBA GMRES Commit

f469ac5

Added check for preconditioner being set

d7716e2

Use Preconditioner Class for ABBAGMRES

7ab3e29

Apply pre-commmit fixes

df69725

Apply pre-commmit fixes

d8c6f39

Update resolve/PreconditionerMatvec.cpp

0104791

Co-authored-by: Kakeru Ueda <kakecoco2580@icloud.com>

Update resolve/LinSolverIterativeFGMRES.cpp

1032c26

Co-authored-by: Kakeru Ueda <kakecoco2580@icloud.com>

Update resolve/Preconditioner.hpp

d2a590e

Co-authored-by: Kakeru Ueda <kakecoco2580@icloud.com>

Apply pre-commmit fixes

4bb1fac

Renamed PreconditionerMatvec to PreconditionerABBA

0588909

Fixed some small bugs

Apply pre-commmit fixes

8c899bd

Fix left preconditioning and add corresponding tests

e620de1

Fix typo and some comments

8e55d7d

Apply precommit and update CHANGELOG.md

2714e54

Update CHANGELOG.md

b0f3ebb

Use preconditioned bnorm in RandFGMRES

8be5465

Apply pre-commmit fixes

68807d0

Remove unnecessary variables

f34279f

kakeueda and others added 3 commits March 19, 2026 15:06

Apply pre-commmit fixes

25d85e8

[skip ci] Update CHANGELOG.md

ca953b9

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

Fix norm evaluation in GMRES

4d36a83

kakeueda force-pushed the kakeru/preconditioner-ABBA branch from e866f59 to 4d36a83 Compare March 19, 2026 06:09

Fix naming and error check

03175f7

ORNL deleted a comment from Copilot AI Mar 19, 2026

kakeueda added 2 commits March 19, 2026 21:45

[skip ci] Move #pragma above

1c28d2e

[skip ci] Fix typo

4238c7f

kakeueda marked this pull request as ready for review March 19, 2026 13:30

kakeueda requested review from pelesh and shakedregev March 19, 2026 13:31

shakedregev reviewed Mar 19, 2026

View reviewed changes

Comment thread resolve/LinSolverIterativeFGMRES.cpp Outdated

shakedregev reviewed Mar 19, 2026

View reviewed changes

Comment thread resolve/LinSolverIterativeRandFGMRES.cpp Outdated

shakedregev reviewed Mar 19, 2026

View reviewed changes

Apply suggestions from code review

258c058

Co-authored-by: Shaked Regev <35384901+shakedregev@users.noreply.github.com>

shakedregev reviewed Mar 19, 2026

View reviewed changes

shakedregev approved these changes Mar 19, 2026

View reviewed changes

Fix naming and comment

581a946

pelesh requested changes Mar 26, 2026

View reviewed changes

kakeueda added 2 commits April 2, 2026 21:48

Refactor preconditioner side to use enum

36cebb1

Apply formatting

d83452c

-    matrix_handler_->matvec(B_, rhs, x, &ONE, &ZERO, memspace_);
-    return 0;
-  }
-  void PreconditionerUserMatrix::setMemorySpace()
-  {
-    bool is_matrix_handler_cuda = matrix_handler_->getIsCudaEnabled();
-    bool is_matrix_handler_hip  = matrix_handler_->getIsHipEnabled();
+    // Ensure that the matrix handler and preconditioner matrix have been set
+    if (matrix_handler_ == nullptr || B_ == nullptr)
+    {
+      // Nonzero return value indicates failure (e.g., missing setup)
+      return 1;
+    }
+    int status = matrix_handler_->matvec(B_, rhs, x, &ONE, &ZERO, memspace_);
+    return status;
+  }
+  void PreconditionerUserMatrix::setMemorySpace()
+  {
+    bool is_matrix_handler_cuda = false;
+    bool is_matrix_handler_hip  = false;
+    if (matrix_handler_ != nullptr)
+    {
+      is_matrix_handler_cuda = matrix_handler_->getIsCudaEnabled();
+      is_matrix_handler_hip  = matrix_handler_->getIsHipEnabled();
+    }

		true_rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_);
		true_rnorm = std::sqrt(true_rnorm);

		vec_v.setData(rhs->getData(memspace_), memspace_);
		preconditioner_->apply(&vec_v, &vec_z);

		real_type true_res_norm = 0.0;
		real_type true_rhs_norm = 0.0;

		// Report the true relative residual norm
		final_residual_norm_ = true_res_norm / true_rhs_norm;

Conversation

kakeueda commented Mar 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Proposed changes

Checklist

Further comments

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Copilot AI Mar 4, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

JeffZ594 commented Mar 9, 2026

Uh oh!

kakeueda commented Mar 10, 2026

Uh oh!

JeffZ594 commented Mar 10, 2026

Uh oh!

pelesh commented Mar 18, 2026

Uh oh!

kakeueda commented Mar 19, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

kakeueda commented Mar 19, 2026

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

shakedregev left a comment

Choose a reason for hiding this comment

Uh oh!

shakedregev left a comment

Choose a reason for hiding this comment

Uh oh!

kakeueda commented Mar 19, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

nkoukpaizan commented Mar 19, 2026

Uh oh!

pelesh left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kakeueda Apr 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kakeueda commented Mar 4, 2026 •

edited

Loading

kakeueda commented Mar 19, 2026 •

edited

Loading

kakeueda commented Mar 19, 2026 •

edited

Loading

kakeueda Apr 2, 2026 •

edited

Loading

kakeueda commented Apr 2, 2026 •

edited

Loading