function [P,L,U] = plumovie(A,speed,tol);
%PLUMOVIE   This function is equivalent to plu(A), except it has 
%           been made into a "movie" to show the progress of the 
%           calculation. Starting with U = A and PinvL = eye(A) 
%           the same row operations are successively applied to 
%           both matrices  until U is an upper echelon matrix. 
%           PinvL is then L^(-1)*P, a unit lower triangular matrix 
%           up to a permutation of columns.
%
%           plumovie(A) produces the factorization PA = LU.
%           plumovie(A,speed) produces the factorization PA = LU. 
%           If the argument speed is equal to 'slow', then the 
%           program pauses after each operation until a key is 
%           pressed.
%           plumovie(A,speed,tol) uses the given tolerance in the 
%           rank tests.

[m,n] = size(A);
P = eye(m);
PinvL = eye(m);
U = A;
clc

if (nargin < 2), speed = 'fast'; end 

home, disp(['Setup: U = A, PinvL = eye' blanks(10)])
U
PinvL
if speed == 'slow', pause, end

% Compute the default tolerance if none was provided.
if (nargin < 3), tol = max(m,n)*eps*norm(A,'inf'); end


% Loop over the entire matrix.
i = 1;
j = 1;
while (i <= m) & (j <= n)
   % Find value and index of largest element in the remainder 
   % of column j.
     [p,k] = max(abs(U(i:m,j))); k = k+i-1;
   %    
   if (p <= tol)
      % The column is negligible, zero it out.
      home, disp(['column ' int2str(j) ' is negligible' blanks(10)])
      U(i:m,j) = zeros(m-i+1,1)
      j = j + 1;
      if speed == 'slow', pause, end
   else
      % put pivot element in i-th row   
      % Swap i-th and k-th rows.
      home 
      disp(['column ' int2str(j) ': swap rows '...
      int2str(i) ' and ' int2str(k) blanks(10)])
      U([i k],:) = U([k i],:)
      P([i k],:) = P([k i],:);
      PinvL([i k],:) = PinvL([k i],:)
      if speed == 'slow', pause, end
      % Subtract multiples of the pivot row from the lower rows.
      for k = [i+1:m]
         home, disp(['eliminate in column ' int2str(j) blanks(10)])
         multiplier = U(k,j)/U(i,j);
         U(k,j:n) = U(k,j:n) - multiplier*U(i,j:n)
         PinvL(k,:) = PinvL(k,:) - multiplier*PinvL(i,:)
      if speed == 'slow', pause, end
      end
      i = i + 1;
      j = j + 1;
      if speed == 'slow', pause, end
   end
end
L = P*inv(PinvL);
if speed == 'slow', pause, end
clc
P
L
U
% end of plumovie.  --SSp 5.2.89