lapack gives different answers in scipy and slepc4py

[Francis Poulin]

We are building some software to compute the stability of ocean currents. This requires discretizing a system of Partial Differential Equations, setting them up as an eigenvalue problem (sparse) and then computing the eigenvalues.

We have had some success with one model but are having problems with our second model. We have build a sparse matrix in numpy, as someone previously suggested, and we are able to compute the eigenvalues using scipy with eig (lapack) but not eigs (arpack). That's not entirely true. We need to find the eigenvalues with the largest imaginary part and it seems that eigs is happen to give us the eigenvalues with the largest magnitude instead. This is confusing but not something I want to ask about here but it's related to what I do want to ask about.

After building the matrix we convert it to a petsc matrix. When we try using the default krylov based method in slepc4py and it fails or gives us zero eigenvalues. We have tried other methods but either they are not installed on our version of slepc or they fail. The one method that does not fail is lapack (which is not idea but thought we'd try it for testing) and it gives answers that are different from what we get in scipy. To make sure the conversion is correct we compute the eigenvalues using eig, convert to petsc, compute eigenvalues using slepc (lapack) and then convert back and compute the eigenvalues using eig. The two calculations with eig yield the same results, which makes us think that petsc has the right matrix.

My two questions are the following:

1) Why does LAPACK in slepc not yield the same eigenvalues as lapack in scipy? Might it have something to do with balancing the matrix and if yes is this easily done in slepc4py?

2) Why is it that the default krylov method fails is do you have any suggestions how we can get convergence?

Thanks in advance.

[Jose E. Roman]

1) The LAPACK solver in SLEPc currently uses DGEHRD+DHSEQR, instead of a driver subroutine. Hence, as you pointed out, it is not balancing the matrix. The LAPACK solver in SLEPc is intended for debugging purposes, it is not supposed to be used in production runs.

2) This question is too generic to be answered without further information. However, since balancing seems to be important in your application, you might want to try setBalance() to see if that helps.

[Francis Poulin]

Thank you for your response.

1) We tried using balancing with LAPACK and it seems to work better but not as robust. We don't plan to use this regularly, it is just for testing as you mentioned.

2) Sorry for my vagueness. Our matrix is real non-Hermitian and we need to find the eigenvalues with the largest imaginary part, which has one of the smallest magnitudes. One issue is that when we use a krylov method to find these eigenvalues it finds other eigenvalues, which have a larger magnitude. I have been told that these are easier to converge to. So I suppose the problem reduces down to finding particular eigenvalues. To do this should we try using shift and invert?

[Jose E. Roman]

The scenario you describe is difficult. You can try shift-and-invert if you know a target value around which the wanted eigenvalues are. Note that using a complex target is possible only if SLEPc has been built with complex scalars. If this strategy is not good for you, contact my personal email and we'll see if there are alternatives.

[Francis Poulin]

We have been trying the shift-and-invert and it does help. Sometimes it converges to the correct values but sometimes not. Unfortunately, our SLEPc was not built with complex scalars. We will look into rebuilding it with the right options. However, I imagine there is not an easy fix to this but I might write to you for some guidance along the way.

Thanks for the advice.