Source

geoalg / geoalg / gauss.py

#!/usr/bin/env python
#coding:utf-8
# Author:   --<>
# Purpose:
# Created: 27.03.2010

from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import absolute_import

__author__ = 'mozman <mozman@gmx.at>'

__all__ = ["gaussian_elimination"]

class GaussianElimination(object):
    def __init__(self, matrix):
        self.matrix = matrix
        self.rank = len(matrix)
        self.x = [0.0] * self.rank

    def run(self):
        self.eliminate()
        self.substitute()

    def eliminate(self):
        a = self.matrix
        for i in range(self.rank):
            max = i
            for j in range( i+1, self.rank):
                if abs(a[j][i]) > abs(a[max][i]):
                    max = j
                tmp = a[i]
                a[i] = a[max]
                a[max] = tmp
            for j in range(i+1, self.rank):
                tmp = a[j][i] / a[i][i]
                for k in reversed(range(self.rank)):
                    a[j][k] -= a[i][k] * tmp

    def substitute(self):
        a = self.matrix
        x = self.x
        for j in reversed(range(self.rank)):
            tmp = 0.0
            for k in range(j+1, self.rank):
                tmp += a[j][k] * x[k]
            x[j] = (a[j][self.rank+1] - tmp) / a[j][j]

def gaussian_elimination(matrix):
    solver = GaussianElimination(matrix)
    solver.run()
    return solver.x