function b = tspdSolve(d, s, b)
# usage: b = tspdSolve(d, s, b)
#
# Solve tridiagonal symmetric and positive definite
# system using the output from tspdChol.m
# Do not call when tspdChol has detected a singularity
# or non positive definiteness
#
# Inputs
#	d	main diagonal of the Cholesky factor L
#	s	subdiagonal of L 
#	b	right-hand side vector, overwritten by solution
# Output
#	b	solution

n = length(b);
#
# forward solve L*y = b
b(1) = b(1)/d(1);
for i=2:n
    b(i) = (b(i) - s(i-1)*b(i-1))/d(i);
endfor
# backward solve L^T*b = y
b(n) = b(n)/d(n);
for i=n-1:-1:1
    b(i) = (b(i) - s(i)*b(i+1))/d(i);
endfor;

endfunction