[SciPy-dev] Ticket #709 (was: Re: Next scipy release 0.7)

[Nils Wagner]

On Tue, 30 Sep 2008 18:36:12 +0000 (UTC)
  Pauli Virtanen <pav at 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 at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev