How to extract intermediate results from GMRES of Krylov.jl

Specific Domains Optimization (Mathematical)
4

Hi @jinml!

FOM, GMRES and GPMR are the only methods in Krylov.jl that compute the solution x at the last iteration.
We can’t update x because x = V * y where y is the solution of a triangular system Ry = z. At each iteration, the coefficients of y are different if we perform the backward substitution.

However, you can create a callback that compute each xk. The initial goal of a callback is to add new stopping conditions but you can also use it to see the current state of the KrylovSolver (workspace).

T = eltype(b)
global X = Vector{T}[]

function jinml_callback(solver)
   z = solver.z
   k = solver.inner_iter
   nr = sum(1:k)
   V = solver.V
   R = solver.R
    y = copy(z)
    for i = k : -1 : 1
      pos = nr + i - k               # position of rᵢ.ₖ
      for j = k : -1 : i+1
        y[i] = y[i] - R[pos] * y[j]  # yᵢ ← yᵢ - rᵢⱼyⱼ
        pos = pos - j + 1            # position of rᵢ.ⱼ₋₁
      end
      y[i] = y[i] / R[pos]  # yᵢ ← yᵢ / rᵢᵢ
    end
    xk = sum(V[i] * y[i] for i = 1:k)
    push!(X, xk)
    return false
end

x, stats = gmres(A, b, itmax=3000, memory=500, history=true, callback=jinml_callback)
stats.residuals
residuals2 = [norm(b - A * X[i]) for i = 1:stats.niter]
4 Likes