Smallest magnitude eigenvalues of the generalized eigenvalue equation for a large sparse matrix

General Usage
, ,
11

Thank you all for the suggestions. I have now successfully figured out how to solve this problem in parallel.

Before I present my program, let me respond to all the suggestions:

  1. @PetrKryslUCSD suggests to use eigs directly with explicittransform=:none. This actually works, although it took me 9 hours to diagonalize 6 eigenvalues and eigenvectors. This is in principle acceptable, but I am hoping for a quicker solution since I have a dozen of similar generalized eigenvalue equations to solve.

  2. @rveltz and @antoine-levitt mentioned preconditioner. Although I read the wikipedia on this topic, I still don’t fully understand how to use and choose preconditioners for this problem.

  3. @antoine-levitt suggested spectrum folding with (H - \mu)^2, but I am solving a H x =\lambda S x type problem and similar tricks will involve S^{-1}, which seems to be very slow.

  4. @antoine-levitt This is indeed a quantum mechanical problem. It is a non-interacting many-particle eigenstate problem in thousands of atoms. Therefore, the Hilbert space is very large. As I said, it is on the order of 10^5. I am looking for the highest occupied state.

My solution is still to use shift-and-invert transformation. The issue is that the default matrix factorization in julia only utilizes one core and is relatively slow in my case (I believe the 9 hours of eig is mostly spent in factorize). Therefore, I choose Pardiso.jl solver, which utilized all 128 cores. The relevant program is

using LinearAlgebra, Serialization
using Pardiso, Arpack, LinearMaps

function construct_linear_map(H, S)
    ps = MKLPardisoSolver()
    set_matrixtype!(ps, Pardiso.COMPLEX_HERM_INDEF)
    pardisoinit(ps)
    fix_iparm!(ps, :N)
    H_pardiso = get_matrix(ps, H, :N)
    b = rand(ComplexF64, size(H, 1))
    set_phase!(ps, Pardiso.ANALYSIS)
    pardiso(ps, H_pardiso, b)
    set_phase!(ps, Pardiso.NUM_FACT)
    pardiso(ps, H_pardiso, b)
    return (
        LinearMap{ComplexF64}(
            (y, x) -> begin
                set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
                pardiso(ps, y, H_pardiso, S * x)
            end,
            size(H, 1);
            ismutating=true
        ),
        ps
    )
end

lm, ps = construct_linear_map(H - μ * S, S)
λs_inv, X = eigs(lm, which=:LM)
set_phase!(ps, Pardiso.RELEASE_ALL)
pardiso(ps)
λs = 1 ./ λs_inv

Now it only takes me roughly 5 minutes to get the eigenvalues and eigenvectors, which is actually much better than what I expected.

2 Likes