# Thomas algorithm for 3D finite difference

For 1D finite difference, the resulting linear system is tri-diagonal and can be solved in O(n) using the Thomas algorithm.

I am trying to solve a finite difference system in 3D. Numerical Methods by Quarteroni et al goes over some theory behind extending the Thomas algorithm to block linear systems. This method seems to be ideal for me because it is still O(n).

From an algorithmic point of view however, I'm not sure how to proceed. Is there a reference that covers the implementation of this algorithm for a linear system such as the one resulting from 3D finite difference?

• The point of Thomas algorithm, which is really just a special-cased GE/LU for a tridiagonal matrix, is that a tridiagonal matrix has tridiagonal LU decomposition. But a 3d finite difference matrix wouldn't have a tridiagonal LU decomposition, it would be $n^2$-banded. Sep 24 '16 at 21:02
• Right, it would be a banded matrix so of course the algorithm wouldn't be identical. The reference I posted (section 3.8) goes over some variations for banded systems. I'm just trying to work out the implementation. Sep 24 '16 at 21:50
• What is your objective? If it is to compute a solution to your equations with minimal CPU time, I suggest you simply call the Lapack routine dgbsv in a high-performance implementation of Lapack/BLAS (e.g. OpenBLAS or MKL) and be done with it. Sep 24 '16 at 22:55
• What PDE do you want to solve? For diffusion problems, you might be able to use alternating-direction-implicit (ADI) methods. It would allow to apply Thomas algorithm to 3D problems (at the cost of being restricted to 2nd order spatial accuracy.) Sep 26 '16 at 19:34

Of course, this fact is a great pity. If there was a tridiagonal decomposition, we would have a way of getting an $O(N)$ algorithm from a purely algebraic perspective (i.e., by just looking at the matrix). Instead, the only way to get $O(N)$ algorithms are ones that employ multigrid methods that take into account where the matrix actually comes from.