%	perform GECP
%
	function [L,U,PL,PU] = gecp2(A);
%
%       input: An n-by-m matrix A. 
%
%       output: Permutations PL and PU (stored in vectors) and L and U matrices such that  
%
%                         A(PL,PU) = L * U.
%
%       To solve a linear system of equations A x = b, where A is a square matrix, 
%       do the following in matlab:
%
%         x      = zeros(size(b));
%         [L,U,PL,PU] = gecp2(A); 
%         x(PU)  = U\(L\b(PL));
%
%       Note: this code is not as efficient as it could be.
%
%       Ming Gu, Fall 2008
%
%       first determine the dimensions of A and set up the necessary arrays.
%
	[n,m] = size(A);
	if (m > n)
	   L = eye(n,n);
	else
	   L = zeros(n,m);
	   L(1:m,1:m) = eye(m);
	end
	PL = (1:n)'; PU = (1:m)';
        Amax = max(max(abs(A)));
        if (Amax == 0) 
           disp('Input matrix is numerically singular; abort computation')
           return
        end
        tol = eps;
%
%       now do min(m,n)-1 steps of elimination.
%
	for k = 1:min(m,n) - 1
%
%          perform complete pivoting and permute all the matrices accordingly.
%
	   [Y,I] = max(abs(A(k:n,k:m)));
	   [y,i] = max(Y);
	   ku = i + k - 1;
	   kl = I(i) + k - 1;
	   temp = A(kl,:);
	   A(kl,:) = A(k,:);
	   A(k,:) = temp;
	   temp = A(:,ku);
	   A(:,ku) = A(:,k);
	   A(:,k) = temp;
	   temp = PL(kl);
	   PL(kl) = PL(k);
	   PL(k) = temp;
	   temp = PU(ku);
	   PU(ku) = PU(k);
	   PU(k) = temp;
%
%          perform the elimination
%
           if (abs(A(k,k)) < (tol * Amax)) 
              disp('Input matrix is numerically singular')
              return
           end
	   A(k+1:n,k) = A(k+1:n,k) /A(k,k);
           A(k+1:n,k+1:m) = A(k+1:n,k+1:m) - A(k+1:n,k) * A(k,k+1:m);
        end
%
%       now perform the m-th step of elimination if A has more rows than columns.
%
	if (n > m)
	   [y,i] = max(abs(A(m:n,m)));
	   kl = i + m - 1;
	   temp = A(kl,:);
	   A(kl,:) = A(m,:);
	   A(m,:) = temp;
	   temp = PL(kl);
	   PL(kl) = PL(m);
	   PL(m) = temp;
           A(m+1:n,m) = A(m+1:n,m) / A(m,m);
        end
%
%       now perform the n-th step of elimination if A has more columns than rows.
%
        if (m > n) 
	   [y,i] = max(abs(A(n,n:m)));
	   ku = i + n - 1;
	   temp = A(:,ku);
	   A(:,ku) = A(:,n);
	   A(:,n) = temp;
	   temp = PU(ku);
	   PU(ku) = PU(n);
	   PU(n) = temp;
        end
%
%       setup the final L and U matrices.
%
        U = A;
        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