# Comparing algorithms for tridiagonal linear systems solution

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}$$

• I looked around a bit at some of the primary resources online, and I couldn't find any papers in the obvious places (LAPACK Working Notes, etc...) discussing this particular routine, and the LAPACK routine for this, xgtsl, still bears Jack's original copyright. Neither of your approaches employs pivoting, which is probably a more important factor than the other differences between them. – Aron Ahmadia Apr 19 '12 at 12:05
• @AronAhmadia: That's true, unless the tridiagonal system is diagonally dominant, in which case no partial pivoting is necessary. – Paul Apr 20 '12 at 4:44

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).
• If you want to ask a new question, that's fine, though you should just look up the definition of pivoting, there are plenty of easily-constructed cases where you create a divide-by-zero or divide-by-near-zero problem. Regarding the article describing the implementation, there may not be much, since it is a fairly trivial extension of LU. – Aron Ahmadia Apr 20 '12 at 14:00