Maple - LU Factorization with Partial Pivoting

[Rupert11]

1. I am asked to write a procedure that will inverse square matices using LU factorization with partial pivoting.



2. I am also told that the procedure should return the inverse matrix and report an error if it cannot do so.



3. So far I've come up with the code below but am not confident it is correct and I am also really struggling to find help with the next step for solving the linear systems i need too proceed with. Any adivce or suggestions:
restart :

n:= n; # all matrices are nxn matrices
A:= array (1..n, 1..n);
L:= array (1..n, 1..n);
U:= array (1..n, 1..n);
LU:= array (1..n, 1..n);

print('A is', A);

# carry out the Gaussian elimination
for j from 1 to n-1 do
for i from j+1 to n do
# calculate the multiplier for the ith row at the jth stage
m:= A[i,j]/A[j,j];
A[i,j]:= m;
for k from j+1 to n do
A[i,k]:= A[i,k]-m*A[j,k];
od:
od:
od:
print ('new A is', A);

# Extract L and U from the new A
for i from 1 to n do
for k from 1 to n do
if i<k then
U[i,k]:= A[i,k]
L[i,k]:= 0;
elif i=k then
U[i,k]:= A[i,j]
L[i,k]:=1;
else
U[i,k]:= 0
L[i,k]:= A[i,j]
fi;
od:
od:

print('L is',L);
print('U is',U);

# verify that A=LU
for i from 1 to n do
for k from 1 to n do
LU[i,k]:= add(L[i,j]*U[i,k], k=1..n);
od:
od:

print('LU is', LU);

 

[Rupert11]

So, I've been working on this problem for a few hours today and have come up with this procedure (below) but need advice on how to change it so that it works for a square matrix only and then finally how to form the inverse matrix to finish the question. Any ideas welcome...
# LU decomposition using Gaussian elimination with partial pivoting
# lufactor does the LU factorization
# lusolve solves a linear system using the LU factorization
lufactor := proc(A)
local n, r, k, p, i, temp, m, j;

n := rowdim(A); # get number of rows in matrix A
r := vector(n); # create pivot vector
for k from 1 to n do r[k] := k od; # initialize pivot vector

for k from 1 to n-1 do # loop over columns

# We find the maximal element in absolute value among the
# entries A[r[k],k], A[r[k+1],k], ... A[r[n],k].
# Interchange entries in pivot vector

p := k;
for i from k+1 to n do
if (abs(A[r,k]) > abs(A[r[p],k])) then p := i; fi;
od;
if A[r[p],k]=0 then
ERROR(`Matrix is singular`); # procedure will exit here
fi;
if (p <> k) then # interchange entries k and p in pivot vector
temp := r[k];
r[k] := r[p];
r[p] := temp;
fi;

# now do a gaussian elimination step

for i from k+1 to n do # loop down column below diagonal
m := A[r,k]/A[r[k],k]; # get multiplier
for j from k+1 to n do # subtract multiple of row i from
A[r,j] := A[r,j] - m*A[r[k],j]; # row k
od;
A[r,k] := m; # save multiplier
od;
od;
RETURN(op(r)); # return pivot vector
end:

# solve the linear system Ax=b given the LU factorization in a returned
# from ludecomp, the pivot vector r returned from ludecomp and the right
# hand side vector b.

lusolve := proc(A,r,b)
local n, y, x, i, s, j;
n := rowdim(A); # get number of rows in matrix A
y := vector(n); # create vector for solution of Ly=Pb
x := vector(n); # create vector for solution of Ux=y

# do forward substitution to solve Ly=Pb

y[1] := b[r[1]];
for i from 2 to n do
s := b[r];
for j from 1 to i-1 do
s := s - A[r,j]*y[j];
od;
y := s;
od;

# do backward substitution to solve Ux=y

x[n] := y[n]/A[r[n],n];
for i from n-1 to 1 by -1 do
s := y;
for j from i+1 to n do
s := s - A[r,j]*x[j];
od;
x := s / A[r,i];
od;
RETURN(op(x)); # return solution vector
end: