Singular exception in eig(A,B)

[ViralBShah]

This seems to work in Octave.

a = rand(10,10)
b = rand(10,10)
eig(a+a',b+b')

In Julia, we get a SingularException(12).

[andreasnoack]

The symmetric solver in LAPACK requires that the second argument is positive definite. The easy fix is to use the general solver, but you'll lose all the usual benefits of a symmetric solver. All the symmetric generalised eigenproblems I meet have positive definite B. Are non-positivedefinite B symmetric generalised eigenproblems common?

[andreasnoack]

But the error message is wrong. It should tell that the problem is definiteness.

[ViralBShah]

I discovered this while testing eigs. The eigs symmetric solver also has a similar requirement. In both cases, we should perhaps test for this and give the right error.

[jiahao]

This is an edge case I missed in the error handling. At the very least, the default error handler in LAPACK.geevx! should be replaced by one that understands an error code >N. The complication is that there are only two such error codes, N+1 (failed to compute QZ rotation in gehqz!) and N+2 (failed to compute generalized eigenvector pairs in tgevc!), but each of these error codes could represent multiple failure modes in the underlying routines.

It also looks like for error codes between 1 and N that the computation is partially successful. Not sure if it's worth trapping and returning the partial computation.

[andreasnoack]

Unless we want to be able to handle symmetric, but non-positive definite generalised eigenproblems, I think it is enough just to use @assertposdef instead of @assertnonsingular. This will give an informative error message for @ViralBShah's example. It will result in wrong error messages for the more exotic errors of (sy/he)gvd, but I think they are quite rare.

If we want to cover problems similar to @ViralBShah's, I think we should return the error code from (sy/he)gvd if it is positive and then fall back to the general solver when info[1]>n.

[stevengj]

@andreasnoackjensen, I know that indefinite real-symmetric generalized eigenproblems appear in wave-propagation problems (it is the form of the eigenproblem for the modes of a waveguide at a given temporal frequency ω, where the eigenvalue is the spatial frequency or "propagation constant" β). The fact that it is indefinite is important because it means that complex β solutions (evanescent modes) can arise, while the symmetry is related to electromagnetic reciprocity and mode orthogonality.

[andreasnoack]

@stevengj Thank you for the example.

For indefinite B problems, the eigfact now defaults to the general solver so this problem should be fixed.