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

# Illustrates the power method to compute the eigenvector
# with the largest eigenvalue.

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

using LinearAlgebra

"""
    power(A::Array{Float64,2},x0::Array{Float64,1},
          tol::Float64=1.0e-8,maxit::Int=10,verbose::Bool=true)

Runs the power method on A, starting at the vector x0.

Does no more than maxit steps.

Stops when the error is less than the tolerance.

By default, the verbose flag is true.
"""
function power(A::Array{Float64,2},x0::Array{Float64,1},
               tol::Float64=1.0e-4,maxit::Int=10,verbose::Bool=true)
    x = deepcopy(x0)
    (L, previousL, err) = (0.0, 0.0, 0.0)
    for i=1:maxit
        if verbose
            show(stdout, "text/plain", transpose(x)); println("")
        end
        x = A*x
        previousL = L
        L = norm(x)
        x = x/L
        err = abs(previousL - L)
        if verbose
            strL = @sprintf("%.16e", L)
            strE = @sprintf("%.2e", err)
            println("|lambda1| = $strL, error = $strE")
        end
        if err < tol
            return (x, L, err, false)
        end
    end
    return (x, L, err, true)
end

"""
Runs the power method on a random 4-by-4 matrix.
"""
function main()
    A = rand(4, 4)
    x0 = rand(4)
    z, L, err, fail = power(A,x0)
    println("The dominant eigenvector :")
    show(stdout, "text/plain", transpose(z)); println("")
    strL = @sprintf("%.16e", L)
    strE = @sprintf("%.2e", err)
    println("The spectral radius : $strL,")
    print("with error $strE.  ")
    if fail
        println("Failed to reach the tolerance.");
    else
        println("Reached the tolerance.");
    end
end

main()