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
end

% 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;
         end
      end
%     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);
      end
   end

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

%  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);
end % for k

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

% 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);
   end
   ek = 1.0;
   if (tmp<0)
      ek = -1.0;
   end
   y(k) = -(ek + tmp)/A(k,k);
end % for k

% 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;
end

% 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;
end