function [A, rcondA, pvt] = decomp(A)
# usage: [A, rcondA, pvt] = decomp(A)
#--------------------------------------------------------------------
# triangularize a real matrix using Gaussian elimination            
# with or without partial pivoting and estimate its condition number.
#
# input:
#   A: matrix to be triangularized, to be overwritten by L and U
#
# outputs:
#   A: contains an upper triangular matrix U and a permuted
#      version of a lower triangular matrix L so that
#      (permutation matrix) P*A = L*U, where P is determined by pvt,
#      L = eye(n) + tril(A, -1), and U = triu(A)
#   rcondA: an estimate of the reciprocal of the condition of A for
#      the linear system
#      A*x = b. If rcondA+1.0 = 1.0, A is singular to working
#      precision. Set to 0 if exact singularity is detected.
#   pvt: the pivot vector
#      pvt(k) is the index of the kth pivot row
#      pvt(n) = (-1)**(number of interchanges)
#      optional argument, without this argument, no pivoting is
#      performed
#
# The determinant of A can be obtained on output by
# pvt(n)*A(1,1)*A(2,2)*...*A(n,n)
#
# Dependency
#   solve.m
#-------------------------------------------------------------------

n = length(A(:,1));

normA = norm(A,1);
if (nargout>2)
#  initial pivot vector
   pvt = zeros(n,1);
   pvt(n) = 1;     # pvt(n) = (-1)**0
endif

# Gaussian elimination
for k=1:n-1
   if (nargout>2)
#     find pivot row
      m = abs(A(k, k));
      pvt(k) = k;
      for i=k+1:n,
         if (abs(A(i,k)) > m)
            m = abs(A(i,k));
            pvt(k) = i;
         endif
      endfor
#     alternatively, use max function
#      [tmp, ind] = max(abs(A(k:n, k)));
#      pvt(k) = ind + k -1;
#     interchange row k and the pivot row pvt(k)
      if (k~=pvt(k))
         A([k, pvt(k)],:) = A([pvt(k),k],:);
         pvt(n) = -pvt(n);
      endif
   endif

   if (A(k,k)==0.0)
      rcondA = 0.0;
      return
   endif

#  compute multipliers
   A(k+1:n,k) = A(k+1:n,k)/A(k,k);
#  update by block, BLAS2 version 
   A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
endfor #k

if (A(n,n)==0.0)
   rcondA = 0.0;
   return
endif

# estimate cond(A)
# this involves solving two systems A^T*y = e
# and A*z = y, where the components of e are 1 or -1
# so that y is heuristically large.
# Applying permutations is unnecessary for cond estimation.
#
# solve A^T*y = e
#
# solve U^T*work = e
y = zeros(n,1);
for k=1:n
   if (k==1)
      tmp = 0.0;
   else
      tmp = A(1:k-1,k)'*y(1:k-1);
   endif
   ek = 1.0;
   if (tmp<0)
      ek = -1.0;
   endif
   y(k) = -(ek + tmp)/A(k,k);
endfor

# solve L^T*y = work
for k=n-1:-1:1
#  update the right hand side
   tmp = A(k+1:n,k)'*y(k+1:n);
   y(k) = y(k) - tmp;
endfor

# solve A*z = y
z = solve(A,y);

rcondA = norm(y,1)/(normA*norm(z,1));
# rcondA is at most 1
if (rcondA>1.0)
   rcondA = 1.0;
endif

endfunction