% lufac.m
% This file demonstrates LU-factorization, equivalently
% Gaussian elimination.
clear

N=4;
A=randn(4,4);
b=randn(4,1);

% solve Ax=b
x=A\b;

% now do it by explicit row operations
A1=A;
b1=b;
pivfac1=A1(2,1)/A1(1,1);	% Pivot is A(1,1)
A1(2,:)=A1(2,:)-pivfac1*A1(1,:);
b1(2)=b1(2)-pivfac1*b1(1);
% norm(x-A1\b1,'fro')	% check that we have the same solution
pivfac2=A1(3,1)/A1(1,1);
A1(3,:)=A1(3,:)-pivfac2*A1(1,:);
b1(3)=b1(3)-pivfac2*b1(1);
% norm(x-A1\b1,'fro')	% check that we have the same solution
pivfac3=A1(4,1)/A1(1,1);
A1(4,:)=A1(4,:)-pivfac3*A1(1,:);
b1(4)=b1(4)-pivfac3*b1(1);
% norm(x-A1\b1,'fro')	% check that we have the same solution
pivfac4=A1(3,2)/A1(2,2);	% Pivot is A(2,2)
A1(3,:)=A1(3,:)-pivfac4*A1(2,:);
b1(3)=b1(3)-pivfac4*b1(2);
% norm(x-A1\b1,'fro')	% check that we have the same solution
pivfac5=A1(4,2)/A1(2,2);
A1(4,:)=A1(4,:)-pivfac5*A1(2,:);
b1(4)=b1(4)-pivfac5*b1(2);
% norm(x-A1\b1,'fro')	% check that we have the same solution
pivfac6=A1(4,3)/A1(3,3);	% Pivot is A(3,3)
A1(4,:)=A1(4,:)-pivfac6*A1(3,:);
b1(4)=b1(4)-pivfac6*b1(3);
norm(x-A1\b1,'fro')	% check that we have the same solution

% now A1 is U by construction, what is L? remember A=L*U
L=A*inv(A1);
% however, we can do it without inverting any matrix
L=eye(N);
L(2,1)=pivfac1;
L(3,1)=pivfac2;
L(4,1)=pivfac3;
L(3,2)=pivfac4;
L(4,2)=pivfac5;
L(4,3)=pivfac6;
L, A1, norm(A-L*A1,'fro')
pause

% Finally let's do the LDU
v=diag(A1);
D=diag(v);
U=diag(1./v)*A1;	% The inversion is now scalar.
L, D, U, norm(A-L*D*U,'fro')


