Block-Krylov processes

[amontoison]

https://jso.dev/Krylov.jl/previews/PR800/block_processes/

• Block Hermitian Lanczos
• Block Non-Hermitian Lanczos
• Block Arnoldi
• Block Golub-Kahan
• Block Saunders-Simon-Yip
• Block Montoison-Orban

[github-actions]

Package name	latest	stable
CaNNOLeS.jl
DCISolver.jl
FletcherPenaltySolver.jl
JSOSolvers.jl
LLSModels.jl
LinearSolve.jl
Percival.jl
RipQP.jl

[codecov]

Codecov Report

Attention: 39 lines in your changes are missing coverage. Please review.

Files	Coverage Δ
src/Krylov.jl	100.00% <ø> (ø)
src/block_krylov_processes.jl	100.00% <100.00%> (ø)
src/krylov_utils.jl	91.17% <85.71%> (-0.20%)	⬇️
src/processes_utils.jl	70.99% <70.99%> (ø)

... and 1 file with indirect coverage changes

📢 Thoughts on this report? Let us know!.

[amontoison]

@frapac
I did a generic function to easily do preliminary tests with block-Krylov methods.
I think that block-minres, block-gmres and block-usymqr are the most relevant methods for our applications.

using Krylov, LinearAlgebra, Printf

function block_krylov(solver::String, A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC
  m, n = size(A)
  t, p = size(B)

  m == n || error("System must be square")
  t == m || error("Inconsistent problem size")

  (itmax == 0) && (itmax = div(n, p) + 1)
  k = 0
  Bnorm = norm(B)
  Brank = rank(B)
  ε = atol + rtol * Bnorm
  solved = Bnorm ≤ ε
  tired = false
  algo = "householder"
  local X

  verbose && @printf("block-%s: system of %d equations in %d variables with %d rhs.\n", solver, m, n, p)
  verbose && @printf("%5s  %7s  %8s\n", "k", "‖Rₖ‖", "rank(Rₖ)")
  verbose && @printf("%5d  %7.1e  %8d\n", k, Bnorm, Brank)

  while !(solved || tired)
    k += 1

    # Block-Krylov processes
    if solver ∈ ("fom", "gmres")
      V, Γ, H = arnoldi(A, B, k; algo)
    elseif solver ∈ ("symmlq", "cg", "minres")
      V, Ψ₁, T = hermitian_lanczos(A, B, k; algo)
    elseif solver ∈ ("usymlq", "usymqr")
      V, Ψ₁, T, U, Φ₁ᴴ, Tᴴ = saunders_simon_yip(A, B, B, k; algo)
    elseif solver ∈ ("bilq", "bicg", "qmr")
      V, Ψ₁, T, U, Φ₁ᴴ, Tᴴ = nonhermitian_lanczos(A, B, B, k)
    else
      error("The method block-$solver is not supported.")
    end

    # Block-Krylov methods
    if solver == "fom"
      rhs = zeros(FC, p*k, p)
      rhs[1:p,1:p] .= Γ
      Y = H[1:k*p,1:k*p] \ rhs
    elseif solver == "gmres"
      rhs = zeros(FC, p*(k+1), p)
      rhs[1:p,1:p] .= Γ
      Y = H \ rhs
    elseif solver ∈ ("symmlq", "bilq", "usymlq")
      if k == 1
        Y = zeros(FC, p*k, p)
      else
        rhs = zeros(FC, p*(k-1), p)
        rhs[1:p,1:p] .= Ψ₁
        T = Matrix(T)
        Y = T[1:(k-1)*p,1:k*p] \ rhs
      end
    elseif solver ∈ ("cg", "bicg")
      rhs = zeros(FC, p*k, p)
      rhs[1:p,1:p] .= Ψ₁
      Y = T[1:k*p,1:k*p] \ rhs
    elseif solver ∈ ("minres", "qmr", "usymqr")
      rhs = zeros(FC, p*(k+1), p)
      rhs[1:p,1:p] .= Ψ₁
      Y = T \ rhs
    end

    # Compute the approximate solution Xₖ and the residual Rₖ
    if solver ∈ ("usymlq", "usymqr")
      X = U[:,1:k*p] * Y[1:k*p,:]
    else
      X = V[:,1:k*p] * Y[1:k*p,:]
    end
    R = B - A * X
    Rnorm = norm(R)
    Rrank = rank(R)
    verbose && @printf("%5d  %7.1e  %8d\n", k, Rnorm, Rrank)
    tired = k ≥ itmax
    solved = Rnorm ≤ ε 
  end
  return X
end

block_fom(A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC = block_krylov("fom", A, B; atol, rtol, itmax, verbose)
block_gmres(A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC = block_krylov("gmres", A, B; atol, rtol, itmax, verbose)
block_symmlq(A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC = block_krylov("symmlq", A, B; atol, rtol, itmax, verbose)
block_cg(A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC = block_krylov("cg", A, B; atol, rtol, itmax, verbose)
block_minres(A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC = block_krylov("minres", A, B; atol, rtol, itmax, verbose)
block_usymlq(A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC = block_krylov("usymlq", A, B; atol, rtol, itmax, verbose)
block_usymqr(A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC = block_krylov("usymqr", A, B; atol, rtol, itmax, verbose)
block_bilq(A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC = block_krylov("bilq", A, B; atol, rtol, itmax, verbose)
block_bicg(A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC = block_krylov("bicg", A, B; atol, rtol, itmax, verbose)
block_qmr(A, B::AbstractMatrix{FC}; atol::FC=√eps(FC), rtol::FC=√eps(FC), itmax::Int=0, verbose::Bool=false) where FC = block_krylov("qmr", A, B; atol, rtol, itmax, verbose)