from LuDecomposition import *
from Numeric import *
from ArrayPrinter import *
import pdb
def getMatrixY(aMatrixLower, b):
matrixY = makeEmptyMatrix(len(b),1)
for n in range(len(b)):
matrixY[n][0] = float(b[n][0] - _minusForGetMatrixY(n, aMatrixLower, matrixY)) / float(aMatrixLower[n][n])
return matrixY
def getMatrixX(aMatrixUpper, y):
limitedMatrix = len(y)
matrixX = makeEmptyMatrix(limitedMatrix,1)
for n in range(limitedMatrix-1, -1,-1):
matrixX[n][0] = float(y[n][0] - _minusForGetMatrixX(n, aMatrixUpper, matrixX)) / float(aMatrixUpper[n][n])
#print "x[%d]: y[%d][0] - minus:[%f] / u[%d][%d]:%f : %f"% (n,n,_minusForGetMatrixX(n, aMatrixUpper, matrixX),n,n, aMatrixUpper[n][n], matrixX[n][0])
return matrixX
def _minusForGetMatrixX(n, aUpperMatrix, aMatrixX):
totalMinus = 0
limitedMatrix = len(aMatrixX)
for t in range(n+1,limitedMatrix):
totalMinus += aUpperMatrix[n][t] * aMatrixX[t][0]
return totalMinus
def _minusForGetMatrixY(n, aLowerMatrix, aMatrixY):
totalMinus = 0
for t in range(n):
totalMinus += aLowerMatrix[n][t] * aMatrixY[t][0]
return totalMinus
def makeEmptyMatrix(aRow, aCol):
emptyMatrix = []
for i in range(0,aRow):
emptyRow = []
for j in range(0,aCol):
emptyRow.append(0.0)
emptyMatrix.append(emptyRow)
return emptyMatrix
def prettyPrintMatrix(aMatrix):
print array2string(array(aMatrix))
def printSeparator():
print "*-" * 35
if __name__=="__main__":
a = [[4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0]]
b = [[5.0], [6.0], [6.0], [6.0], [6.0], [6.0], [6.0], [6.0], [5.0]]
l, u = LuDecomposition(a).perform()
printSeparator()
print "LuDecomposition - Matrix L:"
prettyPrintMatrix(l)
#pdb.set_trace()
printSeparator()
matrixY = getMatrixY(l, b)
print "Calculated - Matrix Y"
prettyPrintMatrix(matrixY)
printSeparator()
print "LuDecomposition - Matrix U:"
prettyPrintMatrix(u)
printSeparator()
print "Final - Matrix X"
matrixX = getMatrixX(u, matrixY)
prettyPrintMatrix(matrixX)
