On Tue, 30 Sep 2008 18:36:12 +0000 (UTC) Pauli Virtanen <pav@iki.fi> wrote:
Tue, 30 Sep 2008 19:56:58 +0200, Nils Wagner wrote:
IMHO,
http://projects.scipy.org/scipy/scipy/ticket/709
should be resolved, too.
Short summary about this one: it's about a special matrix pair for which LAPACK's generalized eigenvalue problem solver DGGEV returns bogus output, losing one eigenvalue from a complex eigenpair. (The eigenvalue problem should nonetheless be well-defined, AFAIK, this looks like a bug in LAPACK.)
In that case the LAPACK developers should be informed as soon as possible. Did you check a FORTRAN implementation of the test case ? Can someone run the test (2dof.py) in Matlab, Octave, Scilab ?
The question here is what we should/can do:
1. Raise LinAlgError if we detect this condition.
2. Try to repair the returned data by filling in the other eigenvalue of the pair, as we know all complex eigenvalues come in pairs.
Well, this holds for real matrices but what could be done in case of complex matrices ? Nils
3. Try to return as much as we can make sense of, but put NaN in place of the missing eigenvalue.
Details:
Instead of expected complex eigenvalue pair (a_R +- i a_I)/beta, the LAPACK routine instead returns the values
a_R = whatever a_I = [0., -2.18645248e-16j] beta = [0., 7.95804395e-17]
The correct result would have been something like
a_R = whatever a_I = [+2.18645248e-16j, -2.18645248e-16j] beta = [7.95804395e-17, 7.95804395e-17]
but one eigenvalue of the pair was lost, apparently due to rounding etc.
Eigenvector data corresponding to zeros in a_I are interpreted differently than for positive entries, so the zero in the wrong place messed up the Scipy routine. We could detect this and raise an error (1), substitute in the missing eigenvalue based on the correctly calculated one (2), or leave the missing eigenvalue as corrupted and return only the correct one (3).
-- Pauli Virtanen
_______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev