function [A,flag,pivots,Cond] = Factor(A)
%  Factor decomposes the matrix A using Gaussian elimination and
%  estimates its condition number.  Factor is used in conjunction 
%  with Solve to solve A*x = b.  REQUIRES Solve.m.
%
%  Input variables:
%     A       = matrix A of neq rows to be decomposed
%  Output variables:
%     A       = the upper triangular matrix U in its upper portion 
%               (by rows) and a permuted version of a lower triangular
%               matrix I-L such that (permutation matrix)*A = L*U.
%               A record of interchanges is kept in pivots.
%    flag     = reports whether the matrix A has a zero pivot.  A value 
%               of zero means no zero pivot was encountered. If positive,
%               the first zero pivot occurred at the flag value's equation 
%               and the decomposition could not be completed. 
%    Cond     = an estimate of the condition number of A (unless flag
%               is nonzero).
%    pivots   = record of row interchanges.  Also,
%                  pivots(neq) = (-1)**(number of interchanges).
%
%  The determinant of A can be obtained on output from
%  det(A) = pivots(neq) * A(1,1) *  ... * A(neq,neq).

[neq,cols] = size(A);
flag = 0;
Cond = 0;
pivots = ones(neq,1);
temp = zeros(neq,1);

if neq == 1 
  if A(1,1) == 0
    flag = 1;
  else
    Cond = 1;
  end
  return
end

%  Compute infinity-norm of A for later condition number estimation.
Anorm = norm(A,inf);

%  Gaussian elimination with partial pivoting.
for k = 1:neq-1
  %  Determine the row m containing the largest element in
  %  magnitude to be used as a pivot.
  [biggest,occurred] = max(abs(A(k:neq,k)));
  m = k + occurred - 1;

  %  Check for a nonzero pivot.  If all possible pivots are 
  %  zero, A is numerically singular.
  if biggest == 0
    flag = k;
    return 
  end
  pivots(k) = m;
  if m ~= k
    pivots(neq) = -pivots(neq);
    %   Interchange the current row k with the pivot row m.
    A([m k],k:neq) = A([k m],k:neq);
  end

  %  Eliminate subdiagonal entries of column k.
  for i = k+1:neq
    t = A(i,k) / A(k,k);
    A(i,k) = -t;
    if t ~= 0
      A(i,k+1:neq) = A(i,k+1:neq) - t * A(k,k+1:neq);
    end
  end
end
if A(neq,neq) == 0
  flag = neq;
  return
end

%  Estimate the condition number of A.
temp(1) = -1 / A(1,1);
for k = 2:neq
  t = A(1:k-1,k)' * temp(1:k-1);
  if t >= 0
    sgnt =  1;
  else
    sgnt = -1;
  end
  temp(k) = -(sgnt + t) / A(k,k);
end
for k = neq-1:-1:1
  temp(k) = sum(temp(k) * A(k+1:neq,k));
  m = pivots(k);
  temp([m k]) = temp([k m]);
end
dnorm = norm(temp,inf);
temp = Solve(A,pivots,temp);
ynorm = norm(temp,inf);
Cond = Anorm * ynorm / dnorm;

