Using MATLAB develop an Mfile to determine LU factorization

Using MATLAB. develop an M-file to determine LU factorization of a square matrix with partial pivoting. That is, develop a function called mylu that is passed the square matrix [A] and returns the triangular matrices [L] and [U] and the permutation P. You are not to use MATLAB built-in function lu in your codes. Test your function by using it to solve a system of equations listed below in part 3. Confirm that your function is working properly by verifying that [L][U]=P[AJ and by using the MATLAB built-in function lu.

Solution

function [L,U,P] = mylu(A)
% LU factorization.
%
%   
% [L,U,P] = Lu(A) returns unit lower triangular matrix L, upper
% triangular matrix U, and permutation matrix P so that P*A = L*U.
%
%
s=length(A);
U=A;
L=zeros(s,s);
PV=(0:s-1)\';
for j=1:s,
% Pivot Voting (Max value in this column first)
[~,ind]=max(abs(U(j:s,j)));
ind=ind+(j-1);
t=PV(j); PV(j)=PV(ind); PV(ind)=t;
t=L(j,1:j-1); L(j,1:j-1)=L(ind,1:j-1); L(ind,1:j-1)=t;
t=U(j,j:end); U(j,j:end)=U(ind,j:end); U(ind,j:end)=t;

% LU
L(j,j)=1;
for i=(1+j):size(U,1)
c= U(i,j)/U(j,j);
U(i,j:s)=U(i,j:s)-U(j,j:s)*c;
L(i,j)=c;
end
end
P=zeros(s,s);
P(PV(:)*s+(1:s)\')=1;

function Ainv=MatrixInverse(A)

% Do Lu decompositon to obtain triangle matrices (can easily be inverted)
[L,U,P] = mylu(A);

% Solve linear system for Identity matrix
I=eye(size(A));
s=size(A,1);
Ainv=zeros(size(A));
for i=1:s
b=I(:,i);
Ainv(:,i)=TriangleBackwardSub(U,TriangleForwardSub(L,P*b));
end

function C=TriangleBackwardSub(U,b)
% Triangle Matrix Backward Substitution
%
% Solves C = U \\ B;
%
% |1| |2 2 1|
% With b = |2| and U = |0 1 4|
% |3| |0 0 3|
%
s=length(b);
C=zeros(s,1);
C(s)=b(s)/U(s,s);
for j=(s-1):-1:1
C(j)=(b(j) -sum(U(j,j+1:end)\'.*C(j+1:end)))/U(j,j);
end

function C=TriangleForwardSub(L,b)
% Triangle Matrix Forward Substitution
%
% Solves C = L \\ b
%
% |1| |1 0 0|
% With b = |2| and L = |2 1 0|
% |3| |3 4 1|
%
s=length(b);
C=zeros(s,1);
C(1)=b(1)/L(1,1);
for j=2:s
C(j)=(b(j) -sum(L(j,1:j-1)\'.*C(1:j-1)))/L(j,j);
end