Can the Thomas algorithm be used to solve banded matrices of any arbitrary size? (ie larger than tridiagonal, penta/septa/+diagonal systems)

Question

I am working a product concerning the Poisson Equation. In this project, I end up with a very large square matrix. It's impossible to solve efficiently with simple Gaussian elimination, but the matrix has interesting features.

First, the matrix is very sparse: the vast majority of elements are a zero. On top of this, the matrix is symmetric and positive-defined. Finally, and most importantly concerning this question, it is banded. However, the length of the bandwidth changes based on problem to problem.

I was looking at the Thomas Algorithm, used for solving tridiagonal banded matrices, and I was wondering if the Thomas Algorithm could be implemented in a way that could solve a banded matrix of any arbitrary bandwidth? If so, could someone please explain how I would go about doing that or maybe link a research paper or article that explains it? And if not, does anyone have any other ways that a matrix like this could be solved efficiently? I currently have Cholesky Decomposition implemented to solve this, but I can be (and need to be) more efficient. Thank you all very much for the assistance.

Answer 1

If you can use an iterative method, conjugate gradient (for symmetric positive definite matrices ) and BICGSTAB (for general matrices) are fairly easy to implement and fast. You can find source code for these methods at math.nist.gov. LAPACK includes linear system solvers for banded matrices. Finally, if you want to implement a linear system solver for banded matrices yourself, take a look at this document.