function [d,s,rcond]=tspdChol(d,s)
# usage: [d,s,rcond]=tspdChol(d,s)
#
# The Cholesky factorization of a tridiagonal symmetric and
# positive definite matrix
#
# Inputs
#	d	main diagonal, overwritten by output
#	s	subdiagonal, overwritten by output
# Outputs
#	d	main diagonal of the Cholesky factor
#	s	subdiagonal of the Cholesky factor 
#       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.
# Dependency
#   tspdSolve.m

n = length(d);
#
if (d(1) == 0.0)
    rcond = 0.0;
    return
endif

if (d(1) < 0)
    disp("The matrix is not positive definite.");
    return
endif

d(1) = sqrt(d(1));

for i=1:n-1,
    s(i) = s(i)/d(i);
    d(i+1) = d(i+1) - s(i)*s(i);
    if (d(i+1) == 0.0)
        rcond = 0.0;
        return
    endif
    if (d(i+1) < 0)
        disp("The matrix is not positive definite.");
        return
    endif
    d(i+1) = sqrt(d(i+1));
endfor

# condition number estimation

# 1-norm
norm1 = norm([d(1);s(1)],1);            # column 1
if (n>2)
    tmp = max(abs(s(1:n-2)) + abs(d(2:n-1)) + abs(s(2:n-1)));
                                        # columns 2 to n-1
    if tmp>norm1
       norm1 = tmp;
    endif
endif
tmp = norm([s(n-1);d(n)],1);            # column n
if norm1<tmp
   norm1 = tmp;
endif

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

# solve A*z = y
z = y;
z = tspdSolve(d,s,z);

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

endfunction