4
$\begingroup$

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?

Improve this question
$\endgroup$
4
  • 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
    Commented Sep 24, 2016 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
    Commented Sep 24, 2016 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
    Commented Sep 24, 2016 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
    Commented Sep 26, 2016 at 19:34

1 Answer 1

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I misunderstood something in the reference. Thank you for clearing that up. $\endgroup$
    – Lukas Bystricky
    Commented Oct 3, 2016 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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