[SciPy-dev] Bug in linalg.eig

[Nils Wagner]

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()