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?

  • 1
    $\begingroup$ 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. $\endgroup$ – Kirill Sep 24 '16 at 21:02
  • $\begingroup$ 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. $\endgroup$ – Lukas Bystricky Sep 24 '16 at 21:50
  • 1
    $\begingroup$ 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. $\endgroup$ – Bill Greene Sep 24 '16 at 22:55
  • $\begingroup$ 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.) $\endgroup$ – Abhilash Reddy M Sep 26 '16 at 19:34

The problem is that 2d and 3d discretizations are not block tridiagonal. The only tridiagonal decomposition would be a 2x2 decomposition that encompasses the entire matrix.

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.

  • $\begingroup$ I misunderstood something in the reference. Thank you for clearing that up. $\endgroup$ – Lukas Bystricky Oct 3 '16 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.