%LUEXAMPLE   script to demonstrate Gaussian elimination as a matrix
%decomposition  PA = LU 

format rat  % makes all entries appear as rational numbers;
            % note "*" denotes a nearly zero entry (i.e. up to 14 digits)

A=[11 -4 -2  7;
   -2  7  5  5; 
    4  3  1  7;
   -1  9 11 -4]

P1=eye(4)

P1*A

m21=-2/11;  m31=4/11;  m41=-1/11;

L1=[1    0 0 0;
    -m21 1 0 0;
    -m31 0 1 0;
    -m41 0 0 1]

L1*P1*A 

P2=[1 0 0 0; 0 0 0 1; 0 0 1 0; 0 1 0 0]

P2*L1*P1*A

u2=P2*L1*P1*A;  m32=u2(3,2)/u2(2,2);  m42=u2(4,2)/u2(2,2);

L2=[1 0    0 0;
    0 1    0 0;
    0 -m32 1 0;
    0 -m42 0 1]

L2*P2*L1*P1*A

P3=eye(4)

P3*L2*P2*L1*P1*A

u3=P3*L2*P2*L1*P1*A;  m43=u3(4,3)/u3(3,3);

L3=[1 0 0    0;
    0 1 0    0;
    0 0 1    0;
    0 0 -m43 1]

U=L3*P3*L2*P2*L1*P1*A

% so   A = (L3*P3*L2*P2*L1*P1)^{-1} *U

Ltilde=inv(L3*P3*L2*P2*L1*P1)

P=P3*P2*P1

L=P*Ltilde

P*A

L*U  %  it equals P*A !