Using GMRES and KrylovPrecond classes

This tutorial demonstrates how to use the GMRES solver with a Krylov preconditioner to solve linear systems efficiently. We’ll walk through constructing a linear operator, creating instances of GMRES and KrylovPrecond, and solving a linear system. For more advanced usage and additional features, please refer to the API in lin-solve.hpp.

Step 1: Constructing a Linear Operator

First, we need to define our linear operator. This operator represents the matrix-vector multiplication in the linear system Ax = b. We’ll use a lambda function to define this operator.

auto LinOp = [&A](Vector<double>* Ax, const Vector<double>& x) {
    // Define the action of A on x and store the result in Ax
};

Note

GMRES passes an Ax that is already sized to x.Dim() and is backed by fixed-size scratch storage. Do not call Ax->ReInit(...) when the size already matches — it asserts. The safe idiom is:

if (Ax->Dim() != x.Dim()) Ax->ReInit(x.Dim());

Step 2: Creating Instances of GMRES and KrylovPrecond

Next, we’ll create instances of the GMRES solver and the Krylov preconditioner.

GMRES<double> solver;
KrylovPrecond<double> krylov_precond;

The GMRES constructor accepts two optional knobs that control the Arnoldi orthogonalization step:

// Classical Gram-Schmidt with one reorthogonalization pass — recommended
// for distributed (MPI) runs: one batched Allreduce per iteration, with
// MGS-equivalent stability.
GMRES<double> solver(comm, /*verbose=*/true,
                     GMRES<double>::GramSchmidt::CGS, /*num_reorth=*/1);

The default (MGS, 0) matches the legacy behavior and is a good choice for serial or shared-memory runs. See lin-solve.hpp for the full latency/stability trade-off.

Step 3: Solving the Linear System

With our linear operator and solver instances ready, we can now solve the linear system Ax = b. We’ll call the GMRES solver with the necessary arguments.

Vector<double> x, b;
Real tol = 1e-10;
solver(&x, LinOp, b, tol, -1, false, nullptr, &krylov_precond);

The above call to the GMRES solver will initialize the Krylov preconditioner and solve the linear system. Subsequent solves with the same linear operator will benefit from the preconditioner and may require fewer iterations.

Note

Ensure that your linear operator and right-hand side vector b are properly defined before calling the solver.

Example Code

The complete example code is shown below:

#include "sctl.hpp"
using namespace sctl;

Matrix<double> LowRankMatrix(const Long N, const Long rank) {
  Matrix<double> M(N, N);
  M.SetZero();

  for (Long r = 0; r < rank; r++) {
    Matrix<double> U(N,1), Vt(1, N);
    for (auto& a : U ) a = drand48();
    for (auto& a : Vt) a = drand48();
    M += U * Vt * exp(log(machine_eps<double>())*r/rank);
  }
  return M;
}

double max_norm(const Vector<double>& v) {
  double max_val = 0;
  for (const auto a : v) max_val = std::max<double>(max_val, fabs(a));
  return max_val;
}

int main(int argc, char** argv) {
  Comm::MPI_Init(&argc, &argv);

  { // Example usage for GMRES and KrylovPrecond
    Long N = 200, rank = 200;

    // Build A = I + <low-rank>
    Matrix<double> A = LowRankMatrix(N, rank);
    for (Long i = 0; i < N; i++) A[i][i] += 1;

    // Build linear operator
    auto LinOp = [&A](Vector<double>* Ax, const Vector<double>& x) {
      const Long N = x.Dim();
      if (Ax->Dim() != N) Ax->ReInit(N);
      Matrix<double> Ax_(N, 1, Ax->begin(), false);
      Ax_ = A * Matrix<double>(N, 1, (Iterator<double>)x.begin(), false);
    };

    // Set exact solution x0 and the RHS b := A * x
    Vector<double> x0(N), b(N), x;
    for (auto& a : x0) a = drand48();
    LinOp(&b, x0);

    // Solve using GMRES
    GMRES<double> solver;
    KrylovPrecond<double> krylov_precond;
    solver(&x, LinOp, b, 1e-10, -1, false, nullptr, &krylov_precond);
    std::cout<<"Solution error = "<<max_norm(x-x0)<<"\n\n";

    // Solve a new problem, reusing the Krylov preconditioner
    for (auto& a : x0) a = drand48();
    LinOp(&b, x0);
    solver(&x, LinOp, b, 1e-10, -1, false, nullptr, &krylov_precond);
    std::cout<<"Solution error = "<<max_norm(x-x0)<<"\n\n";
  }

  Comm::MPI_Finalize();
  return 0;
}

This generates the following output, where the first solve requires 23 GMRES iterations, while the second solve converges just 8 iterations due to preconditioning.

KSP Residual norm 2.099176367385e+03
KSP Residual norm 9.003181335456e+00
KSP Residual norm 4.229781314809e+00
KSP Residual norm 3.014948102699e+00
KSP Residual norm 2.848098976816e+00
KSP Residual norm 2.353313340801e+00
KSP Residual norm 1.317737385363e+00
KSP Residual norm 1.299136900499e+00
KSP Residual norm 5.830786318937e-01
KSP Residual norm 4.104966218497e-01
KSP Residual norm 2.296154200197e-01
KSP Residual norm 1.132160334702e-01
KSP Residual norm 7.357733656165e-02
KSP Residual norm 6.190744358418e-02
KSP Residual norm 5.755616763787e-02
KSP Residual norm 5.741807896259e-02
KSP Residual norm 1.680949998420e-02
KSP Residual norm 1.816481776180e-03
KSP Residual norm 6.359376010342e-04
KSP Residual norm 1.981085110869e-04
KSP Residual norm 5.229997256823e-05
KSP Residual norm 1.046989208963e-05
KSP Residual norm 1.443750846938e-06
KSP Residual norm 1.212803427478e-07
Solution error = 1.67549e-06

KSP Residual norm 1.834616605684e+01
KSP Residual norm 7.976544634403e+00
KSP Residual norm 1.159404821738e+00
KSP Residual norm 2.388777265664e-02
KSP Residual norm 4.215249431341e-03
KSP Residual norm 1.516499740199e-04
KSP Residual norm 1.152629619068e-05
KSP Residual norm 1.254922886814e-06
KSP Residual norm 1.542127320799e-08
Solution error = 5.41102e-08