# [SciPy-dev] Urgent: Bug in linalg.lu_factor linalg.lu_solve

Nils Wagner nwagner at mecha.uni-stuttgart.de
Fri Feb 20 07:13:25 CST 2004

```Hi all,

I have observed a very strange behaviour within linalg.lu_factor /
linalg.lu_solve.
Any idea or suggestion ?

Nils

from scipy import *

def F(A,l):

tmp=zeros((n,n),Complex)
tmp[:n,:n] = A+l*identity(n)+exp(-l)*identity(n)
return tmp

def Fs(l):

tmp=zeros((n,n),Complex)
tmp[:n,:n] = identity(n)-l*exp(-l)*identity(n)
return tmp

def T(n):

tmp=zeros((n,n),Complex)
for i in arange(0,n):
tmp[i,i] = -2.0
if i < n-1:
tmp[i+1,i] = 1.0
tmp[i,i+1] = 1.0
return tmp

n = 4
A = T(n)

gamma = 2.0+0j
rho = 2.0

ns = 180
#for i in arange(0,ns):
for i in arange(0,2):
l  = gamma+rho*exp(2.*pi*i*1j/ns)
D1 = F(A,l)
Ds = Fs(l)

Dinv = linalg.inv(D1)
lu, piv = linalg.lu_factor(D1,overwrite_a=0)

tr = 0.0
for k in arange(0,n):
rhs = Ds[:,k]
g1  = linalg.lu_solve((lu,piv),rhs,trans=0,overwrite_b=0)
tr  = tr + g1[k]

#
#  tr and trace(dot(Dinv,Ds)) should be the same !!!!!
#
print tr, trace(dot(Dinv,Ds))

```