## Block-GMRES

Note

block_gmres works on GPUs with Julia 1.11.

If you want to use block_gmres on previous Julia versions, you can overload the function Krylov.copy_triangle with the following code:

using KernelAbstractions, Krylov

@kernel function copy_triangle_kernel!(dest, src)
i, j = @index(Global, NTuple)
if j >= i
@inbounds dest[i, j] = src[i, j]
end
end

function Krylov.copy_triangle(Q::AbstractMatrix{FC}, R::AbstractMatrix{FC}, k::Int) where FC <: Krylov.FloatOrComplex
backend = get_backend(Q)
ndrange = (k, k)
copy_triangle_kernel!(backend)(R, Q; ndrange=ndrange)
KernelAbstractions.synchronize(backend)
end
Krylov.block_gmresFunction
(X, stats) = block_gmres(A, b::AbstractMatrix{FC};
memory::Int=20, M=I, N=I, ldiv::Bool=false,
restart::Bool=false, reorthogonalization::Bool=false,
atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0,
timemax::Float64=Inf, verbose::Int=0, history::Bool=false,
callback=solver->false, iostream::IO=kstdout)

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

(X, stats) = block_gmres(A, B, X0::AbstractMatrix; kwargs...)

Block-GMRES can be warm-started from an initial guess X0 where kwargs are the same keyword arguments as above.

Solve the linear system AX = B of size n with p right-hand sides using block-GMRES.

Input arguments

• A: a linear operator that models a matrix of dimension n;
• B: a matrix of size n × p.

Optional argument

• X0: a matrix of size n × p that represents an initial guess of the solution X.

Keyword arguments

• memory: if restart = true, the restarted version block-GMRES(k) is used with k = memory. If restart = false, the parameter memory should be used as a hint of the number of iterations to limit dynamic memory allocations. Additional storage will be allocated if the number of iterations exceeds memory;
• M: linear operator that models a nonsingular matrix of size n used for left preconditioning;
• N: linear operator that models a nonsingular matrix of size n used for right preconditioning;
• ldiv: define whether the preconditioners use ldiv! or mul!;
• restart: restart the method after memory iterations;
• reorthogonalization: reorthogonalize the new matrices of the block-Krylov basis against all previous matrix;
• atol: absolute stopping tolerance based on the residual norm;
• rtol: relative stopping tolerance based on the residual norm;
• itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to 2 * div(n,p);
• timemax: the time limit in seconds;
• verbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;
• history: collect additional statistics on the run such as residual norms;
• callback: function or functor called as callback(solver) that returns true if the block-Krylov method should terminate, and false otherwise;
• iostream: stream to which output is logged.

Output arguments

• X: a dense matrix of size n × p;
• stats: statistics collected on the run in a SimpleStats structure.
source