Power method

powm.jl

using LinearAlgebra

mutable struct State{Tv,V<:AbstractVector{Tv}, T}
    x::V
    r::V
    Ax::V
    θ::T
end

function power_method(A::AbstractMatrix{T}) where {T}
    x = randn(complex(T), size(A, 1))
    x ./= norm(x)

    r = similar(x)
    Ax = similar(x)

    state = State(x, r, Ax, zero(complex(T)))

    if run_power_method!(A, state, 0.001, 100)
        return (state.θ, state.x)
    else
        throw("Algorithm did not converge")
    end
end

function run_power_method!(A, state::State{Tv}, tol::Real, maxiter = 100) where {Tv}
    iter = 0

    while iter < maxiter
        mul!(state.Ax, A, state.x)

        # Approx eigenvalue
        state.θ = dot(state.Ax, state.x)

        # Residual (r = Ax - θx)
        copyto!(state.r, state.Ax)
        state.r .-= state.θ * state.x

        # Normalize eigenvector
        copyto!(state.x, state.Ax)
        rmul!(state.x, inv(norm(state.x)))

        if norm(state.r) < tol
            return true
        end

        iter += 1
    end

    return false
end

function example()
    A = rand(10, 10)
    θ, x = power_method(A)
    @show norm(A * x - θ * x)
end