#format python
from Numeric import *
from ArrayPrinter import *

class LuDecomposition:
	def __init__(self, aSquareArray):
		assert len(aSquareArray) == len(aSquareArray[0])
		self.array = aSquareArray
		self.n = len(aSquareArray)

		self.l = self.makeEmptySquareMatrix(self.n)
		self.u = self.makeEmptySquareMatrix(self.n)

		for i in range(self.n):
			self.u[i][i] = 1.0

	def makeEmptySquareMatrix(self, n):
		matrix = []
		for i in range(0,n):
			row = []
			for j in range(0,n):
				row.append(0.0)
			matrix.append(row)
		return matrix

	def perform(self):
		for i in range(self.n):
			self._writeL(i)
			self._writeU(i)
		return self.l, self.u

	def _writeL(self, aCol):
		for eachRow in range(aCol, self.n):
			self.l[eachRow][aCol] = self.array[eachRow][aCol] - self._minusForWriteL(eachRow, aCol)

	def _minusForWriteL(self, aRow, aCol):
		totalMinus = 0
		for n in range(0, aCol):
			totalMinus += self.l[aRow][n] * self.u[n][aCol]
		return totalMinus

	def _writeU(self, aRow):
		for eachCol in range(aRow+1, self.n):
			self.u[aRow][eachCol] = float(self.array[aRow][eachCol] - self._minusForWriteU(eachCol, aRow)) / float(self.l[aRow][aRow]) 

	def _minusForWriteU(self, aCol, aRow):
		totalMinus = 0
		for n in range(0, aRow):
			totalMinus += self.l[aRow][n] * self.u[n][aCol]
		return totalMinus

	def _getL(self): 
		return self.l

	def _getU(self):
		return self.u

if __name__=="__main__":
		a = [[1.0,2.0,3.0],[4.0,5.0,6.0],[7.0,8.0,9.0]]
		l, u = LuDecomposition(a).perform()
		print "array : "
		print array2string(array(a))
		print "Lu Decomposition - L" 
		print array2string(array(l))
		print "Lu Decomposition - U"
		print array2string(array(u))