Algorithm Implementation/Linear Algebra/Tridiagonal matrix algorithm

All the provided implementations of the tridiagonal matrix algorithm assume that the three diagonals, a (below), b (main), and c (above), are passed as arguments.

$$\begin{bmatrix} b_0 & c_0 & 0 & 0 & 0 & \cdots & 0 & 0 \\ a_1 & b_1 & c_1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & a_2 & b_2 & c_2 & 0 & \cdots & 0 & 0 \\ 0 & 0 & a_3 & b_3 & c_3 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & a_{n-3} & b_{n-3} & c_{n-3} & 0 \\ 0 & 0 & 0 & \cdots & 0 & a_{n-2} & b_{n-2} & c_{n-2} \\ 0 & 0 & 0 & \cdots & 0 & 0 & a_{n-1} & b_{n-1} \\ \end{bmatrix} \cdot \begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \\ \vdots \\ x_{n-3} \\ x_{n-2} \\ x_{n-1} \\ \end{bmatrix} = \begin{bmatrix} d_0 \\ d_1 \\ d_2 \\ d_3 \\ \vdots \\ d_{n-3} \\ d_{n-2} \\ d_{n-1} \\ \end{bmatrix}$$

Haskell
note that the function takes in a subdiagonal list (as) whose first element (index 0) is on row 1, which is why we shift the index in the convenience accessors so we can call it "a(1)" anyway.

C
The following is an example of the implementation of this algorithm in the C programming language.

The following variant preserves the system of equations for reuse on other inputs. Note the necessity of library calls to allocate and free scratch space - a more efficient implementation for solving the same tridiagonal system on many inputs would rely on the calling function to provide a pointer to the scratch space.

C++
The following C++ function will solve a general tridiagonal system (though it will destroy the input vector c and d in the process). Note that the index $$i$$ here is zero based, in other words $$i = 0, 1, \dots, N - 1$$ where $$N$$ is the number of unknowns. Note that, save for the printing of text in the main function, this code is valid C as well as C++.

Python
Note that the index $$i$$ here is zero-based, in other words $$i = 0, 1, \dots, N-1$$ where $$N$$ is the number of unknowns.

Another way

Fortran 90
Note that the index $$i$$ here is one based, in other words $$i = 1, 2, \dots, n$$ where $$n$$ is the number of unknowns.

Sometimes it is undesirable to have the solver routine overwrite the tridiagonal coefficients (e.g. for solving multiple systems of equations where only the right side of the system changes), so this implementation gives an example of a relatively inexpensive method of preserving the coefficients.

This subroutine offers an option of overwriting d or not.