Solution of singular linear systems.

Question

Hello,

I am solving a linear system for a fluid dynamics problem of which I need to extract even some eigensolutions. I used solve(A, x, b, "lu") in c++ to solve it and I got the right solution. When I invoked SLEPcEigensolver with shift-and-invert option I received PETSc error messages. It complained about zero pivots in LU factorization.

The matrix could be singular, but in this case I can't understand how Fenics solved the system, maybe a different kind of pivoting. I could not find the objects where LU factorization and system solution is actually implemented in Fenics, to have a better insight in what is happening.
Thank you a lot.

Answer 1

LU and other techniques quite often work, but not always.
The LU factorization happens outside FEniCS,

There are several different LU solvers and libraries that can be used
from dolfin, see:
https://bitbucket.org/fenics-project/dolfin/src/12a8c744503fb50a51e7444bd23d22bb0de4ef7e/dolfin/la/?at=master

Comments

Hello
Thank you very much. I read the scripts listed in your link. In particular PETScLusolver.cpp, which I think is used in the direct solution when using problems defined with PETScMatrix. In the end I realized my problem was due to pivoting. I found that there is a flag in PETScLUsolver.cpp on MAT_SOLVER_MUMPS as a method for performing LU factorization. This explains why the program was able to solve the problem but not to extract the eigenvalues. In SLEPcEigensolver.cpp there is no reference to such methods for performing factorization, or maybe I was not able to find them.

I found easier to do it directly from SLEPc. After creating the eigensolver I accessed the inthernal PETSc's KSP object, the linear system solver, and then set the parameters of its preconditioner, PC object, with MAT_SOLVER_MUMPS. Maybe all these procedures could be done directly using Dolfin's methods.
Thank you
Stefano

---

I am having the same problem with eigenvalue solver. Is it possible to set the PC object in python ? Can you show how to do that ?

---

Hello
Unfortunately I'm not using python interface and I'm not expert of it. I think it is quite difficult even using the lower level c++ procedures. I was able to do it just using directly SLEPc. Sorry.
Hello
Stefano