Symmetric, Tridiagonal, Positive Definite System Solver
(from Golub-Van Loan)
Given an n-by-n symmetric, positive definite matrix A and a vector b in R^n,
the following algorithm overwrites b with the solution to Ax=b. It is assumed
that the diagonal of A is stored in d(1:n) and the superdiagonal in e(1:n-1).
for k=2:n
t=e(k-1);
e(k-1)=t/d(k-1);
d(k)=d(k)-te(k-1);
end
for k=2:n
b(k)=b(k)-e(k-1)b(k-1)
end
b(n)=b(n)/d(n)
for k=n-1:-1:1
b(k)=b(k)/d(k)-e(k)b(k+1)
end