# "Numerical Recipes" ludcmp() C code on page 46 translated into Python
# I make it as simple as possible, disregard efficiency.
# Here a is a list of list, n is integer (size of the matrix)
# index is a list, and d is also a list of size 1
# Python list index starts from 0. So matrix index is from 0 to n-1.
#
import math
def ludcmp(a, n, indx, d):
d[0] = 1.0
# looking for the largest a in each row and store it in vv as inverse
# We need a new list same size as indx, for this we use .copy()
vv=indx.copy()
for i in range(0,n):
big = 0.0
for j in range(0,n):
temp = math.fabs(a[i][j])
if (temp > big):
big = temp
vv[i]=1.0/big
#
# run Crout's algorithm
for j in range(0,n):
# top half & bottom part are combined
# but the upper limit l for k sum is different
big = 0.0
for i in range(0,n):
if(i=j):
dum = vv[i]*math.fabs(sum)
if (dum >= big):
big = dum
imax = i
# end if (i>= ...)
# end for i
# pivoting part, swap row j with row imax, a[j] is a whole row
if (j != imax):
dum = a[imax]
a[imax] = a[j]
a[j] = dum
d[0] = - d[0]
vv[imax] = vv[j]
# end if (j != ...)
# divide by the beta diagonal value
indx[j] = imax
dum = 1.0/a[j][j]
for i in range(j+1,n):
a[i][j] *= dum
# end for i
# end for j
# end of def ludcmp
# We do backward substitution in lubksb() take the row swapped LU decomposed
# a, size n, and swapping indx, and b vector as input. The output is
# in b after calling.
def lubksb(a, n, indx, b):
ii = -1
# forward
for i in range(0,n):
ip = indx[i]
sum = b[ip]
b[ip] = b[i]
if(ii != -1):
for j in range(ii,i):
sum -= a[i][j]*b[j]
elif(sum != 0):
ii = i
b[i] = sum
# bote alpha_{ii} is 1 above
# backward
for i in range(n-1,-1,-1):
sum = b[i]
for j in range(i+1,n):
sum -= a[i][j]*b[j]
b[i] = sum/a[i][i]
# end lubksb()
# unfortunately a is destroyed (become swapped LU)
def linearsolver(a,n,b):
indx = list(range(n))
d =[1]
ludcmp(a,n,indx,d)
x = b.copy()
lubksb(a,n,indx,x)
print("x=",x)
return x
# end linearsolver
a=[[1,3,3,-5],[2,-4,7,-1],[7,0.5,3,-6],[9,-2,3,8]]
n=4
b=[0,2,3,-10]
linearsolver(a,n,b)