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

[Nils Wagner]

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