수치해석 레포트로 나온 TriDiagonal 문제에 대한 나름대로의 (--;) TFP 식 접근 풀이. 오히려 다른 사람들보다 소스가 커지긴 했지만, 소스를 원한다면 더 줄일 수 있고 (단, 코드를 알아보기 어려워질까봐 묶을 수 있는 부분에 대해서도 풀어 씀), LU 분해 부분에 대해서는 어느정도 일반화를 시켰다고 생각하기에 그리 나쁘진 않다고 생각하는중.

LU 분해 뒤 해당 계산이 제대로 되었는지를 확인하기 위해 Numeric Python 모듈을 이용했다. [http://sourceforge.net/project/showfiles.php?group_id=1369 여기] 에서 받을 수 있다.

 * 느낀점 - LU 분해를 일반화시키는 부분에 대해서는 연습장에서 계산을 시키고 대입을 해보고 패턴을 찾아낸 뒤 일반화시켰다. 만일 이 부분을 코드를 짜면서 ["TestFirstProgramming"] * ["Refactoring"] 의 방법으로 풀었더라면 어떠했을까. (두 개의 과정이 비슷하다고 여겨진다. 코드가 줄어들고 OAOO 의 원칙을 지킬 수 있을테니까.)  -- 석천

 * 해결한 시간 : 대략 6시간 10분 (수학 일반화 계산 1시간 30분, LU 분해 2시간 30분, tri-diagonal 2시간 12분)
=== test_lu.py - LU Decomposition 을 위한 test code ===
{{{~cpp 
import unittest
from Numeric import *
from Matrix import *
from LuDecomposition import *

class TestLuDecomposition(unittest.TestCase):
	def setUp(self):
		self.a = [[1.0,2.0,3.0],[4.0,5.0,6.0],[7.0,8.0,9.0]]
		self.matrixA = Matrix(array(self.a))

	def testLwrite(self):
		ld = LuDecomposition(self.a)
		col = 0	
		ld._writeL(col)
		self.assertEquals(ld._getL(), [[1,0,0],[4,0,0],[7,0,0]])

	def testUwrite(self):
		ld = LuDecomposition(self.a)
		col = 0
		ld._writeL(col)
		row = 0
		ld._writeU(row)
		self.assertEquals(ld._getU(), [[1,2.0,3.0],[0,1,0], [0,0,1]])

	def testLuDecMore(self):
		ld = LuDecomposition(self.a)
		ld._writeL(0)
		ld._writeU(0)
		self.assertEquals(ld._getL(), [[1,0,0],[4,0,0],[7,0,0]])
		self.assertEquals(ld._getU(), [[1,2,3],[0,1,0], [0,0,1]])
		ld._writeL(1)
		ld._writeU(1)
		self.assertEquals(ld._getL(), [[1,0,0],[4,-3,0],[7,-6,0]])
		self.assertEquals(ld._getU(), [[1,2,3],[0,1,2], [0,0,1]])
		ld._writeL(2)
		ld._writeU(2)
		self.assertEquals(ld._getL(), [[1,0,0],[4,-3,0],[7,-6,0]])
		self.assertEquals(ld._getU(), [[1,2,3],[0,1,2], [0,0,1]])

	def testLuDecomposition(self):
		l, u = LuDecomposition(self.a).perform()
		matrixL = Matrix(array(l))
		matrixU = Matrix(array(u))
		actual = matrixL * matrixU
		self.assertEquals(str(actual), str(self.matrixA))


if __name__=="__main__":
	unittest.main()
}}}

=== LuDecomposition.py - LU Decomposition ===
{{{~cpp 
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))
}}}

=== test_tridiagonal.py - tri-diagonal 부분에 대한 test-code ===
{{{~cpp 
import unittest
from Numeric import *
from Matrix import *
from LuDecomposition import *
from TriDiagonal import *

class TestTridiagonal(unittest.TestCase):
	def testGetMatrixY(self):
		a = [[1,1,-1],[2,6,-4],[1,-1,-1]]
		b = [[3],[6],[8]]
		l, u = LuDecomposition(a).perform()
		expected = [[3],[0],[-5]]
		actual = getMatrixY(l, b)
		self.assertEquals(actual, expected)
	
	def testGetMatrixX(self):
		a = [[1,1,-1],[2,6,-4],[1,-1,-1]]
		b = [[3],[6],[8]]
		l, u = LuDecomposition(a).perform()
		expectedY = [[3],[0],[-5]]
		matrixY = getMatrixY(l, b)
		self.assertEquals(matrixY, expectedY)
		expectedX = [[0.5],[-2.5],[-5]]
		matrixX = getMatrixX(u, matrixY)
		self.assertEquals(matrixX, expectedX)
		
if __name__=="__main__":
	unittest.main()
}}}

=== tri-diagonal.py - 실질적인 tri-diagonal 을 계산하는 모듈 ===
{{{~cpp 
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)
}}}