Complex Exponential of Hermitian

expim-tests.jl

using LinearAlgebra, BenchmarkTools
using LinearAlgebra: RealHermSymComplexHerm
import Base: cis

cis(H::RealHermSymComplexHerm) = cis!(Matrix{complex(eltype(H))}(H),copy(H))
cis!(H::RealHermSymComplexHerm) = cis!(Matrix{complex(eltype(H))}(H),H)
function cis!(R::Matrix{<:Complex},H::RealHermSymComplexHerm)
    # First decompose H into U*Λ*Uᴴ
    Λ,U = eigen!(H)
    # Calculate the imaginary exponential of each eigenvalue and multiply
    # each column of U to obtain B = U*exp(iΛ)
    B = similar(R)
    @inbounds for j = 1:length(Λ)
        a = cis(Λ[j])
        @simd for i = 1:length(Λ)
            B[i,j] = a * U[i,j]
        end
    end
    # Finally multiply B by Uᴴ to obtain U*exp(iΛ)*Uᴴ = exp(i*H)
    mul!(R,B,adjoint(U))
    return R
end

H = Hermitian(rand(ComplexF64,100,100))
exp(1im*H) ≈ cis(H)
@benchmark exp(1im*$H)
@benchmark cis($H)

sizes = [10,20,50,64,100,200,300,512,1000,2000,4096]
res = Array{Float64}(undef,length(sizes),2)
for (i,N) in enumerate(sizes)
    H = Hermitian(rand(ComplexF64,N,N))
    res[i,1] = @belapsed exp(1im*$H)
    res[i,2] = @belapsed cis($H)
end