##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
import matplotlib.pylab as mplot
from numpy import *
from pylab import *
from LU_decomposition import *
import pdb
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')

TL = 4
##Lower band
a=frange(9*(TL-1)-1)
a.shape=(TL-1,3,3)
a[2,:] = 0

##Middle Band
b=frange(9*TL-1)
b.shape=(TL,3,3)
b[0,:] = .25*eye(3)
b[1,:] = 5*eye(3)
b[2,:] = .25*eye(3)
b[2,1,2]=2
b[1,2,2]=5
b[3,:] = 5*eye(3)

##Upper Band
c=4.2*frange(9*(TL-1)-1)
c.shape=(TL-1,3,3)
c[0,:]=0

##Right hand side
f=4.9*frange(3*(TL)-1)
f.shape=(TL,3)

##LU Decomposition 
d_f = LU_solve_3x3_1d(0,TL-1,a,b,c,f)

##Check to see if the rhs is correct
print dot(b[0],d_f[0])+dot(c[0],d_f[1])
print dot(a[0],d_f[0])+dot(b[1],d_f[1])+dot(c[1],d_f[2])
print dot(a[1],d_f[1])+dot(b[2],d_f[2])+dot(c[2],d_f[3])
print dot(a[2],d_f[2])+dot(b[3],d_f[3])