If Cut/Copy and Paste fails, then click here for download.
/*
*
* Solves the linear system `transpose(U)*U*x=L*transpose(L)*x=b' where `U'
* and `L' are nonsingular bidiagonal matrices. The diagonals of `U' and
* `L', respectively, are stored in the vector `d', and any of the secondary
* diagonals is stored in the vector `e'. On input, `b' is the RHS. On
* output, `b' contains the solution `x'.
*
*/
void
cholsoltridiag (unsigned int n, double *d, double *e, double *b)
{
unsigned int i;
if (n == 0) {
return;
}
*(b + 0) /= *(d + 0);
for (i = 1; i < n; i++) {
*(b + i) = (*(b + i) - *(e + i - 1) * *(b + i - 1)) / *(d + i);
}
*(b + n - 1) /= *(d + n - 1);
for (i = n - 1; i > 0; i--) {
*(b + i - 1) = (*(b + i - 1) - *(e + i - 1) * *(b + i)) / *(d + i - 1);
}
return;
}