function main_LU_pivot

DRIVER_partial_pivot
DRIVER_complete_pivot

return

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

function DRIVER_partial_pivot

A = [ 1, 0, 1, 2;...
      2, 2, 3, 4;...
     -1, 4, 2, 1;...
      1, 2, 1, 1,];
A = rand(7);

%%% Initialize L, U, and I.
m = size(A,2);
U = A;
L = eye(m);
I = 1:m; % Note that I is a VECTOR. It encodes the permutation P.

%%% Perform LU factorization.
for k = 1:(m-1)
    
  %%% Find the pivot index.
  [~,ipiv] = max(abs(U(k:m,k)));
  
  %%% Swap the rows, if necessary
  if ~(ipiv == 1)
    ipiv                = ipiv + k - 1;
    U([ipiv,k],:)       = U([k,ipiv],:);
    L([ipiv,k],1:(k-1)) = L([k,ipiv],1:(k-1));
    I([ipiv,k])         = I([k,ipiv]);
  end
  
  %%% Compute the k'th column of L and update U.
  J      = (k+1):m;
  lvec   = (1/U(k,k))*U(J,k);
  L(J,k) = lvec;
  U(J,k) = zeros(m-k,1);
  U(J,J) = U(J,J) - lvec*U(k,J); % This is where "all" flops are spent!!
  
  %%% Print results
  fprintf(1,'U = \n')
  disp(U)
  fprintf(1,'L = \n')
  disp(L)
  fprintf(1,'|A(I,:) - L*U| = %10.3e\n',max(max(abs(A(I,:) - L*U))))
  keyboard
  
end

fprintf(1,'Done!!\n')
keyboard

return

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

function DRIVER_complete_pivot

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

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

%%% Perform LU factorization.
for k = 1:(m-1)
    
  %%% Find the pivot index.
  [ipiv,jpiv] = find_max(U(k:end,k:end));
  ipiv        = ipiv + k - 1;
  jpiv        = jpiv + k - 1;
  
  %%% Swap the rows, if necessary
  if ~(ipiv == k)
    U([ipiv,k],:) = U([k,ipiv],:);
    L([ipiv,k],1:(k-1)) = L([k,ipiv],1:(k-1));
    I([ipiv,k]) = I([k,ipiv]);
  end
  
  %%% Swap the columns, if necessary
  if ~(jpiv == k)
    U(:,[jpiv,k]) = U(:,[k,jpiv]);
    J([jpiv,k]) = J([k,jpiv]);
  end
  
  %%% Compute the k'th column of L and update U.
  ind        = (k+1):m;
  lvec       = (1/U(k,k))*U(ind,k);
  L(ind,k)   = lvec;
  U(ind,k)   = zeros(m-k,1);
  U(ind,ind) = U(ind,ind) - lvec*U(k,ind); % This is where "all" flops are spent!!
  
  %%% Print results
  fprintf(1,'U = \n')
  disp(U)
  fprintf(1,'L = \n')
  disp(L)
  fprintf(1,'|A(I,J) - L*U| = %10.3e\n',max(max(abs(A(I,J) - L*U))))
  keyboard
  
end

fprintf(1,'Done!!\n')
keyboard

return

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

function [i,j] = find_max(A)

aa    = max(abs(A));
[~,j] = max(aa);
[~,i] = max(abs(A(:,j)));

return
