[SciPy-dev] Bug in linalg.eig

Nils Wagner nwagner at iam.uni-stuttgart.de
Fri Mar 24 01:30:19 CST 2006


Norbert Nemec wrote:
> Could you please try to isolate the problem into a standalone script?
> Without knowing the content of the files G.mtx and H.mtx, it is a bit
> hard to understand what G and H are and what the problem actually is.
>
>
>   
Hi Norbert,

I have submitted a bug report.
http://projects.scipy.org/scipy/scipy/ticket/41
You will find a file mat.tar.gz that contains useful information.

Please let me know if you can reproduce the bug.

Thanks in advance
                                  Nils

> Nils Wagner wrote:
>
>   
>> Hi all,
>>
>> It turns out that linalg.eig is erroneous !
>>
>> Please find attached my latest findings ...
>>
>> Unfortunately, I have no idea how to resolve this problem within linalg.eig.
>>
>>
>> As a workaround one can use the functions
>>
>> linalg.flapack.zgegv
>> linalg.flapack.zggev
>>
>> Nils
>>
>>     
> >from scipy import *
> >from pylab import *
>   
>> G = io.mmread('G.mtx')
>> H = io.mmread('H.mtx')
>> #w = linalg.eigvals(G,H)
>> n = G.shape[0]
>> w = zeros(n,Complex)
>> w1= zeros(n,Complex)
>> #
>> # http://www.netlib.org/lapack/patch/src/zggev.f
>> #
>> alpha,beta,vl,vr,work,info =
>> linalg.flapack.zggev(G,H,compute_vl=0,compute_vr=0,lwork=8*n,overwrite_a=0,overwrite_b=0)
>> #
>> # http://www.netlib.org/lapack/patch/src/zgegv.f
>> #
>> alpha1,beta1,vl1,vr1,info1 =
>> linalg.flapack.zgegv(G,H,compute_vl=0,compute_vr=0,lwork=2*n,overwrite_a=0,overwrite_b=0)
>>
>>
>> for i in arange(0,n):
>>    if beta[i] <> 0.0:
>>       w[i] = alpha[i]/beta[i]
>>       w1[i] = alpha1[i]/beta1[i]
>>
>> #w0 = linalg.eigvals(G,H)  # Doesn't work
>>
>> subplot(211)
>> plot(w.real,w.imag,'b.')
>> subplot(212)
>> plot(w1.real,w1.imag,'r.')
>> show()
>>
>> _______________________________________________
>> Scipy-dev mailing list
>> Scipy-dev at scipy.net
>> http://www.scipy.net/mailman/listinfo/scipy-dev
>>
>>
>>  
>>
>>     
>
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-dev
>   




More information about the Scipy-dev mailing list