GitHub

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_gmres — Function
(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
Krylov.block_gmres! — Function
solver = block_gmres!(solver::BlockGmresSolver, B; kwargs...)
solver = block_gmres!(solver::BlockGmresSolver, B, X0; kwargs...)

where kwargs are keyword arguments of block_gmres.

See BlockGmresSolver for more details about the solver.

source