Solving a large non-hermitian generalised eigenvalue problem from a linear stability analysis using SLEPc
I have a generalised matrix problem: $A x = \lambda B x$ from a spectral method on a linear stability analysis problem. My matrix B is diagonal and positive semi-definite. A is non-hermitian and complex.
My problem is essentially that when using the SLEPc generalised eigenvalue solver I get the error "zero pivot in LU factorisation". The rest of the below is details about the problem and things I have tried so far. Thanks for the help!
Details of the problem
The matrix will be at its largest about 48000 by 48000, and I want to find the eigenvalues. The eigenvalues I am interested in are ones with the largest real part near 0+0i. Ideally, I want to be able to find them even if they are internal (i.e when there are other eigenvalues with larger positive real part in the spectrum). However, I would be happy if I could get it to work for problems where all eigenvalues have real parts < 0 apart from the eigenvalue of interest.
At the moment I have used the scipy linalg.eig and sparse.eigs functions. As far as I know, these use LAPACK and ARPACK respectively to do the heavy lifting. I have decided to see if I can achieve better performance through using the SLEPc library. 768
If this is a bad decision, let me know!
I want to move onto using PETSc with the SLEPc eigenvalue solvers. I have been trying out SLEPc using the examples provided as part of the tutorial. Exercise 7 (http://www.grycap.upv.es/slepc/handson/handson3.html) reads matricies A and B from a file and outputs the solutions. I got this to work fine using the matrices provided. However, if I substitute a smaller sized test version of my problem (6000x6000), I get a variety of errors depending on the command line arguments I supply.
The main problem I have is the error: "zero pivot in LU factorisation!" when I use the default settings.
I think this might be related to the fact that B contains rows of zeros, although my understanding of linear algebra is somewhat basic. Is this true?
I have tried setting the options suggested on the petsc website, -pc_factor_shift_type NONZERO etc but all I get is an additional warning that these options were not used
I assumed that this was a problem with the preconditioner, so I tried setting -eps_target to 0.1 and both with and without specifying -st_type sinvert and shift. Still I get the same error.
Then I tried -st_pc_type jacobi and st_pc_type bjacobi. jacobi runs, but does not produce any eigenvalues. Block jacobi does an LU factorisation and gives me the same error again.
The default method is krylov-schur, so I have experimented with the -eps_type gd and -eps_type jd options. Unfortunately these seem to produce nonsense eigenvalues, which do not appear on the spectrum at all when I solve using LAPACK in scipy.
I know my matrix problem is not singular, because I can solve it using scipy.
Do you know of any books/guides I might need to read besides the PETSC and SLEPC manuals to understand the behaviour of all these different solvers?
The output from the case with no command line options is given below.
Thanks a lot for taking the time to read my first post!
Kind Regards, <PERSON>
Terminal output from SLEPc
tobymac:SLEPC toby$ mpiexec ./ex7 -f1 LHS-N7-M40-Re0.0-b0.1-Wi5.0-amp0.02.petsc -f2 RHS-N7-M40-Re0.0-b0.1-Wi5.0-amp0.02.petsc -eps_view
Generalized eigenproblem stored in file.
[0]PETSC ERROR: --------------------- Error Message ------------------------------------ [0]PETSC ERROR: Detected zero pivot in LU factorization: see http://www.mcs.anl.gov/petsc/documentation/faq.html#ZeroPivot! [0]PETSC ERROR: Empty row in matrix: row in original ordering 2395 in permuted ordering 3600! 768