function [L,U,pv,qv] = lugui(A,strategy,fast)
%LUGUI Gaussian elimination demonstration.
%
% from http://www.stanford.edu/class/cs138/mlmath.html
%
%   LUGUI(A) shows the steps in LU decomposition by Gaussian elimination.
%   At each step of the elimination, the pivot that would be chosen by
%   MATLAB's partial pivoting algorithm is shown in magenta.  You can use
%   the mouse to pick any pivot.  The pivot is shown in red, the emerging
%   columns of L in green, and the emerging rows of U in blue.
%
%   Popup menus allow the pivot strategy and speed to be changed dynamically.
%   lugui(A,'mouse'), the default, allows you to use the mouse to pick pivots.
%   lugui(A,'diagonal') uses diagonal elements as pivots.
%   lugui(A,'partial') uses the largest element in the unreduced column.
%   lugui(A,'complete') uses the largest element in the unreduced matrix.
%   lugui(A,'...','slow'), the default, pauses between steps.
%   lugui(A,'...','fast') executes faster.
%
%   [L,U,p,q] = lugui(A,...) returns a lower triangular L, an upper
%   triangular U and permutation vectors p and q so that L*U = A(p,q).

% Initialize

if nargin < 1
   n = 2 + ceil(6*rand);
   A = ceil(10*rand(n,n));
end
if nargin < 2
   strategy = 'mouse';
end
if nargin < 3
  fast = 'slow';
end

[m,n] = size(A);
shg
clf
dx = 100;
dy = 30;
set(gcf,'double','on','name','LU Gui', ...
   'menu','none','numbertitle','off','color','white', ...
   'pos',[480-(dx/2)*min(9,n), 320, (n+1)*dx, (m+3)*dy], ...
   'windowbuttonupfcn','set(gcf,''tag'',''pivot'')')
axes('pos',[0 0 1 1])
axis off
Lcolor = [0 .65 0];
Ucolor = [0  0 .90];
Acolor = [0 0 0];
Grey = [.5 .5 .5];
PartialPivotColor = [1 0 1];
PivotColor = [1 0 0];
TempColor = [1 1 1];

% Each element has its own handle

for j = 1:n
   for i = 1:m
      t(i,j) = text('units','pixels','string',spf(A(i,j)), ...
         'fontname','courier','fontweight','bold','fontsize',14, ...
         'horiz','right','color',Acolor, ...
         'pos',[20+j*dx 20+(m+2-i)*dy],'userdata',[i j], ...
         'buttondownfcn','set(gcf,''userdata'',get(gco,''userdata''))');
   end
end

% Menus

switch lower(strategy(1))
   case 'd', val = 2;
   case 'p', val = 3;
   case 'c', val = 4;
   otherwise, val = 1;
end
pivotstrat = uicontrol('pos',[20+(dx/2)*(n-2) 20 160 20],'style','pop', ...
   'val',val,'fontsize',12,'back','white','string',{'Pick a pivot', ...
   'Diagonal pivoting','Partial pivoting','Complete pivoting'});

switch lower(fast(1))
   case 'f', val = 2;
   otherwise, val = 1;
end
speed = uicontrol('pos',[200+(dx/2)*(n-2) 20 100 20],'style','pop', ...
   'val',val,'fontsize',12,'back','white','string', ...
   {'Slow','Fast','Pause','Close'});

% Elimination

pv = 1:m;
qv = 1:n;
for k = 1:min(m,n)
   if all(all(~finite(A(k:m,k:n))))
      break
   end
   if (m == n) & (k == n)
      p = n;
      q = n;
   else
      pp = min(find(abs(A(k:m,k)) == max(abs(A(k:m,k)))))+k-1;
      set(t(pp,k),'color',PartialPivotColor)
      p = 0;
      q = 0;
      while p < k | q < k
         switch get(pivotstrat,'val')
   
            case 1  % Pick a pivot with mouse
               pq = get(gcf,'userdata');
               if isequal(get(gcf,'tag'),'pivot') & ~isempty(pq)
                  p = pq(1);
                  q = pq(2);
                  set(gcf,'tag','','userdata',[])
               else
                  drawnow
               end
   
            case 2 % Diagonal pivoting
               p = k;
               q = k;
   
            case 3 % Partial pivoting
               p = pp;
               q = k;
   
            case 4 % Complete pivoting
               [p,q] = find(abs(A(k:m,k:n)) == max(max(abs(A(k:m,k:n)))));
               p = p(1)+k-1;
               q = q(1)+k-1;
         end

         switch get(speed,'val')
            case 1
               paws = 0.05;
            case 2
               paws = 0;
            case 3
               p = 0;
               q = 0;
               drawnow
            case 4
               close
               return
         end
      end
   end
   set(t(pp,k),'color',Acolor)
   set(t(p,q),'color',PivotColor)

   if get(speed,'val')==2, paws = 0; else, paws = 0.05; end 
   pause(10*paws)

%  Swap columns

   A(:,[q,k]) = A(:,[k,q]);
   qv([q,k]) = qv([k,q]);
   for s = .05:.05:1
      for i = 1:m
         set(t(i,k),'pos',[20+(k+s*(q-k))*dx 20+(m+2-i)*dy])
         set(t(i,q),'pos',[20+(q+s*(k-q))*dx 20+(m+2-i)*dy])
      end
      drawnow
   end
   t(:,[q,k]) = t(:,[k,q]);
   for i = 1:m
      set(t(i,k),'string',spf(A(i,k)),'userdata',[i k])
      set(t(i,q),'string',spf(A(i,q)),'userdata',[i q])
   end
   pause(10*paws)

%  Swap rows

   A([p,k],:) = A([k,p],:);
   pv([p,k]) = pv([k,p]);
   for s = .05:.05:1
      for j = 1:n
         set(t(k,j),'pos',[20+j*dx 20+(m+2-(k+s*(p-k)))*dy])
         set(t(p,j),'pos',[20+j*dx 20+(m+2-(p+s*(k-p)))*dy])
      end
      drawnow
   end
   t([p,k],:) = t([k,p],:);
   pause(10*paws)

   for j = k:n
      set(t(k,j),'string',spf(A(k,j)),'userdata',[k j])
      set(t(p,j),'string',spf(A(p,j)),'userdata',[p j])
   end
   pause(10*paws)

%  Skip step if column is all zero

   if all(A(k:m,k) == 0)
      for i = k+1:m
         set(t(i,k),'string',spf(A(i,k)),'color',Lcolor)
         drawnow
      end
      for j = k:n
         set(t(k,j),'string',spf(A(k,j)),'color',Ucolor)
         drawnow
      end
   else

%     Compute multipliers in L

      for i = k+1:m
         A(i,k) = A(i,k)/A(k,k);
         set(t(i,k),'string',spf(A(i,k)),'color',Lcolor)
         pause(paws)
         drawnow
      end

%     Elimination    

      for j = k+1:n
         for i = k+1:m
            set(t(i,j),'color',TempColor)
            drawnow
            pause(paws)
            A(i,j) = A(i,j) - A(i,k)*A(k,j);
            set(t(i,j),'string',spf(A(i,j)),'color',Acolor)
            drawnow
            pause(paws)
         end
         if get(speed,'val')==2, paws = 0; else, paws = 0.05; end 
      end

      for j = k:n
         set(t(k,j),'string',spf(A(k,j)),'color',Ucolor)
         drawnow
      end
      pause(paws)
   end
   if k < min(m,n), pause(10*paws), end
end

% Seperate L and U into two matrices

delete(pivotstrat)
delete(speed)

for s = .1:.1:1.5
   for j = 1:n
      for i = 1:m
         if i <= j
            set(t(i,j),'pos',[20+(j+.10*s)*dx 20+(m+2-i)*dy])
         else
            set(t(i,j),'pos',[20+(j-.10*s)*dx 20+(m+2-s-i)*dy])
         end
      end
   end
   drawnow
end

% Insert ones on diagonal of L
 
r = min(m,n);
for j = 1:r
   text('units','pixels','string',spf(1.0), ...
      'fontname','courier','fontweight','bold','fontsize',14, ...
      'horiz','right','color',Lcolor, ...
      'pos',[20+(j-0.15)*dx 20+(m+.5-j)*dy]);
end
drawnow

if nargout > 0
   L = tril(A(:,1:r),-1) + eye(m,r);
   U = triu(A(1:r,:));
end

%------------------------------------------------------------

function s = spf(aij);
% Subfunction to format text strings
if aij == 0
   f = '%10.0f';
elseif (abs(aij) < 1.e-4) | (abs(aij) >= 1.e4) 
   f = '%10.1e';
else
   f = '%10.4f';
end
s = sprintf(f,aij);