function [u,d,l,rcond]=decompt(u,d,l)
% usage: [u,d,l,rcond]=decompt(u,d,l)
%
% The LU decomposition of a tridiagonal matrix
% using Gaussian elimination without pivoting
%
% Inputs
%	u	the upper diagonal, overwritten by output
%	d	the main diagonal, overwritten by output
%	l	the subdiagonal, overwritten by output
% Outputs
%	u	the upper diagonal of U
%	d	the main diagonal of U
%	l	the subdiagonal of L 
%       rcond   estimate of 1/cond 
%               If 1.0+rcond=1.0, the matrix is singular to
%               working precision. Set to 0 if exact singularity
%               is detected.

n = length(d);
%
if (d(1) == 0.0)
    rcond = 0.0;
    return
end

for i=1:n-1,
    l(i) = l(i)/d(i);			% multiplier
    d(i+1) = d(i+1) - l(i)*u(i);	% update
    if (d(i+1) == 0.0)
        rcond = 0.0;
        return
    end
end

% condition number estimation

% 1-norm
norm1 = norm([d(1);l(1)],1);            % column 1
tmp = max(abs(l(2:n-1)) + abs(d(2:n-1)) + abs(u(1:n-2)));
                                        % columns 2 to n-1
if tmp>norm1
   norm1 = tmp;
end
tmp = norm([u(n-1);d(n)],1);            % column n
if norm1<tmp
   norm1 = tmp;
end

% solving A^T*y=e, where the components of e are 1 or -1
% so that y is heuristically large
%
% solve U^T*work = e
y = zeros(n,1);
y(1) = -1.0/d(1);
for k=2:n
   tmp = u(k-1)*y(k-1);
   ek = 1.0;
   if (tmp<0)
      ek = -1.0;
   end
   y(k) = -(ek + tmp)/d(k);
end % for k
% solve L^T*y = work
for k=n-1:-1:1
   y(k) = y(k) - l(k)*y(k+1);
end;

% solve A*z = y
z = y;
z = solvet(u,d,l,z);

rcond = norm(y,1)/(norm1*norm(z,1));
% rcond is at most 1
if (rcond>1.0)
   rcond = 1.0;
end