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


  • 3
    $\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 Jan 14 '17 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 Jan 14 '17 at 14:47
  • $\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$ – Brian Borchers Jan 15 '17 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 Jan 15 '17 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$ – Wolfgang Bangerth Jan 16 '17 at 21:44

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.