After a bit more thought, I think it's actually fairly simple to modify TDMA to deal with arbitrary diagonal offsets:
Forward sweep:
$$ \begin{align} c'i & = \begin{cases} \cfrac{c_i}{b_i} & ; & i=1,2,\dots,k \ \cfrac{c_i}{b_i-c'{i-k} a_i} & ; & i=k+1,k+2,\dots,n-k \end{cases} \ \ d'i &= \begin{cases} \cfrac{d_i}{b_i} & ; & i=1,2,\dots,k \ \cfrac{d_i-d'{i-k}a_i}{b_i-c'_{i-k}a_i} & ; & i=k+1,k+2,\dots,n. \ \end{cases} \end{align} $$
Back substitution:
$$ \begin{align} x_n &= d'n\ x_i &= d'_i-c'_ix{i+k} & ; {}& i=n-1,n-2,\dots,1. \end{align} $$
I've written a quick implementation in Cython:
cimport numpy as cnp
# NB: `b` and `d` are modified in place
cpdef TDMA_offset(const double[:] a, double[:] b, const double[:] c, double[:] d):
cdef size_t i, k, m, n
m = a.shape[0]
n = b.shape[0]
k = n - m
# Forward sweep
for i in range(m):
d[i+k] -= d[i]*a[i]/b[i]
b[i+k] -= c[i]*a[i]/b[i]
# Back substitution
for i in range(m-1, -1, -1):
d[i] -= d[i+k]*c[i]/b[i+k]
for i in range(n):
d[i] /= b[i]
This approach still doesn't make use of the fact that $M$ is symmetric, and a parallel implementation would also be very nice.
Following the comment by <PERSON>, I made a simple Python test script:
import numpy as np
from cymodules import _solvers
from scipy import sparse
def test_trisolve(n=10, k=5):
a = np.random.randn(n - k)
b = np.random.randn(n)
c = np.random.randn(n - k)
x = np.random.randn(n)
M = sparse.diags((a, b, c), (-k, 0, k), format='csc')
d = M.dot(x)
# solve for x in place
x_hat = d.copy()
_solvers.TDMA_offset(a, b, c, x_hat)
print "||x - x_hat|| = ", np.linalg.norm(x - x_hat)
I've tested my solver for a variety of values of $n$ and $k$, including cases where $k \geq (n - 1) / 2$:
test_trisolve(n=10, k=3)
# ||x - x_hat|| = 1.67369803253e-15
test_trisolve(n=10, k=6)
# ||x - x_hat|| = 6.85231143806e-16
It seems to be working OK. 262