Julia is slower than matlab when it comes to matrix diagonalization?

General Usage
, ,
21

The code is below, reproduced from this post:

using ArnoldiMethod, LinearAlgebra, LinearMaps

PQ = P*Q

# Factorizes A and builds a linear map that applies inv(A) to a vector.
function construct_linear_map(A)
    F = factorize(A - ev*I)
    LinearMap{eltype(A)}((y, x) -> ldiv!(y, F, x), size(A,1), ismutating=true)
end

function rerun_ArnoldiMethod(n=1, PQ=PQ)
    Threads.nthreads() == n || error("Restart Julia with -t$n")
    BLAS.set_num_threads(n)
    println("ArnoldiMethod $n Threads:")
    @time begin
        decomp, history = partialschur(construct_linear_map(PQ), nev=NMODES, tol=1e-14, restarts=300, which=LM())
        D2, Exy = partialeigen(decomp)
    end
    println(history)
    D2 = inv.(D2) .+ ev
    @show D2
end

I’m sorry but I don’t have the time to dig up and test the code to recreate the matrix right now. As stated in my original post in that thread, the matrix PQ is a sparse, nonsymmetric, Float64 matrix of order 65274 with 1103319 stored entries.
ev is a scalar with value -12.25.

1 Like