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

Pauli Virtanen pav@iki...
Tue Sep 30 14:43:10 CDT 2008


Tue, 30 Sep 2008 21:21:38 +0200, Nils Wagner wrote:

> 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 ?

I didn't contact LAPACK devs yet. I'll try to see if I can reduce this to 
a simple Fortran testcase.

But to my understanding, LAPACK's DGEEV has also some other known issues, 
see eg. http://icl.cs.utk.edu/lapack-forum/viewtopic.php?
f=2&t=317&p=1060&hilit=dggev

> Can someone run the test (2dof.py) in Matlab, Octave, Scilab ?

I got similar suspect result from Octave's QZ when I tried it (Octave 
dropped one of the eigenvalues, too).

Matlab AFAIK does not use Lapack's DGGEV, so the testcase probably will 
not fail in the same way with it.

>> 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 ?

For complex matrices the eigenvectors are returned by LAPACK as complex 
data directly (and don't need to be separately constructed from real 
data, as in the real case), so they are simply passed throught by Scipy. 
If the results are wrong, we probably can't do anything to fix or detect 
this.

-- 
Pauli Virtanen



More information about the Scipy-dev mailing list