# L-20 MCS 471 Fri 8 Oct 2021 : gmres.jl

# Illustrates the Generalized Minimum Residual Method.

using Printf
Base.show(io::IO, f::Float64) = @printf(io, "%.3e", f)

using LinearAlgebra

"""
    gmres(A::Array{Float64,2}, b::Array{Float64,1},
          x::Array{Float64,1}, N::Int=8)

runs N steps of the Generalized Minimum Residual Method,
to solve the system with matrix A with right hand side vector b,
starting at x.

The number N of iterations should be not be larger
than the number of columns of A.
"""
function gmres(A::Array{Float64,2}, b::Array{Float64,1},
               x0::Array{Float64,1}, N::Int, verbose::Bool=true)
    nrows = size(A,1)
    ncols = size(A,2)
    x = zeros(ncols)
    Q = zeros(nrows, ncols+1)
    H = zeros(nrows+1, ncols+1)
    r = b - A*x0[1:ncols]
    q = r/norm(r)
    Q[:,1] = q
    stop = false
    for k=1:N
        y = A*Q[1:ncols,k]
        for j=1:k
            H[j,k] = transpose(Q[:,j])*y
            y = y - H[j,k]*Q[:,j]
        end
        H[k+1,k] = norm(y)
        if H[k+1,k] == 0
            stop = true
        else
            Q[:,k+1] = y/H[k+1,k]
        end
        rhs = zeros(k+1)
        rhs[1] = norm(r)
        c = H[1:k+1,1:k]\rhs  # least squares solving
        if verbose
            println("H at step ", k, " : ")
            show(stdout, "text/plain", H[1:k+1,1:k]); println("")
        end
        dx = Q[:,1:k]*c
        x = x0 + dx
        if verbose
            println("x at step ", k, " : ")
            show(stdout, "text/plain", transpose(x)); println("")
            err = @sprintf("%.2e", norm(dx))
            println("|dx| : ", err);
        end
        if stop
            return x
        end
    end
    return x
end

"""
Runs gmres on a random example.
"""
function main()
    rowdim = 4
    coldim = 3
    A = rand(rowdim, coldim)
    sol = rand(coldim)
    b = A*sol
    x0 = rand(rowdim)
    x = gmres(A,b,x0,3)
    r = b - A*x[1:3]
    err = @sprintf("%.2e", norm(r'*A))
    println("the residual : ", err)
end

main()