#!/usr/bin/env python
# coding=gb2312


# 以上的資訊隨自己的需要改動吧
def print_matrix( info, m ): # 輸出矩陣
    i = 0; j = 0; l = len(m)
    print info
    for i in range( 0, len( m ) ):
        for j in range( 0, len( m[i] ) ):
            if( j == l ):
                print ' |',
            print '%6.4f' % m[i][j],
        print
    print


def swap( a, b ):
    t = a; a = b; b = t


def solve( ma, b, n ):
    global m; m = ma # 這裡主要是方便最後矩陣的顯示
    global s; 
    
    i = 0; j = 0; row_pos = 0; col_pos = 0; ik = 0; jk = 0
    mik = 0.0; temp = 0.0
    
    n = len( m )


    # row_pos 變數標記行迴圈, col_pos 變數標記列迴圈


    print_matrix( "一開始 de 矩陣", m )


    while( ( row_pos < n ) and( col_pos < n ) ):
        print "位置：row_pos = %d, col_pos = %d" % (row_pos, col_pos)
        # 選主元
        mik = - 1
        for i in range( row_pos, n ):
            if( abs( m[i][col_pos] ) > mik ):
                mik = abs( m[i][col_pos] )
                ik = i




        if( mik == 0.0 ):
            col_pos = col_pos + 1
            continue


        print_matrix( "選主元", m )


        # 交換兩行
        if( ik != row_pos ):
            for j in range( col_pos, n ):
                swap( m[row_pos][j], m[ik][j] )
                swap( m[row_pos][n], m[ik][n] );     # 區域之外？


        print_matrix( "交換兩行", m )


        try:
            # 消元
            m[row_pos][n] /= m[row_pos][col_pos]
        except ZeroDivisionError:
            # 除零異常 一般在無解或無窮多解的情況下出現……
            return 0;
            


        j = n - 1
        while( j >= col_pos ):
            m[row_pos][j] /= m[row_pos][col_pos]
            j = j - 1


        for i in range( 0, n ):
            if( i == row_pos ):
                continue
            m[i][n] -= m[row_pos][n] * m[i][col_pos]


            j = n - 1
            while( j >= col_pos ):
                m[i][j] -= m[row_pos][j] * m[i][col_pos]
                j = j - 1


        print_matrix( "消元", m )
        row_pos = row_pos + 1; col_pos = col_pos + 1


    for i in range( row_pos, n ):
        if( abs( m[i][n] ) == 0.0 ):
            return 0
    return 1


if __name__ == '__main__':
    matrix = [[2.0,   0.0, - 2.0,   0.0], 
              [0.0,   2.0, - 1.0,   0.0], 
              [0.0,   1.0,   0.0,  10.0]]


    i = 0; j = 0; n = 0
    # 輸出方程組
    print_matrix( "一開始的矩陣", matrix )


    # 求解方程組, 並輸出方程組的可解資訊
    ret = solve( matrix, 0, 0 )
    
    if( ret!= 0 ):
        print "方程組有解\n"
    else:
        print "方 程組無唯一解或無解\n"


    # 輸出方程組及其解
    print_matrix( "方程組及其解", matrix )
    for i in range( 0, len( m ) ):
        print "x[%d] = %6.4f" % (i, m[i][len( m )])


 

>>> ================================ RESTART ================================
>>> 
一開始的矩陣
2.0000 0.0000 -2.0000  | 0.0000
0.0000 2.0000 -1.0000  | 0.0000
0.0000 1.0000 0.0000  | 10.0000

一開始 de 矩陣
2.0000 0.0000 -2.0000  | 0.0000
0.0000 2.0000 -1.0000  | 0.0000
0.0000 1.0000 0.0000  | 10.0000

位置：row_pos = 0, col_pos = 0
選主元
2.0000 0.0000 -2.0000  | 0.0000
0.0000 2.0000 -1.0000  | 0.0000
0.0000 1.0000 0.0000  | 10.0000

交換兩行
2.0000 0.0000 -2.0000  | 0.0000
0.0000 2.0000 -1.0000  | 0.0000
0.0000 1.0000 0.0000  | 10.0000

消元
1.0000 0.0000 -1.0000  | 0.0000
0.0000 2.0000 -1.0000  | 0.0000
0.0000 1.0000 0.0000  | 10.0000

位置：row_pos = 1, col_pos = 1
選主元
1.0000 0.0000 -1.0000  | 0.0000
0.0000 2.0000 -1.0000  | 0.0000
0.0000 1.0000 0.0000  | 10.0000

交換兩行
1.0000 0.0000 -1.0000  | 0.0000
0.0000 2.0000 -1.0000  | 0.0000
0.0000 1.0000 0.0000  | 10.0000

消元
1.0000 0.0000 -1.0000  | 0.0000
0.0000 1.0000 -0.5000  | 0.0000
0.0000 0.0000 0.5000  | 10.0000

位置：row_pos = 2, col_pos = 2
選主元
1.0000 0.0000 -1.0000  | 0.0000
0.0000 1.0000 -0.5000  | 0.0000
0.0000 0.0000 0.5000  | 10.0000

交換兩行
1.0000 0.0000 -1.0000  | 0.0000
0.0000 1.0000 -0.5000  | 0.0000
0.0000 0.0000 0.5000  | 10.0000

消元
1.0000 0.0000 0.0000  | 20.0000
0.0000 1.0000 0.0000  | 10.0000
0.0000 0.0000 1.0000  | 20.0000

方程組有解

方程組及其解
1.0000 0.0000 0.0000  | 20.0000
0.0000 1.0000 0.0000  | 10.0000
0.0000 0.0000 1.0000  | 20.0000

x[0] = 20.0000
x[1] = 10.0000
x[2] = 20.0000