Could you please check my code and help me to make the right one? mine is not working correctly. Please do not put yours just modified mine. And I tried side codes it didn't work. The GE part can be wrong, I am not sure. Thanks!

The main code is:

function [L,U,p,cri]=lu_fac_pp(A)

cri=1;

n=size(A,1);

p=zeros(n,1);

L=zeros(n);

U=zeros(n);

for k=1:n

p(k)=k;

end

for k=1:n-1

[mav,r]=max(abs(A(k:n,k)));

r=r+k-1;

if mav <(n*10^(-15))

disp('A is not invertible')

cri=0;

return

end

L=triL(A,-1)+eye(n);

U=triU(A);

if(r ~= k)

b([r k])=b([k r]);

A([r k],:)=A([k r],:);

end

p=A(k+1:n,k)/A(k,k);

A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-p*A(k,k+1:n);

b(k+1:n)=b(k+1:n)-p*b(k);

A(k+1:n,k)=0;

end

if abs(A(n,n)) < n*10^(-15)

disp('A is not invertible')

cri=0;

return

end

L=triL(A,-1)+eye(n);

U=triU(A);

cri=0;

return

end

Command window :

>> A=[1,1,0,3;2,1,-1,1;3,-1,-1,2;-1,2,6,-1]

>> b=[4;1;-3;4]

>>[ L,U,p,cri]=lu_fac_pp(A)

L=[1,0,0,0;-1/3,1,0,0;2/3,1,1,0;1/3,4/5,3/10,1]

U=[3,-1,-1,2;0,5/6,17/3,-1/3;0,0,-6,0;0,0,0,13/5]

p=[3;4;2;1]

>>c=perm_b(b,p)

c=[-3;4;1;4]

>>d=ult_sys(L,c)

d=[-3;3;0;2.6]

>> x=ut_sys(U,d)

x=[-1;2;0;1]

Update: I am trying to write a code for LU factorization with partial pivot of A matrix.

L= lower triangle matrix

U= upper

p=permutation

cri= criteria for checking A is invertible or not

mav(abs)=the max possible pivot entry

Inside suppose to include Gauss elimination and I need to write side codes to find d, c, and x.