function main_LU

%DRIVER_forwards
DRIVER_general
DRIVER_efficient
DRIVER_how_to_build_funny_example

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_forwards

A = [ 1, 0, 1, 2;...
      2, 2, 3, 4;...
     -1, 4, 2, 1;...
      1, 2, 1, 1,];

U = A;
fprintf(1,'A = U = \n')
disp(U)
keyboard

L1        = eye(4);
L1(2:4,1) = -(1/U(1,1))*U(2:4,1); % This is the vector "l1"/
U         = L1*U;
fprintf(1,'L1 = \n')
disp(L1)
fprintf(1,'U = L1*U = \n')
disp(U)
keyboard

L2        = eye(4);
L2(3:4,2) = -(1/U(2,2))*U(3:4,2);
U         = L2*U;
fprintf(1,'L2 = \n')
disp(L2)
fprintf(1,'U = L2*U = \n')
disp(U)
keyboard

L3        = eye(4);
L3(4,3)   = -(1/U(3,3))*U(4,3);
U         = L3*U;
fprintf(1,'L3 = \n')
disp(L3)
fprintf(1,'U = L3*U = \n')
disp(U)
keyboard

L = inv(L1)*inv(L2)*inv(L3);
fprintf(1,'L = inv(L1)*inv(L2)*inv(L3) = \n')
disp(L)
fprintf(1,'max|A - L*U| = %10.3e\n',max(max(abs(A - L*U))))
keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_general

A = [ 1, 0, 1, 2;...
      2, 2, 3, 4;...
     -1, 4, 2, 1;...
      1, 2, 1, 1,];
%A = randn(7);
fprintf(1,'A = \n')
disp(A)
keyboard
  
m = size(A,2);
U = A;
L = eye(m);
for i = 1:(m-1)
  LL            = eye(m);
  LL((i+1):m,i) = -(1/U(i,i))*U((i+1):m,i);
  U             = LL*U;
  L             = L*inv(LL);
  fprintf(1,'U = \n')
  disp(U)
  fprintf(1,'L = \n')
  disp(L)
  fprintf(1,'|A - L*U| = %10.3e\n',max(max(abs(A - L*U))))
  keyboard
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_efficient

A = [ 1, 0, 1, 2;...
      2, 2, 3, 4;...
     -1, 4, 2, 1;...
      1, 2, 1, 1,];
%A = randn(1000);

%%% Initialize L and U.
m = size(A,2);
U = A;
L = eye(m);

tic
%%% Perform LU factorization.
for i = 1:(m-1)
  J      = (i+1):m;
  lvec   = (1/U(i,i))*U(J,i);
  L(J,i) = lvec;
  U(J,i) = zeros(m-i,1);
  U(J,J) = U(J,J) - lvec*U(i,J);
  %fprintf(1,'U = \n')
  %disp(U)
  %fprintf(1,'L = \n')
  %disp(L)
  %fprintf(1,'|A - L*U| = %10.3e\n',max(max(abs(A - L*U))))
  %keyboard
end
fprintf(1,'Time = %0.5f\n',toc)

keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_how_to_build_funny_example

U  = [1,0,1,2;...
      0,2,1,0;...
      0,0,1,3;...
      0,0,0,2];
  
L3        = eye(4);
L3(4,3)   = 1;
L2        = eye(4);
L2(3:4,2) = [-2,-1];
L1        = eye(4);
L1(2:4,1) = [-2,1,-1];

A         = inv(L1)*inv(L2)*inv(L3)*U;

disp(inv(L3)*U)
disp(inv(L2)*inv(L3)*U)
disp(inv(L1)*inv(L2)*inv(L3)*U)
keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

