It can be used to solve linear equation systems or to invert a matrix.高斯消元法用于解決線性代數求多元方程組的解，或者用于求可逆矩陣的逆
嘿嘿，python代碼現成可用：
def gauss_jordan(m, eps = 1.0/(10**10)):
"""Puts given matrix (2D array) into the Reduced Row Echelon Form.
Returns True if successful, False if 'm' is singular.
NOTE: make sure all the matrix items support fractions! Int matrix will NOT work!
Written by J. Elonen in April 2005, released into Public Domain"""
(h, w) = (len(m), len(m[0]))
for y in range(0,h):
maxrow = y
for y2 in range(y+1, h): # Find max pivot
if abs(m[y2][y]) > abs(m[maxrow][y]):
maxrow = y2
(m[y], m[maxrow]) = (m[maxrow], m[y])
if abs(m[y][y]) <= eps: # Singular?
return False
for y2 in range(y+1, h): # Eliminate column y
c = m[y2][y] / m[y][y]
for x in range(y, w):
m[y2][x] -= m[y][x] * c
for y in range(h-1, 0-1, -1): # Backsubstitute
c = m[y][y]
for y2 in range(0,y):
for x in range(w-1, y-1, -1):
m[y2][x] -= m[y][x] * m[y2][y] / c
m[y][y] /= c
for x in range(h, w): # Normalize row y
m[y][x] /= c
return True
使用方法 :
If your matrix is of form [A:x] (as is usual when solving systems), items of A and x both have to be divisible by items of A but not the other way around. Thus, you could, for example, use floats for A and vectors for x. Example:
mtx = [[1.0, 1.0, 1.0, Vec3(0.0, 4.0, 2.0), 2.0],
[2.0, 1.0, 1.0, Vec3(1.0, 7.0, 3.0), 3.0],
[1.0, 2.0, 1.0, Vec3(15.0, 2.0, 4.0), 4.0]]
if gauss_jordan(mtx):
print mtx
else:
print "Singular!"
# Prints out (approximately):
#
# [[1.0, 0.0, 0.0, ( 1.0, 3.0, 1.0), 1.0],
# [0.0, 1.0, 0.0, ( 15.0, -2.0, 2.0), 2.0],
# [0.0, 0.0, 1.0, (-16.0, 3.0, -1.0), -1.0]]
Auxiliary functions contributed by Eric Atienza (also released in Public Domain):
def solve(M, b):
"""
solves M*x = b
return vector x so that M*x = b
:param M: a matrix in the form of a list of list
:param b: a vector in the form of a simple list of scalars
"""
m2 = [row[:]+[right] for row,right in zip(M,b) ]
return [row[-1] for row in m2] if gauss_jordan(m2) else None
def inv(M):
"""
return the inv of the matrix M
"""
#clone the matrix and append the identity matrix
# [int(i==j) for j in range_M] is nothing but the i(th row of the identity matrix
m2 = [row[:]+[int(i==j) for j in range(len(M) )] for i,row in enumerate(M) ]
# extract the appended matrix (kind of m2[m:,...]
return [row[len(M[0]):] for row in m2] if gauss_jordan(m2) else None
def zeros( s , zero=0):
"""
return a matrix of size `size`
:param size: a tuple containing dimensions of the matrix
:param zero: the value to use to fill the matrix (by default it's zero )
"""
return [zeros(s[1:] ) for i in range(s[0] ) ] if not len(s) else zero
算法偽代碼:
i := 1j := 1
while (i ≤ m and j ≤ n) do
Find pivot in column j, starting in row i:
maxi := i
for k := i+1 to m do
if abs(A[k,j]) > abs(A[maxi,j]) then
maxi := k
end if
end for
if A[maxi,j] ≠ 0 then
swap rows i and maxi, but do not change the value of i
Now A[i,j] will contain the old value of A[maxi,j].
divide each entry in row i by A[i,j]
Now A[i,j] will have the value 1.
for u := i+1 to m do
subtract A[u,j] * row i from row u
Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
end for
i := i + 1
end if
j := j + 1
end while