function [p,L,U]=goodfact(A)
%
%  function [p,L,U]=goodfact(A)
%    
% This is how a standard code for Gaussian elimination with
% partial pivoting might be written in the storage for A
%  Usually the permutation matrix P is stored in an index vector p.
[m,n]=size(A);
U=zeros(n,n);
L=eye(m,n);
if m== n
   s=n-1;
else
   s=min([m n]);
end;
p=1:m;
for k=1:s
 [z,kmax]=max(abs(A(k:m,k)));
 kmax=kmax+k-1;
 if kmax > k    
  ptemp=p(k);
  p(k) = p(kmax);
  p(kmax)=ptemp;
  x=A(k,:);
  A(k,:)=A(kmax,:);
  A(kmax,:)=x;
 end;
% Compute new rows of U and L
 U(k,k:n)=A(k,k:n);
 L(k,1:k-1)=A(k,1:k-1);
% This line computes the multipliers
% Here it is one vector statement
 A(k+1:m,k)=A(k+1:m,k)/U(k,k);
% This line modifies that matrix. It is an outer product of
% a column vector A(k+1:n,k) and a row vector A(k,k+1:n)
% that modifies the active (n-k)x(n-k) part of A
 A(k+1:m,k+1:n)=A(k+1:m,k+1:n) -A(k+1:m,k)*U(k,k+1:n);
end;
if m == n
U(n,n)=A(n,n);
L(n,1:n-1)=A(n,1:n-1);
end;
%
% If A is rectangular, we need these lines.
if m > n 
    L(n+1:m,1:n)=A(n+1:m,1:n);
end;
    
% This lines give us the MATLAB L and U with the permutation
% included in L
%
% To get the L that the MATLAB routine lu outputs add this line.
% L(p,:)=L;