# solving tridiagonal system with multiple right hand sides

I need to solve a tridiagonal system (positive definite, diagonally dominant) $Ax = b$ in a time stepping loop. $A \in \mathbb{R}^{N \times N}$ remains constant but $b$ changes during each time step.

I was thinking which of the following would be faster :

1. Compute and store $\mathbf{z_i}$ that solves $A\mathbf{z_i} = \mathbf{e_i} \; \forall 1\leq i\leq N$. Then any $b \in \mathbb{R}^{N}$ is a linear combination of $\mathbf{z_i}$ and $x = b_i \mathbf{z_1} + \dotsb b_N \mathbf{z_N}$

2. Naively solve $A x = b$ for each given $b$ by tridiagonal matrix algorithm.

I thought (1) will save the time taken for computing the tridiagonal solves since the $z$ essentially gives the inverse of $A$. However, I calculated that one solve for $\mathbf{e_i}$ takes $8N$ flops, so that is $8N \times N$ to compute all the $\mathbf{z_i}$. Then for each given $b$, it takes $N^2$ flops to compute $x$. If I take $n$ time-steps, that is a total of $8N^2 + nN^2$.

In contrast (2) needs $n\times 8N$.

Am I correct ? Where is (2) gaining the advantage ?

• Since $A$ remains constant, the usual preference is to store $A$ in its factored form $LL^T$, since the back-solve and forward-solve are $O(N)$ complexity.
– hardmath
Apr 2, 2016 at 15:42
• That $N^2$ may be missing a coefficient 2 in front of it, but your analysis looks correct to me. As usual, computing inverses explicitly is the wrong choice -- not at all surprising: the same happens for dense matrices. Apr 2, 2016 at 15:47

## 1 Answer

The most efficient way would likely be to perform one step of Gauss-Elimination on A to eliminate the sub-diagonal and then store the inverse of the new main diagonal, the new modified super-diagonal and the constants to preform forward substitution on b.

This way each rhs will require one forward substitution to modify b followed by one back-solve. The forward operations consist of N additions and N multiplications and the backward operations consist of 2N multiplications and N additions. This adds up to 5N operations and no divisions (which are slower).