% Gaussian elimination with complete pivoting
% function [Pr,Pc,L,U] = gecp(A)
% Input:
%    A = matrix to factor as Pr'*A*Pc = L*U, or A = Pr*L*U*Pc'
%        A may have more rows than columns
% Output:
%    Pr = permutation matrix, i.e. matrix whose rows are the
%         rows of the identity in a permuted order.
%    Pc = permuations matrix. 
%         This means Pr'*A*PC is the same as A with its rows
%         and columns permuted.
%    L  = lower triangular matrix with ones on its diagonal;
%         each subdiagonal |L(i,j)| <= 1
%    U  = upper triangular matrix;
%         each superdiagonal |U(i,j)| <= |U(i,i)|
%
function [Pr,Pc,L,U] = gecp(A)
[n,m]=size(A);
% assume n >= m
Pr=eye(n);                             % Initialize Pr to the identity matrix
Pc=eye(m);                             % Initialize Pc to the identity matrix
for i=1:(min(m,n-1)),                  % For each column
   Am=max(max(abs(A(i:n,i:m))));          % Find largest value among |A(i:n,i:m)|
   [I,J]=find(abs(A(i:n,i:m)) == Am);     % Find location(s) of largest value
   imax=I(1)+i-1;                         % Shift (first) location to be in rows i:n
   jmax=J(1)+i-1;                         % Shift (first) location to be in cols i:n
   if (imax ~= i),                        % Swap rows of A and Pr if necessary
      temp = Pr(:,i);
      Pr(:,i) = Pr(:,imax);
      Pr(:,imax) = temp;
      temp = A(i,:);
      A(i,:) = A(imax,:);
      A(imax,:) = temp;
   end
   if (jmax ~= i),                        % Swap columns of A and Pc if necessary
      temp = Pc(:,i);
      Pc(:,i) = Pc(:,jmax);
      Pc(:,jmax) = temp;
      temp = A(:,i);
      A(:,i) = A(:,jmax);
      A(:,jmax) = temp;
   end
%  At this point |A(i,i)| is largest entry among |A(i:n,i:m)|
%
   A(i+1:n,i) = A(i+1:n,i)/A(i,i);        % Compute column of L
   if i+1 <= m,                           % Update trailing submatrix
      A(i+1:n,i+1:m) = A(i+1:n,i+1:m) - A(i+1:n,i)*A(i,i+1:m);
   end 
end
L = tril(A,-1) + eye(n,m);             % Extract L from A
U = triu(A);                           % Extract U from A