Eigensolver

[yianzeng]

Description

This PR introduces LinearEigenproblem and LinearEigensolver as analogous classes to LinearVariationalProblem and LinearVariationalSolver.

This enables high-level programming of Finite Element Eigenproblems. The API is documented and the existing Eigenproblem demo has been ported to the new API.

Dirichlet Boundary conditions are supported, currently by shifting the consequential nullspace to a user-specified eigenvalue.

[ksagiyam]

Looks good. Left some comments. Thanks!

[UZerbinati]

Why the bc_shift by default, not zero? This would give an Inf eigenvalue for each bc node.

[dham]

You answered your own question. If bc_shift is zero then the boundary eigenvalues are mathematically undefined. Expecting a numerical algorithm to do something sensible in those circumstances is optimistic, to say the least. In particular, the way SLEPc will solve the generalised eigenvalue problem is $M^{-1}Au = \lambda u$. If bc_shift==0 then $M$ has rows which are zero and the solve $M^{-1}(Au)$ is highly likely to fail. Using a finite shift makes the problem non-singular.

[UZerbinati]

Let me be more precise, my understanding is that if the bc_shift is zero only on the mass matrix there are no boundary eigenvalues, since the only $\lambda$ for which the eigenvalue equation can be verified is $\infty$ and $\infty$ is not a real number. In particular, we have a problem (i.e. a singular pencil) if eigenvectors belong to the kernel of $M$ and $A$ (hence why the eigenvector can not be the zero vector when $M$ is the identity).
If a matrix in the pencil of the generalised eigenvalue problem is singular this doesn't mean that the pencil itself is singular.
Indeed SLEPc will solve the generalised eigenvalue problem as $M^{-1}Au=\lambda u$ by default, but my understanding is that this is incorrect if your $M$ is singular (for obvious reasons). In my opinion, the correct thing to do is to use a spectral transformation and my understanding is that the standard choice in this case is the SINVERT spectral transform, i.e. SLEPc will solve the generalised eigenvalue problem as $(A-\sigma M)^{-1}M u = \lambda u$. If we have a zero shift is equivalent to inverting the role of the mass and stiffness matrix, this is the usual hand-waving solution when $M$ is singular (for example saddle point problems) which if we use a spectral transform can be viewed from a more formal point of view.

[dham]

I've made this change, however you still don't get a particularly nice system, because of the 0 eigenvalues in the inverted system. If you want to see this, change the test to not shift the eigenvalues and you'll see that the error in the non-boundary eigenvalues increases by a lot.

[UZerbinati]

This is very specific. One might be interested in using the eigensolver to compute the eigenvalues of a saddle point problem, which usually has a different mass matrix with a 0 block on the diagonal.

[dham]

It's only a default, and one that is correct in the straightforward case. The user can always provide a different M.