(x, y, stats) = gpmr(A, B, b::AbstractVector{FC}, c::AbstractVector{FC}; memory::Int=20,
                     C=I, D=I, E=I, F=I, atol::T=√eps(T), rtol::T=√eps(T),
                     gsp::Bool=false, reorthogonalization::Bool=false,
                     itmax::Int=0, λ::FC=one(FC), μ::FC=one(FC),
                     verbose::Int=0, history::Bool=false,
                     ldiv::Bool=false, callback=solver->false)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

GPMR solves the unsymmetric partitioned linear system

[ λI   A ] [ x ] = [ b ]
[  B  μI ] [ y ]   [ c ],

where λ and μ are real or complex numbers. A can have any shape and B has the shape of Aᵀ. A, B, b and c must be all nonzero.

This implementation allows left and right block diagonal preconditioners

[ C    ] [ λM   A ] [ E    ] [ E⁻¹x ] = [ Cb ]
[    D ] [  B  μN ] [    F ] [ F⁻¹y ]   [ Dc ],

and can solve

[ λM   A ] [ x ] = [ b ]
[  B  μN ] [ y ]   [ c ]

when CE = M⁻¹ and DF = N⁻¹.

By default, GPMR solves unsymmetric linear systems with λ = 1 and μ = 1. If gsp = true, λ = 1, μ = 0 and the associated generalized saddle point system is solved. λ and μ are also keyword arguments that can be directly modified for more specific problems.

GPMR is based on the orthogonal Hessenberg reduction process and its relations with the block-Arnoldi process. The residual norm ‖rₖ‖ is monotonically decreasing in GPMR.

GPMR stops when itmax iterations are reached or when ‖rₖ‖ ≤ atol + ‖r₀‖ * rtol. atol is an absolute tolerance and rtol is a relative tolerance.

Full reorthogonalization is available with the reorthogonalization option.

Additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations.

GPMR can be warm-started from initial guesses x0 and y0 with the method

(x, y, stats) = gpmr(A, B, b, c, x0, y0; kwargs...)

where kwargs are the same keyword arguments as above.

The callback is called as callback(solver) and should return true if the main loop should terminate, and false otherwise.

Reference

• A. Montoison and D. Orban, GPMR: An Iterative Method for Unsymmetric Partitioned Linear Systems, Cahier du GERAD G-2021-62, GERAD, Montréal, 2021.

solver = gpmr!(solver::GpmrSolver, A, B, b, c; kwargs...)
solver = gpmr!(solver::GpmrSolver, A, B, b, c, x0, y0; kwargs...)

where kwargs are keyword arguments of gpmr.

Note that the memory keyword argument is the only exception. It's required to create a GpmrSolver and can't be changed later.

See GpmrSolver for more details about the solver.