Using the Gauss algorithm for $LU$-factorization on banded matrices will result in banded $L$ and $U$ factors with the same bandwith. Thus, it suffices to work on the diagonals. If there are only 3 diagonals, then the algorithm reduces to the Thomas algorithm.
Here is a Mathematica implementation that uses Compile to produce runtime optimized libraries for machined real-valued tridiagonal matrices:
cThomasLUDecomposition = Compile[{{l0, _Real, 1}, {d0, _Real, 1}, {u0, _Real, 1}},
Block[{n, u, uk, d, dold, l, lk},
n = Length[d0];
d = Table[0., {n}];
u = Table[0., {n}];
l = Table[0., {n}];
d[[1]] = dold = Compile`GetElement[d0, 1];
Do[
u[[k]] = uk = Compile`GetElement[u0, k];
l[[k]] = lk = Compile`GetElement[l0, k] / dold;
d[[k + 1]] = dold = Compile`GetElement[d0, k + 1] - lk uk;
,
{k, 1, n - 1}];
{l, d, u}
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
cThomasLUSolve = Compile[{{l, _Real, 1}, {d, _Real, 1}, {u, _Real, 1}, {y, _Real, 1}},
Block[{n, x, xold},
n = Length[d];
x = Table[0., {n}];
(*Solving x = LinearSolve[L,y];*)
x[[1]] = xold = Compile`GetElement[y, 1];
Do[
x[[k]] =
xold = Compile`GetElement[y, k] - Compile`GetElement[l, k - 1] xold
, {k, 2, n}];
(*Solving x = LinearSolve[U,x];*)
x[[n]] = xold = Compile`GetElement[x, n]/Compile`GetElement[d, n];
Do[
x[[k]] =
xold = (Compile`GetElement[x, k] - Compile`GetElement[u, k] xold)/Compile`GetElement[d, k]
, {k, n - 1, 1, -1}];
x
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];
Here a usage example:
A = {{-10.4151, 9.4098, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {105.526, -104.76,
94.4741, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 106.107, -104.511, 93.893,
0, 0, 0, 0, 0, 0, 0}, {0, 0, 106.749, -104.396, 93.2507, 0, 0, 0,
0, 0, 0}, {0, 0, 0, 107.459, -104.413, 92.5409, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 108.244, -104.564, 91.7564, 0, 0, 0, 0}, {0, 0,
0, 0, 0, 109.111, -104.854, 90.8894, 0, 0, 0}, {0, 0, 0, 0, 0, 0,
110.069, -105.29, 89.9312, 0, 0}, {0, 0, 0, 0, 0, 0, 0,
111.128, -105.883, 88.8723, 0}, {0, 0, 0, 0, 0, 0, 0, 0,
112.298, -106.643, 87.702}, {0, 0, 0, 0, 0, 0, 0, 0, 0, -11.1525,
3.02486}};
a = Diagonal[A, -1];
b = Diagonal[A, 0];
c = Diagonal[A, +1];
n = Length[b];
y = RandomReal[{-1, 1}, n];
(*Compute the factorization only one.*)
{l, d, u} = cThomasLUDecomposition[a, b, c];
(*Now you can use that factorization to solve with as many right-hand sides as you like.*)
x = cThomasLUSolve[l, d, u, y];
Max[Abs[A . x - y]]
2.84217*10^-14
By the way, the factors $L$ and $U$ can be obtained as follows:
U = DiagonalMatrix[SparseArray[d]] + DiagonalMatrix[SparseArray[Most[u]], +1];
L = DiagonalMatrix[SparseArray[ConstantArray[1., Length[d]]]] + DiagonalMatrix[SparseArray[Most[l]], -1];
Max[Abs[L . U - A]]
2.84217*10^-14
