Julia translation of MATLAB expm function (matrix exponential) and test code in Julia and MATLAB

expm.jl

using LinearAlgebra
import Base.exp
import LinearAlgebra.exp!

function expmchk()
# EXPMCHK Check the class of input A and
#    initialize M_VALS and THETA accordingly.
    m_vals = [3 5 7 9 13];
    # theta_m for m=1:13.
                theta = [#3.650024139523051e-008
                         #5.317232856892575e-004
                          1.495585217958292e-002  # m_vals = 3
                         #8.536352760102745e-002
                          2.539398330063230e-001  # m_vals = 5
                         #5.414660951208968e-001
               
                         9.504178996162932e-001  # m_vals = 7
                         #1.473163964234804e+000
                          2.097847961257068e+000  # m_vals = 9
                         #2.811644121620263e+000
                         #3.602330066265032e+000
                         #4.458935413036850e+000
                          5.371920351148152e+000];# m_vals = 13
    
    return m_vals, theta
end

function getPadeCoefficients(m)
# GETPADECOEFFICIENTS Coefficients of numerator P of Pade approximant
#    C = GETPADECOEFFICIENTS returns coefficients of numerator
#    of [M/M] Pade approximant, where M = 3,5,7,9,13.
    if (m==3)
        c = [120, 60, 12, 1];
    elseif (m==5)
        c = [30240, 15120, 3360, 420, 30, 1];
    elseif (m==7)
        c = [17297280, 8648640, 1995840, 277200, 25200, 1512, 56, 1];
    elseif (m==9)
        c = [17643225600, 8821612800, 2075673600, 302702400, 30270240, 
             2162160, 110880, 3960, 90, 1];
    elseif (m==13)
        c = [64764752532480000, 32382376266240000, 7771770303897600, 
             1187353796428800,  129060195264000,   10559470521600, 
             670442572800,      33522128640,       1323241920,
             40840800,          960960,            16380,  182,  1];
    end
        
    return c      
end

function PadeApproximantOfDegree(A,m)
#PADEAPPROXIMANTOFDEGREE  Pade approximant to exponential.
#   F = PADEAPPROXIMANTOFDEGREE(M) is the degree M diagonal
#   Pade approximant to EXP(A), where M = 3, 5, 7, 9 or 13.
#   Series are evaluated in decreasing order of powers, which is
#   in approx. increasing order of maximum norms of the terms.

    n = maximum(size(A));
    c = getPadeCoefficients(m);

    # Evaluate Pade approximant.
    if (m == 13)
        # For optimal evaluation need different formula for m >= 12.
        A2 = A*A
        A4 = A2*A2 
        A6 = A2*A4
        U = A * (A6*(c[14]*A6 + c[12]*A4 + c[10]*A2) + c[8]*A6 + c[6]*A4 + c[4]*A2 + c[2]*Matrix{Float64}(I,n,n) )
        V = A6*(c[13]*A6 + c[11]*A4 + c[9]*A2) + c[7]*A6 + c[5]*A4 + c[3]*A2 + c[1]*Matrix{Float64}(I,n,n)
                
        F = (V-U)\(V+U)

else # m == 3, 5, 7, 9
        Apowers = Array{Matrix{Float64}}(undef,ceil(Int,(m+1)/2))
        
        Apowers[1] = Matrix{Float64}(I,n,n)
        Apowers[2] = A*A
        
        for j = 3:ceil(Int,(m+1)/2)
            Apowers[j] = Apowers[j-1]*Apowers[2]
        end
    
        U = zeros(n,n) 
        V = zeros(n,n)
        
        for j = m+1:-2:2
            U = U + c[j]*Apowers[convert(Int,j/2)]
        end
    
        U = A*U
        
        for j = m:-2:1
            V = V + c[j]*Apowers[convert(Int,(j+1)/2)]
        end
    
        F = (V-U)\(V+U)
    end
    
    return F
end

function expm(A)
# EXPM   Matrix exponential.
#   EXPM(X) is the matrix exponential of X.  EXPM is computed using
#   a scaling and squaring algorithm with a Pade approximation.
#
# Julia implementation closely based on MATLAB code by Nicholas Higham
#
    # Initialization
    m_vals, theta = expmchk()

    normA = norm(A,1)

    if normA <= theta[end]
        # no scaling and squaring is required.
        for i = 1:length(m_vals)
            if normA <= theta[i]
                F = PadeApproximantOfDegree(A,m_vals[i])
                break
            end
        end
    else
        t,s = frexp(normA/theta[end])
        s = s - (t == 0.5) # adjust s if normA/theta(end) is a power of 2.
        A = A/2^s          # Scaling
        F = PadeApproximantOfDegree(A,m_vals[end])
        
        for i = 1:s
            F = F*F   # Squaring
        end
    end

    return F
end # expm
exp(x::Matrix{T}) where T<:Real = expm(x)
function exp!(x::Matrix{T}) where T<:Real
    x[:] = expm(x)
end

test_expm.jl

L=[1 2 1 2 3 4 5 6 7 8 16 32 64 128 160 192 256 384 512]

t=zeros(size(L))

for j=1:length(L) 
    tic()
    for i=1:100 
        A=rand(L[j],L[j])
        expm(A)
    end
    print("N = ",L[j]," ")
    t[j]=toc()
end

test_expm.m

L=[1 2 1 2 3 4 5 6 7 8 16 32 64 128 160 192 256 384 512];

t=zeros(size(L));

for j=1:length(L) 
    tic();
    for i=1:100 
        A=rand(L(j),L(j)); 
        expm(A); 
    end
    t(j)=toc();
    disp(['N = ' num2str(L(j)) ' elapsed time: ' num2str(t(j)) ' seconds']);
end

testexpm.jl

using Test
using LinearAlgebra

@testset "Matrix exponential" begin
    @testset "Tests for $elty" for elty in (Float32, Float64, ComplexF32, ComplexF64,Real)
        A1  = convert(Matrix{elty}, [4 2 0; 1 4 1; 1 1 4])
        eA1 = convert(Matrix{elty}, [147.866622446369 127.781085523181  127.781085523182;
                                     183.765138646367 183.765138646366  163.679601723179;
                                      71.797032399996  91.8825693231832 111.968106246371]')
        @test exp(A1) ≈ eA1
        A2  = convert(Matrix{elty},
                      [29.87942128909879    0.7815750847907159 -2.289519314033932;
                       0.7815750847907159 25.72656945571064    8.680737820540137;
                       -2.289519314033932   8.680737820540137  34.39400925519054])
        eA2 = convert(Matrix{elty},
                      [  5496313853692458.0 -18231880972009236.0 -30475770808580460.0;
                       -18231880972009252.0  60605228702221920.0 101291842930249760.0;
                       -30475770808580480.0 101291842930249728.0 169294411240851968.0])
        @test exp(A2) ≈ eA2
        A3  = convert(Matrix{elty}, [-131 19 18;-390 56 54;-387 57 52])
        eA3 = convert(Matrix{elty}, [-1.50964415879218 -5.6325707998812  -4.934938326092;
                                      0.367879439109187 1.47151775849686  1.10363831732856;
                                      0.135335281175235 0.406005843524598 0.541341126763207]')
        @test exp(A3) ≈ eA3

        A4 = convert(Matrix{elty}, [0.25 0.25; 0 0])
        eA4 = convert(Matrix{elty}, [1.2840254166877416 0.2840254166877415; 0 1])
        @test exp(A4) ≈ eA4

        A5 = convert(Matrix{elty}, [0 0.02; 0 0])
        eA5 = convert(Matrix{elty}, [1 0.02; 0 1])
        @test exp(A5) ≈ eA5
        # Hessenberg
        @test hessenberg(A1).H ≈ convert(Matrix{elty},
                                                 [4.000000000000000  -1.414213562373094  -1.414213562373095
                                                  -1.414213562373095   4.999999999999996  -0.000000000000000
                                                  0  -0.000000000000002   3.000000000000000])
    end
end