ガウス消去法をPythonで書いてみた ref: http://qiita.com/IshitaTakeshi/items/a89c700769148f1d4749

file1.txt

[       1       2       0      -1]
[       3       7      -1       3]
[       2       2       1       2]

[       1       2       0      -1]
[       0       1      -1       6]
[       2       2       1       2]

[       1       2       0      -1]
[       0       1      -1       6]
[       0      -2       1       4]

[       1       2       0      -1]
[       0       1      -1       6]
[       0       0      -1      16]

[       1       2       0      -1]
[       0       1       0     -10]
[       0       0      -1      16]

[       1       2       0      -1]
[       0       1       0     -10]
[       0       0      -1      16]

[       1       0       0      19]
[       0       1       0     -10]
[       0       0      -1      16]

[       1       0       0      19]
[       0       1       0     -10]
[       0       0      -1      16]

[       1       0       0      19]
[       0       1       0     -10]
[       0       0      -1      16]

[       1       0       0      19]
[       0       1       0     -10]
[       0       0       1     -16]

gaussian_elimination.py

from fractions import Fraction
import sys
import decimal

import numpy as np


def print_matrix(x):
    for i in range(len(x)):
        sys.stdout.write("[")
        for j in range(len(x[i])):
            m = Fraction(x[i][j]).limit_denominator()
            sys.stdout.write("{!s:>8}".format(m))
        sys.stdout.write("]\n")
    sys.stdout.write("\n")


def gaussian_elimination(x):
    x = np.float64(x)
    N = x.shape[0]

    print_matrix(x)

    for i in range(N - 1):
        for j in range(i + 1, N):
            if(x[i][i] == 0):
                continue
            r = x[j][i] / x[i][i]
            x[j] = x[j] - r * x[i]

            print_matrix(x)

    for i in reversed(range(N)):
        for j in reversed(range(i)):
            if(x[i][i] == 0):
                continue
            r = x[j][i] / x[i][i]
            x[j] = x[j] - r * x[i]

            print_matrix(x)

    for i in range(N):
        if(x[i][i] == 0):
            continue
        x[i] = x[i] / x[i][i]

        print_matrix(x)

    return x

if(__name__ == '__main__'):
    x = [[1,  2,  0, -1],
         [3,  7, -1,  3],
         [2,  2,  1,  2]]

    x = gaussian_elimination(x)