% Solves Ax=b via Lz=b and Ux=z
% A=LU decomposition (with partial pivot)
% where A is triangular nxn matrix,
% L is a unit lower triangular (diagonal entries are 1's)
% and U is upper triangular

function [L,U,z,x]=triLUpvt(A,b)

n=length(A);

%set L=I and U=A

L=eye(n);  %identity matrix size n
U=A;       %u_ij = a_ij for any ij


for k=1:n-1
    
    %partial pivoting
    pvrow=U(k,k:n);
    for i=k:n
        umax=max(U(k:n),k);
        if (U(i,k) == umax)
           U(k,k:n)=U(i,k:n); %new pvrow
           U(i,k:n)=pvrow;    %replace by old pvrow
         end
    end
    
    %gaussian elimination
    if (U(k,k)==0)
    disp('no solution exists')
    end

    for j=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);

    end %for j loop
end %for k loop

%solves Lz=b for z

z(1)=b(1);
for j=2:n
    z(j)=b(j) - L(j,j-1)*z(j-1);
end % j-loop

%solves Ux=z for x

x(n)=z(n)/U(n,n);
for i=n-1:-1:1
    x(i)=(z(i) - U(i,i+1)*x(i+1) )/U(i,i);
end % j-loop

z=z';
x=x';






