function [L,U] =  luquick(A)
%  LUQUICK Step-by-step LU factorization.
%          [L U] = LUQUICK(A) computes and displays the step-by-step
%          the LU factorization of its argument A, 
%          which must be a square, non-singular matrix. 

%  Karen E. Donnelly, Saint Joseph's College, 7/94
%  (c) Saunders College Publishing 1995
%  For use with MATLAB Manual: Computer Laboratory Exercises, Section 7.1
%  with  MATLAB Version 4.0.

%clc; disp(' '); disp(' ')
disp(' Demonstration of LU decomposition for square nonsingular matrix A')
disp(' For each step, shows')
disp('     Row operation required')
disp('     Adjusted (possibly permuted) lower triangular L so far')
disp('     Upper triangular matrix U so far')
disp(' Upon completion of the decomposition the following summary is shown:')
disp('     The final L and U for A = L*U')

[m n] = size(A);
if m ~=n 
  error (' Input argument must be square matrix')
end
tol = m*eps*norm(A,'inf');
if abs(det(A)) < tol
  error(' Matrix is singular')
end

format compact;
U = A;
L = eye(n);
col = 1;
numops = 0;
invElem = [];
while (col < n)
    
    % adjust if zero on the diagonal
   
    nzcol = col;
    while  (nzcol <= n) && (abs(U(nzcol,col)) < tol)
      nzcol = nzcol + 1;
    end

    if nzcol > col
      disp('     Press any key to continue'); pause; clc
      disp(' Upper triangularization of A before next transformation is ')
      U
      numops = numops + 1;
      str = [' Step ',num2str(numops),':'];
      disp(str)
      disp(' Zero along the diagonal -- swap with row below which is not')
      disp(' Note: This means L will be a permutation of lower triangular')
      U([col nzcol],1:n) = U([nzcol col],1:n)
      disp(' Swap corresponding columns of L')
      L(1:n,[col nzcol]) = L(1:n,[nzcol col])
      
    else  % not zero on the diagonal

      % get zeros in remaining positions below
      row = col+1;
      while ( row <= n)
         
        if abs(U(row, col)) > tol
           disp('      Press any key to continue ...'); pause; clc   
           disp(' Upper triangularization of A before next transformation is')
           U
           numops = numops + 1;
           str = [' Step ',num2str(numops),':'];
           disp(str)
           c = U(row,col)/U(col,col);
           if c < 0 
              str = [' Add ',num2str(-c), ' times row ',num2str(col),' to row ',num2str(row)];
           else
              str = [' Subtract ', num2str(c), ' times row ', num2str(col), ' from row ', num2str(row)];  
           end
           disp(str)
           U(row,:) = U(row,:) -c*U(col,:)
           str = [' Put ', num2str(c),' in row ',num2str(row),', column ',num2str(col),' of L to get'];
           disp(str)
           L(row,col) = c
                           
        end 

        row = row+1;
      end % while 

      % zero out negligible values.
      for row = 1:n
        if abs(U(row,col))< tol
           U(row,col) = 0;
        end
      end % for

      col = col + 1;
    end % if zero on diagonal -- else ---
    
end; % while

disp('      Press any key to continue ...'); pause; clc
disp('   Decomposition is complete -- Summary:' )
disp('   The matrix ');
A
disp('   is decomposed as L*U where')
L
U
disp('       Press any key to end'); pause; clc
format loose