[SciPy-user] Urgent: Bug in linalg.lu_factor linalg.lu_solve
Nils Wagner
nwagner at mecha.uni-stuttgart.de
Fri Feb 20 08:13:25 EST 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))
More information about the SciPy-User
mailing list