# Gaussian Elimination with row pivoting¶

In :
import numpy as np


This function performs LU decomposition of given matrix A with pivoting to obtain $$PA = LU$$ where $P$ is a permutation matrix. Here we dont store $P$ as a matrix, but only as a vector.

In :
def LU(A):
n = A.shape
L = np.identity(n)
P = np.arange(n,dtype=int) # Permutation matrix
U = np.array(A)
for k in range(n-1):
i = np.argmax(abs(U[k:n,k]))
i = i+k
U[[k,i],k:n] = U[[i,k],k:n] # swap row i and k
L[[k,i],0:k] = L[[i,k],0:k] # swap row i and k
P[[k,i]] = P[[i,k]]         # swap row i and k
for j in range(k+1,n):
L[j,k] = U[j,k]/U[k,k]
U[j,k:n] = U[j,k:n] - L[j,k]*U[k,k:n]
return L,U,P


This performs solution of $$LUx = Pb$$ in two steps. In first step, solve $$Ly = Pb$$ In second step, solve $$Ux = y$$

In :
def LUSolve(L,U,P,b):
n = L.shape
# solve Ly = Pb
pb = b[P]
y = 0.0*pb
for i in range(n):
y[i] = (pb[i] - L[i,0:i].dot(y[0:i]))/L[i,i]
#solve Ux = y
x = 0.0*pb
for i in range(n-1,-1,-1):
x[i] = (y[i] - U[i,i+1:n].dot(x[i+1:n]))/U[i,i]
return x


Now we test the above function for LU decomposition.

In :
n = 3
A = np.random.rand(n,n)
L,U,P = LU(A)
print A
print L
print U
print P
print P.dot(A) - L.dot(U)

[[ 0.41300456  0.43021602  0.38608235]
[ 0.41796129  0.86040814  0.0505991 ]
[ 0.09585373  0.27291233  0.11822924]]
[[ 1.          0.          0.        ]
[ 0.98814069  1.          0.        ]
[ 0.22933638 -0.17997989  1.        ]]
[[ 0.41796129  0.86040814  0.0505991 ]
[ 0.         -0.41998827  0.33608332]
[ 0.          0.          0.16711326]]
[1 0 2]
[[ 0.18675073  0.11563255  0.57194172]
[ 0.19170746  0.54582466  0.23645847]
[ 0.50885829  0.70312835  0.50431159]]


Solve the linear system.

In :
b = np.random.rand(n)
x = LUSolve(L,U,P,b)
print A.dot(x) - b

[  0.00000000e+00  -1.11022302e-16  -1.11022302e-16]