Unitary propagator speed test

SimSpeedTest.jl

require("profile")
using Profile

@profile begin
function expm_eigen(A::Matrix{Complex128}, t)
    #Calculates expm(t*A) via eigenvalue decomposition and assuming Hermitian matrix
    D, V = eig(A, true)

    # V * diagm(exp(t*D)) * V'
    return diagmm(V, exp(t*D)) * V'
end

function evolution_unitary(Hnat::Matrix{Complex128}, controlHams::Array{Complex128, 3}, controlFields::Matrix{Float64}, controlFreqs::Vector{Float64})

    timeStep = 0.01
    Uprop = eye(Complex128, size(Hnat,1))

    for timect = 1:size(controlFields,2)
        tmpH = copy(Hnat)
        for controlct = 1:size(controlFields,1)
            tmpH += controlFields[controlct, timect]*cos(2*pi*timeStep*timect*controlFreqs[controlct])*controlHams[:, :, controlct]
        end
        Uprop *= expm_eigen(tmpH, -1im*2*pi*timeStep)
    end

    return Uprop
end

function parallel_evolution_unitary(Hnat::Matrix{Complex128}, controlHams::Array{Complex128, 3}, controlFields::Matrix{Float64}, controlFreqs::Vector{Float64})

    timeStep = 0.01

    Uprop = @parallel (*) for timect = 1:size(controlFields,2)
        tmpH = copy(Hnat)
        for controlct = 1:size(controlFields,1)
            tmpH += controlFields[controlct, timect]*cos(2*pi*timeStep*timect*controlFreqs[controlct])*controlHams[:, :, controlct]
        end
        expm_eigen(tmpH, -1im*2*pi*timeStep)
    end

    return Uprop
end
end # end of @profile begin

function sim_setup(dimension::Int64, numTimeSteps::Int64, numControls::Int64)
    #Create a random natural hamiltonian 
    tmpMat = randn(dimension, dimension) + 1im*randn(dimension, dimension)
    Hnat = tmpMat+tmpMat'

    #Create random control Hamiltonians
    controlHams = zeros(Complex128, (dimension, dimension, numControls))
    for ct = 1:numControls
        tmpMat = randn(dimension, dimension) + 1im*randn(dimension, dimension)
        controlHams[:,:,ct] = tmpMat+tmpMat'
    end
    #Create random controlfields
    controlFields = randn(numControls, numTimeSteps)
        
    #Control frequencies
    controlFreqs = randn(numControls)

    return Hnat, controlHams, controlFields, controlFreqs

end

function run_sim()
    evolution_unitary(Hnat, controlHams, controlFields, controlFreqs)
end

function run_parallel_sim()
    parallel_evolution_unitary(Hnat, controlHams, controlFields, controlFreqs)
end

Hnat, controlHams, controlFields, controlFreqs = sim_setup(16, 2000, 4)