Tensor dynamical system for Z-eigenvector computation

tensor_dynsys_alg.jl

using LinearAlgebra
function tensor_apply(T::Array{Float64,3}, x::Vector{Float64})
    n = length(x)
    y = zeros(Float64, n)
    for k in 1:n; y += T[:, :, k] * x * x[k]; end
    return y
end
function tensor_collapse(T::Array{Float64,3}, x::Vector{Float64})
    n = length(x)
    Y = zeros(Float64, n, n)
    for k in 1:n; Y += T[:, :, k] * x[k]; end
    return Y
end

function dynsys_forward_euler(T::Array{Float64,3},
                              h::Float64, niter::Int64=20)
    function Λ(u::Vector{Float64})  # Derivative                                                                                                      
        F = eigen(tensor_collapse(T, u))
        ind = sortperm(abs.(real(F.values)))[1]
        v = F.vectors[:, ind]
        return sign(v[1]) * v  # sign consistency                                                                                                     
    end
    x = normalize(ones(Float64, size(T, 1)), 1)  # starting point                                                                                     
    λ_hist = [x' * tensor_apply(T, x)]
    for _ = 1:niter
        x += h * (Λ(x) - x)  # forward Euler                                                                                                          
        push!(λ_hist, x' * tensor_apply(T, x))  # Rayleigh quotient                                                                                   
    end
    return (x, λ_hist)  # guess at evec and history of evals                                                                                          
end

This code is now updated for Julia 1.0.