Comparing algorithms for tridiagonal linear systems solution

Question

Below there are two algorithms for solving tridiagonal linear systems of the form $$ \left[ \begin{array}{ccccc|c} b_1 & c_1 & & & &d_1\\ a_2 & b_2 & c_2 & & & d_2\\ & \ddots & \ddots & \ddots & & \vdots\\ & & a_{n-1} & b_{n-1} & c_{n-1} & d_{n-1}\\ & & & a_n & b_n & d_n \end{array} \right]. $$ I called them Algorithms A and B. Both of them are equivalent to Gaussian elimination, but with important difference in the form of the resulting triangular (bidiagonal) matrix.

My main question is: which one of them is more preferrable?

Algorithm A is the one that described in Wikipedia and many textbooks, it is called Thomas algorithm and is implemented, for example, in Numerical Recipes in some tricky form. Algorithm B is more straightforward and, in my opinion, is more numerically stable in cases when $|b_i|\gg|a_i|+|c_i|$ . Though I haven't seen Algorithm B in texbooks, note that exactly this algorithm is implemented in the mentioned Wikipedia article, see "Implementation in Fortran 90", while "Implementation in Matlab" deals with Algorithm A ("Implementation in C" in its current state is a mess that does not seem to work at all).

---

$$ \begin{array}{|c|c|}\hline \mathbf{Algorithm\ A} & \mathbf{Algorithm\ B}\\\hline \textit{% Elimination}&\textit{% Elimination}\\ \begin{array}{l} \tilde c_1=c_1/b_1\\ \tilde d_1=d_1/b_1\\ \mathbf{for }\quad i=2 \quad \mathbf{to}\quad n-1 \quad \textbf{do}\\ \quad q=b_i-a_i c_{i-1}\\ \quad \tilde c_i=c_i/q\\ \quad \tilde d_i=(d_i-a_i \tilde d_{i-1})/q\\ \mathbf{end do}\\ \tilde d_n=(d_n-a_n \tilde d_{n-1})/(b_{n}-a_n \tilde c_{n-1})\\ \\ \end{array} & \begin{array}{l} \\ \\ \hat b_1=b_1\\ \mathbf{for }\quad i=2 \quad \mathbf{to}\quad n \quad \textbf{do}\\ \quad q=a_i/\hat b_{i-1}\\ \quad \hat b_i=b_i-q c_{i-1}\\ \quad \hat d_i=d_i-q \hat d_{i-1}\\ \mathbf{end do}\\ \\ \\ \end{array}\\ \hline \textit{% Resulting system} & \textit{% Resulting system}\\ \left[ \begin{array}{ccccc|c} 1 & \tilde c_1 & & & &\tilde d_1\\ & 1 & \tilde c_2 & & & \tilde d_2\\ & & \ddots & \ddots & & \vdots\\ & & & 1 & \tilde c_{n-1} & \tilde d_{n-1}\\ & & & & 1 & \tilde d_n \end{array} \right] & \left[ \begin{array}{ccccc|c} \hat b_1 & c_1 & & & &\hat d_1\\ & \hat b_2 & c_2 & & & \hat d_2\\ & & \ddots & \ddots & & \vdots\\ & & & \hat b_{n-1} & c_{n-1} & \hat d_{n-1}\\ & & & & \hat b_n & \hat d_n \end{array} \right]\\ \hline \textit{% Backsubtitution} & \textit{% Backsubtitution}\\ \begin{array}{l} \\ x_n=\tilde d_n\\ \mathbf{for }\quad i=n-1 \quad \mathbf{downto}\quad 1 \quad \textbf{do}\\ \quad x_i=\tilde d_i-\tilde c_i x_{i+1}\\ \\ \end{array} & \begin{array}{l} \\ x_n=\hat d_n/\hat b_n\\ \mathbf{for }\quad i=n-1 \quad \mathbf{downto}\quad 1 \quad \textbf{do}\\ \quad x_i=(\hat d_i-c_i x_{i+1})/\hat b_i\\ \\ \end{array}\\\hline \end{array} $$

Answer 1

Both algorithms compute $LU$ decompositions (solving against $L$ while it is being formed) and then solve against the resulting $U$. The difference is that Algorithm A forces $U$ to have a diagonal of all ones (we say that $U$ is unit-diagonal), while Algorithm B forces $L$ to have a unit diagonal (this is the usual convention).

Regardless of whether or not one is more stable than the other, both are a bad idea; as @AronAhmadia mentioned, you should use an algorithm which performs partial pivoting. I would go with the LAPACK routine dgtsv.