function [U,p] = PGaussElim(A) %Gauss Elimination for Nonsingular Matrix
        % This is certainly NOT the right way to do this with Matlab, which
        % has sophisticated matrix manipulation routines. This is provided
        % merely to illustrate the algorithm in Olver
        %
        %  To use, create a matrix 
        %   > A = rand(5,10)
        % and type
        %  > [U,p] = GaussElim(A)
        % p is a vector of pointers that gives the permutation (see Olver)
        % U is the upper triangular matrix from Row Reducing A
[m,n] = size(A);
U = A;
p = 1:m;                 % Permutation Pointers
for j = 1:n              % Choose jth column
    i = j;
    while (U( i,j) ==0)
        i = i+1;
      if(i > m)
         error('Matrix is singular')
      end
    end
    ptemp = p(j);
    p(j) = p(i);
    p(i) = ptemp;
    for i = j+1:m         % Pivot on element U(p(j),j)
        l = U(p(i),j)/U(p(j),j);  
      for k = j:n         % Row p(i) -> Row p(i) - l*Row p(j)
            U(p(i),k) = U(p(i),k) -l*U(p(j),k);
      end
    end
end