1 # coding:utf8
 2 import numpy as np
 3 
 4 def lu(mat):
 5     r,c=np.shape(mat)
 6     s=min(r,c)
 7     for k in range(s):
 8         x=1.0/mat[k][k]  # 将后续除法变成乘法
 9         for i in range(k+1,r):
10             mat[i][k]=mat[i][k]*x  # L[1:][0]*U[0][0]=A[1:][0]；A[0][:]=mat[0][:]
11         for i in range(k+1,r):
12             for j in range(k+1,c):
13                 # U[1][2]*L[1][1]=A[1][2]-U[0][2]*L[1][0]；L[1][1]=1
14                 # L[2][1]*U[1][1]=A[2][1]-U[0][1]*L[2][0]；下一个k时mat[i][j]/mat[k][k](i>j)
15                 mat[i][j]=mat[i][j]-mat[k][j]*mat[i][k]
16     return mat,c
17 
18 def solve(A,b):
19     mat,n=lu(A)  # LU合并
20     print mat  # [[16, 4, 8], [0.25, 4.0, -6.0], [0.5, -1.5, 9.0]]
21     Z= np.zeros(n)  # L*Z=b U*X=Z
22     X= np.zeros(n)
23     Z[0]=b[0]
24     for i in range(1,n):
25         sumup=0
26         for tmp in range(0,i):
27             sumup+=mat[i][tmp]*Z[tmp]
28         Z[i]=(b[i]-sumup)
29     X[n-1]=Z[n-1]/mat[n-1][n-1]
30     for i in range(n-2,-1,-1):
31         sumup=0
32         for tmp in range(i+1,n):
33             sumup+=mat[i][tmp]*X[tmp]
34         X[i]=(Z[i]-sumup)/mat[i][i]
35     return X
36 
37 A=[[16,4,8],[4,5,-4],[8,-4,22]]
38 y=[-4,3,10]
39 print "The result of the fomula is:"+str(solve(A,y))  # [-2.25  4.    2.  ]