ALGLIB 를 써서 행렬을 LU 분해하기

ALGLIB is a cross-platform open source numerical analysis and data processing library. It is written in specially designed pseudocode which is automatically translated into several target programming languages (C++, C# and other). ALGLIB is relatively young project - active development started only in 2008, while GSL, for example, has 14 years long history. However, it is actively developed with new releases every 1–2 months.

참조: http://en.wikipedia.org/wiki/ALGLIB
ALGLIB 홈페이지: http://www.alglib.net/
ALGLIB 다운로드: http://www.alglib.net/download.php


 * ALGLIB 를 이용하여 LU 분해하는 파이썬 언어 소스

#!/usr/bin/env python

"""
 Filename: testLU2.py

    Find a LU decomposition of a given square matrix by using
    the function 'xalglib.rmatrixlu' of the library ALGLIB
    (www.alglib.net), Sergey Bochkanov and Vladimir Bystritsky.

    And find also the inverse of a given invertible square matrix
    by using the function 'xalglib.rmatrixinverse' of the library ALGLIN.

 Liscense: GPLv2
 Author: Copyright (c) 2011/09/19 (Mon)  PKim  pkim01 (AT) paran (DOT) com
"""

import math
import xalglib

class RVector:
    def __init__(self, arr):
        a = []
        for i in range(len(arr)):
            a.append(arr[i])
        self.arr = a
    def __str__(self):
        return str(self.arr)
    def __getitem__(self, key):
        return self.arr[key]
    def __setitem__(self, key, value):
        self.arr[i] = value
    def __add__(self, other):
        m = self.size()
        a = []
        for i in range(m):
           a.append(self[i] + other[i])
        return RVector(a)
    def __neg__(self):
        m = self.size()
        a = []
        for i in range(m):
           a.append(-self[i])
        return RVector(a)
    def __sub__(self, other):
        m = self.size()
        a = []
        for i in range(m):
           a.append(self[i] - other[i])
        return RVector(a)
    def __rmul__(self, other):
        if type(other) == int or type(other) == long or type(other) == float:
            m = self.size()
            a = []
            for i in range(m):
               a.append(other*self[i])
            return RVector(a)
        else:
            return None
    def angle(self, other):
        return math.acos( self.dot(other) / (self.dot(self) * other.dot(other)) )
    def cross(self, other):
        a = []
        value = self[1]*other[2] - self[2]*other[1]
        a.append(value)
        value = self[2]*other[0] - self[0]*other[2]
        a.append(value)
        value = self[0]*other[1] - self[1]*other[0]
        a.append(value)
        return RVector(a)
    def norm(self):
        n = self.size()
        sum = 0
        for i in range(n):
            sum = sum + self[i]*self[i]
        return math.sqrt(sum)
    def dot(self, other):
        n = self.size()
        sum = 0
        for i in range(n):
            sum = sum + self[i]*other[i]
        return sum
    def size(self):
        return len(self.arr)
    def printVector(self, msg):
        if msg != None and len(msg) > 0:
            print msg,
        print self