GCR (Generalised Conjugate Residual)#

gko::solver::Gcr<ValueType> is a Krylov method for general non-symmetric systems that, like GMRES, minimises the 2-norm of the residual over the Krylov subspace. Where GMRES does so through an Arnoldi basis and a Hessenberg least-squares solve, GCR maintains an explicit set of search directions \(p_0, p_1, \ldots\) chosen so that \(\{A p_k\}\) is mutually orthogonal,

\[ \langle A p_i,\, A p_j \rangle = 0 \quad \text{for } i \ne j. \]

That orthogonality lets each iteration update the residual with a single projection along \(A p_k\). The trade-off is that the search directions all have to be remembered to keep them orthogonal, so memory grows linearly with the iteration count, the same as GMRES.

When to use GCR#

  • Non-symmetric systems where you want GMRES-style residual minimisation but prefer working with the search directions explicitly, e.g. when pairing with a Krylov-aware preconditioner, or when an analysis pass needs the \(A p_k\) vectors directly.

  • A restartable alternative to GMRES with comparable convergence on smooth problems. Set krylov_dim (default 100) to bound the working memory.

In most situations GMRES is faster; its Arnoldi basis and the small Hessenberg solve are typically cheaper per iteration than GCR’s explicit search-direction orthogonalisation. Reach for GCR when you have a specific reason to want \(A p_k\)-orthogonality rather than \(V_k\)-orthogonality.

For non-symmetric systems where bounded memory matters more than residual minimisation, BiCGSTAB is a short-recurrence alternative.

Construction#

auto gcr = gko::solver::Gcr<double>::build()
               .with_criteria(
                   gko::stop::Iteration::build()
                       .with_max_iters(1000u)
                       .on(exec),
                   gko::stop::ResidualNorm<double>::build()
                       .with_reduction_factor(1e-10)
                       .on(exec))
               .with_krylov_dim(100u)
               .with_preconditioner(
                   gko::preconditioner::Jacobi<double, gko::int32>::build()
                       .on(exec))
               .on(exec);

auto solver = gcr->generate(system_matrix);
solver->apply(b, x);

Factory parameters#

GCR inherits the iterative-solver surface from EnablePreconditionedIterativeSolver and adds one tuning knob, krylov_dim:

Parameter

Type

Purpose

criteria

vector of stop::CriterionFactory

Stopping criteria. Required; see Stopping criteria.

preconditioner

LinOpFactory

Preconditioner factory; rebuilt each generate() call. Defaults to identity.

generated_preconditioner

LinOp

An already-generated preconditioner. Bypasses the factory and is used directly.

krylov_dim

size_type

Maximum number of search directions kept before restart. Default 100. Larger values give better convergence per restart at the cost of linear-in-krylov_dim memory and orthogonalisation work.

Memory footprint#

GCR stores krylov_dim search directions \(p_k\) and krylov_dim preconditioned-and-applied vectors \(A p_k\): two vector arrays of length \(n\) each, both growing with the iteration count up to the restart limit. After a restart the buffers are reused.

Implementation notes#

The orthogonalisation uses modified Gram-Schmidt. Each iteration costs one SpMV, one preconditioner apply, \(O(k)\) inner products for the re-orthogonalisation against the previous search directions, and an AXPY update of the iterate and residual.

See also