In advance I am sorry for my noobish question. I am a physics PHD student and basically I use python for my math/physics problems. But now I have a problem which requires more computing capacity and I intend to use Fortran or C and supercomputer. Those things are very new to me :).

So my problem is:

  • I have a system of 17576 atoms and real symmetric weighted hessian m atrix which describes interatomic interactions and has 52728x52728=2780241984 elements.
  • This matrix is very sparse: 1.56% of elements are nonzero
  • The minimum and maximum of eigenvalues and density of states are known
  • I need to diagonalise this matrix and find all eigenvalues and eigenvectors with high precision

I heard that there is FEAST and SLEPc. However SLEPc is not suitable for whole eigenspectrum calculations. Is it ok to use FEAST for this kind of problem or I should use something different? In FEAST you can define search interval and in documentation it is written that high accuracy can be obtained only for up to 1000 eigenpairs. So my idea is to calculate whole eigenspectrum by splitting this problem into multiple intervals with less than 1000 eigenvalues. Or my intuition is wrong about how should I approach this problem?

So in summary my questions are:

  • Is it OK to use FEAST for this kind of problem or I should use something different? Any other recommendations are very acceptable
  • If FEAST is suitable for this problem, how it scales and how should I distribute resources?
  • Also suggestions about what should I read or learn is very appreciated, because I don't know from where to start and I don't have anyone to consult

Thanks for Your attention


  • 7
    $\begingroup$ could I ask why you want to compute all eigen-pairs with high accuracy? this may not be the best idea for (very) large and sparse matrices. if you are willing to settle for an approximation of the exact spectrum, then this old post may potentially be useful. $\endgroup$
    – GoHokies
    Commented Jan 14, 2017 at 10:05
  • $\begingroup$ I am working on a project where we want to find specific to an isolated defect in a crystal vibrational modes and use them for luminescence band calculations. Convergence of luminescence band is slow when we increase the size of a system so we need to tackle extremely big models. We need almost all vibrational modes and frequencies for precise calculations. $\endgroup$
    – Baranas
    Commented Jan 14, 2017 at 14:47
  • 5
    $\begingroup$ Have you considered just throwing LAPACK/BLAS at it and ignoring the sparsity? You're going to need a lot of storage for the eigenvectors anyway. This computation shouldn't be too unbearable on shared memory machine with say 128 gigabytes of RAM and 32 cores. $\endgroup$ Commented Jan 15, 2017 at 5:50
  • $\begingroup$ Thanks Brian in worst case I will try this approach, but it would be nice to find something a little bit cheaper and faster. I just found this article of shift and invert parallel spectral transformation SIPs. It sounds promising $\endgroup$
    – Baranas
    Commented Jan 15, 2017 at 8:35
  • $\begingroup$ When you say that your matrix has $2,780,241,984$ entries, does that include the zeros? In other words, your matrix is really the given $52728\times 52728$? $\endgroup$ Commented Jan 16, 2017 at 21:44

1 Answer 1


To emphasize what the comments have said, if you want to produce the dense 52728x52728 matrix containing the eigenvectors, there is no point to computing using the sparse input. Use a dense input matrix and call (Sca)LAPACK, ELPA, Elemental or other dense eigensolver library.

Diagonalizing a dense 52728x52728 will take a few minutes or less on GPUs, depending on how much hardware you have, according to GPU Acceleration of Large-Scale Full-Frequency GW Calculations. It will take longer with a big CPU system, but probably not more than 4 hours, assuming you use Intel MKL on an Intel server.

Dense symmetric eigensolver performance

  • 2
    $\begingroup$ Thanks, Jeff. It is an old question, and I already solved this issue using ScaLAPACK and SLEPc. With SLEPc, I used the spectral slicing technique (together with the MUMPS package) to iteratively find the whole eigenspectrum as I had issues with available RAM for very large systems. But the idea of using GPUs is very tempting, at least from looking at your table. I will definitely try this. $\endgroup$
    – Baranas
    Commented Sep 26, 2022 at 9:26
  • $\begingroup$ @Baranas, youn can also create a post and describe your own solution. $\endgroup$
    – nicoguaro
    Commented Sep 26, 2022 at 12:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.