Large Sparse Hermitian Matrix - Intermediate Eigenvalue Problem

General Usage
, , ,
7

I have given a bit of time to testing a few of the different ideas you have mentioned.

  1. ArnoldiMethod.jl (as elucidated in the documentation) provides similar performance to ARPACK.jl, in my testing it seems ArnoldiMethod is a little bit faster.

  2. I tested out the shift invert method you provided (thank you for that by the way). While it worked well, I ran into memory issues quite quickly.

  3. Using the Pardiso.jl library I was able to get 10 eigenvalues in ~10 mins from a 500,000x500,000 sparse matrix test case. I seem to now be running into memory issues again when I increase the size of my system further (now ~4,000,000x4,000,000). Max memory required by Pardiso is 440GB (HPC node only has 400GB available).

A stripped back version of my current set up is seen below:

using DelimitedFiles
using LinearAlgebra
using SparseArrays
using LinearMaps
using ArnoldiMethod
using Pardiso

struct ShiftAndInvertPardiso
    ps::MKLPardisoSolver
    A_pardiso::Any
end

"""
    shiftinvert_pardiso_operator(A, σ; matrixtype=:COMPLEX_HERM_INDEF)

Constructs a LinearMap for (A - σI)⁻¹ using Pardiso factorization.
"""
function shiftinvert_pardiso_operator(A::SparseMatrixCSC, σ;
                                      matrixtype=:COMPLEX_HERM_INDEF)

    n = size(A, 1)
    # Init Pardiso solver
    ps = MKLPardisoSolver()
    set_matrixtype!(ps, getfield(Pardiso, matrixtype))
    pardisoinit(ps)
    fix_iparm!(ps, :N)

    # Shifted matrix
    A_shifted = A - σ * I

    # Register and factorize
    A_pardiso = get_matrix(ps, A_shifted, :N)
    btmp = zeros(eltype(A), n)

    set_phase!(ps, Pardiso.ANALYSIS)
    pardiso(ps, A_pardiso, btmp)

    println("Estimated mem (KB): ", get_iparm(ps, 15))
    println("Permanent mem (KB): ", get_iparm(ps, 16))
    println("Peak mem (KB): ", get_iparm(ps, 17))
    flush(stdout)

    set_phase!(ps, Pardiso.NUM_FACT)
    pardiso(ps, A_pardiso, btmp)

    # Return LinearMap + solver handle
    linmap = LinearMap{eltype(A)}(
        ShiftAndInvertPardiso(ps, A_pardiso),
        n; ismutating=true
    )
    return linmap, ps
end

# Apply operator y = (A - σI)⁻¹ x
function (op::ShiftAndInvertPardiso)(y, x)

    set_phase!(op.ps, Pardiso.SOLVE_ITERATIVE_REFINE)
    pardiso(op.ps, y, op.A_pardiso, x)
end


# Start of Solver Routine
dE = 0.2
levels = 8

S=readdlm("matrix.txt");

A=sparse(Int.(S[:,1]), Int.(S[:,2]), complex.(S[:,3], S[:,4]))

# Pardiso Commands
T, ps = shiftinvert_pardiso_operator(A, dE; matrixtype=:COMPLEX_HERM_INDEF)
decomp, history = partialschur(T; which=:LM, nev=levels, tol=1e-8)
history.converged || error("The eigenvalue solver did not converge.")
μ, V = partialeigen(decomp)

# λ of the original problem are λ = 1 / μ
D = 1 ./ μ

# Cleanup
set_phase!(ps, Pardiso.RELEASE_ALL)
pardiso(ps)

Am I better off swapping to an iterative method rather than a direct method like this?