function [LU, p] = gauss( A )
  %--------------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % February 28, 1999
  %
  % Usage: [LU, p] = gauss( A )
  %
  % Returns:
  %	LU:	The permuted LU factorization of A
  %	p:	Permutation vector
  %
  % Parameters:
  %     A:      Matrix to be LU factored via Gaussian elimination.
  %
  %
  % THIS IS INEFFICIENT CODE - IT WAS WRITTEN FOR DEMONSTRATION PURPOSES.
  %
  % This function performs Gaussian elimination with scaled partial pivoting 
  % on a square matrix a.  It returns the LU factorization of a and
  % a permutation vector that indicates the true position of the rows and 
  % columns in the result.
  %
  % The companion function solve() can be used to solve Ax=b once A has been
  % factored by this function.
  %--------------------------------------------------------------------------

  [n, m] = size( A );
  if ( n ~= m )
    fprintf( 'ERROR: input matrix is not square\n' );
    return;
  end

  % Generate initial index vector, filled with number 1 to n.

  p = [ 1 : n ];

  % Determine the largest (in magnitude) element in each row.  These factors
  % are used to scale the pivot elements for comparison purposes when deciding
  % which row to use as a pivot row.

  s  = zeros( n , 1 );                     % Note: As an alternative to this
  for i = 1 : n                            % double loop, the vector s can be
    smax = 0.0;                            % constructed with the line
    for j = 1 : n                          %   s = max(abs(A'))';
      smax = max( smax, abs( A(i,j) ) );   % The "max" operator finds the
    end                                    % maximum in each column of a
    s(i) = smax;                           % matrix, so we must examine a'
  end                                      % to get the maximum in each row.

  % Begin Gaussian elimination.

  for k = 1 : n - 1
    % Find the remaining row with the largest scaled pivot.
    rmax = 0.0;
    for i = k : n
      r = abs( A(p(i),k) / s(p(i)) );
      if ( r > rmax )
	rmax = r;
	j = i;
      end
    end

    % Row j has the largest scaled pivot, so "swap" that row with the
    % current row (row k).  The swap is not actually done by copying rows,
    % but by swapping two entries in an index vector.

    temp = p(j);
    p(j) = p(k);
    p(k) = temp;

    % Now carry out the next elimination step as usual, except for the
    % added complication of the index vector.

    for i = k + 1 : n
      xmult = A(p(i),k) / A(p(k),k);
      A(p(i),k) = xmult;
      for j = k + 1 : n
	A(p(i),j) = A(p(i),j) - xmult * A(p(k),j);
      end
    end
  end

  LU = A;

end
