%
%	performs GEPP (Gaussian Elimination and row interchanges)
%
%       input: An n-by-m matrix A. 
%
%       output: Permutation PL (stored in a vector) and L and U matrices such that 
%
%                         A(PL,:) = L * U.
%
%       This code uses much more storage than necessary, and is not as efficient as 
%       it could be.
%
%       To solve a linear system of equations A x = b, where A is a square matrix, 
%       do the following in matlab:
%
%         [L,U,PL] = gepp(A); 
%          x  = U\(L\b(PL));
%
%       Ming Gu, Fall 2008
%
	function [L,U,PL] = gepp(A); 
%
%
%       first determine the dimensions of A and set up the necessary arrays.
%
	[n,m] = size(A);
	U = A; 
	if (m > n)
	   L = eye(n,n);
	else
	   L = zeros(n,m);
	   L(1:m,1:m) = eye(m);
	end
	PL = (1:n)'; 
        Amax = max(max(abs(A)));
%
%       now do n-1 steps of elimination.
%
	for k = 1:min(m,n) - 1
%
%          perform partial pivoting and permute all the matrices accordingly.
%
	   [y,i] = max(abs(U(k:n,k)));
	   kl = i + k - 1;
	   temp = U(kl,:);
	   U(kl,:) = U(k,:);
	   U(k,:) = temp;
	   temp = PL(kl);
	   PL(kl) = PL(k);
	   PL(k) = temp;
%
%          perform the elimination
%
           if (abs(U(k,k)) < (eps * Amax)) 
              disp('Input matrix is numerically singular; abort computation')
              return
           end
	   U(k+1:n,k) = U(k+1:n,k) /U(k,k);
           U(k+1:n,k+1:m) = U(k+1:n,k+1:m) - U(k+1:n,k) * U(k,k+1:m);
        end
%
%       now perform the n-th step of elimination if A has more rows than columns.
%
   	if (n > m)
	   [y,i] = max(abs(U(m:n,m)));
	   kl = i + m - 1;
	   temp = U(kl,:);
	   U(kl,:) = U(m,:);
	   U(m,:) = temp;
	   temp = PL(kl);
	   PL(kl) = PL(m);
	   PL(m) = temp;
           if (abs(U(m,m)) < (eps * Amax)) 
              disp('Input matrix is numerically singular; abort computation')
              return
           end
           U(m+1:n,m) = U(m+1:n,m) / U(m,m);
        end
%
%       setup the final L and U matrices.
%
        if (n > m) 
           L = L + tril(U,-1);
	   U = triu(U); 
           U = U(1:m,1:m);
        else
           U1 = tril(U,-1);
           L  = L + U1(1:n,1:n);
           U  = triu(U);
           clear U1;
        end
        return