O(n) TDMA in Python, C, and MATLAB with Stability Checks for Engineers
The Thomas algorithm, also called TDMA, solves tridiagonal linear systems in O(n) time and O(n) storage, far faster than the O(n^3) cost of dense Gaussian elimination. Use it when your coefficient matrix is strictly tridiagonal and either diagonally dominant or symmetric positive definite. When those conditions are not guaranteed, switch to a pivoted library routine such as LAPACK’s DGTSV instead of trusting plain elimination.
TL;DR:
The Thomas algorithm is most effective for strictly tridiagonal matrices that are diagonally dominant or symmetric positive definite, offering linear O(n) solve time.
It is not guaranteed stable for arbitrary matrices; pivoting routines like LAPACK’s DGTSV are recommended when diagonal dominance cannot be assured.
Handling periodic boundary conditions requires solving two standard systems and combining solutions with the Sherman-Morrison formula due to the matrix’s cyclic couplings.
Parallel algorithms such as cyclic reduction or hybrid GPU-optimized methods can outperform serial TDMA when solving many small systems or on hardware optimized for parallel processing.
Always verify sentinel values, perform pivot and residual checks, and consider stability conditions before deploying TDMA in critical applications.
Table of Contents
What is a tridiagonal system and where does it come from?
A tridiagonal system has nonzero entries only on the main diagonal and the two bands adjacent to it. Engineers typically store it as four vectors rather than a full matrix: a for the subdiagonal, b for the main diagonal, c for the superdiagonal, and d for the right-hand side. Row i of the system reads a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i, which is why the method is sometimes called a tridiagonal matrix algorithm rather than a general solver.
Coding conventions matter here because off-by-one errors are the most common bug in TDMA implementations. With zero-based indexing, the first subdiagonal entry a[0] and the last superdiagonal entry c[n-1] don’t correspond to real matrix entries, so they’re set to zero as sentinels. Some references use one-based indexing instead, which shifts every loop bound by one and trips up translations between Python, C, and MATLAB if you’re not careful.
Tridiagonal systems show up constantly in numerical computing, which is why the Thomas algorithm earns a permanent place in most engineers’ toolkits:
1D diffusion and Poisson equations: finite-difference discretization of a second derivative naturally produces three nonzero entries per row.
Implicit time-stepping schemes: methods like backward Euler or Crank-Nicolson for parabolic PDEs generate a new tridiagonal system at every time step.
Natural cubic spline interpolation: the continuity conditions on spline segments reduce to a tridiagonal system for the second derivatives at each knot.
Alternating direction implicit (ADI) methods: each sweep along one spatial direction solves a tridiagonal system per grid line, which is why TDMA is often called in a tight inner loop thousands of times per simulation.
Recognizing this pattern early saves you from reaching for a dense solver on a problem that doesn’t need one.
How do you derive the forward sweep and backward substitution?
TDMA is Gaussian elimination specialized to a matrix with only three nonzero bands per row, so fill-in outside those bands never occurs, which is exactly what keeps the arithmetic linear instead of cubic. The derivation proceeds in two passes over the data.
Forward sweep, modified coefficients. Starting from c’_1 = c_1 / b_1 and d’1 = d_1 / b_1, compute for i = 2 to n: m = b_i - a_i c’{i-1}; c’_i = c_i / m; d’i = (d_i - a_i d’{i-1}) / m.
Pivot check. Each step divides by m, the modified diagonal entry. If m is zero or extremely small relative to the other coefficients in that row, the algorithm is about to divide by an unreliable number and the result should not be trusted without further checks.
Backward substitution. Set x_n = d’_n, then for i = n-1 down to 1: x_i = d’_i - c’i x{i+1}.
Each step in the forward sweep depends only on the result of the previous step, which is both the source of the method’s efficiency and the reason it resists naive parallelization: you cannot compute row i’s modified coefficients before row i-1’s are known. That sequential dependency chain is the central trade-off that shows up again when we look at parallel alternatives later in this guide.
Pseudocode and implementations in Python, C, and MATLAB
The forward sweep and backward substitution translate almost line for line into any language. The pattern below assumes zero-based arrays with a[0] = 0 and c[n-1] = 0 as sentinels, matching the convention most production code uses.
Compact pseudocode:
for i = 1 to n-1:
m = a[i] / b[i-1]
b[i] = b[i] - m * c[i-1]
d[i] = d[i] - m * d[i-1]
x[n-1] = d[n-1] / b[n-1]
for i = n-2 downto 0:
x[i] = (d[i] - c[i] * x[i+1]) / b[i]
A Python implementation using NumPy, written to avoid mutating the caller’s arrays:
import numpy as np
def thomas(a, b, c, d):
n = len(b)
cp, dp = c.copy().astype(float), d.copy().astype(float)
bp = b.copy().astype(float)
x = np.zeros(n)
for i in range(1, n):
m = a[i] / bp[i-1]
bp[i] -= m * cp[i-1]
dp[i] -= m * dp[i-1]
x[-1] = dp[-1] / bp[-1]
for i in range(n-2, -1, -1):
x[i] = (dp[i] - cp[i] * x[i+1]) / bp[i]
return x
Casting to float before dividing avoids silent integer truncation if the caller passes integer arrays, a mistake that produces plausible-looking but wrong output.
A C implementation, written in place since embedded and HPC code rarely has memory to spare:
void thomas(int n, double *a, double *b, double *c, double *d, double *x) {
for (int i = 1; i < n; i++) {
double m = a[i] / b[i-1];
b[i] -= m * c[i-1];
d[i] -= m * d[i-1];
}
x[n-1] = d[n-1] / b[n-1];
for (int i = n-2; i >= 0; i--)
x[i] = (d[i] - c[i] * x[i+1]) / b[i];
}
Keeping the loop in a single forward pass, then a single backward pass, preserves cache-friendly sequential memory access, which matters more than any micro-optimization for a method that’s already O(n).
MATLAB/Octave users should remember that native array indexing starts at 1, so the loop bounds shift accordingly:
function x = thomas(a, b, c, d)
n = length(b);
for i = 2:n
m = a(i) / b(i-1);
b(i) = b(i) - m * c(i-1);
d(i) = d(i) - m * d(i-1);
end
x(n) = d(n) / b(n);
for i = n-1:-1:1
x(i) = (d(i) - c(i) * x(i+1)) / b(i);
end
end
Copy b and d before modifying them whenever the caller might reuse the originals, since the algorithm overwrites both during elimination.
Test every implementation against a known solution: build a small tridiagonal matrix, multiply by a chosen x, and confirm your solver recovers it within floating-point tolerance.
When is plain Thomas elimination unsafe to trust?
Plain Thomas elimination performs no pivoting, so its stability is not guaranteed for an arbitrary tridiagonal matrix. It behaves well when the matrix is diagonally dominant or symmetric positive definite, conditions common in discretized diffusion and Poisson problems, but outside those cases a tiny pivot can blow up rounding error across the whole solution.
Three checks catch most trouble before it reaches production:
Absolute pivot threshold: flag any modified diagonal entry whose magnitude falls below a small multiple of the matrix’s typical coefficient size.
Residual check: after solving, compute the residual Ax minus d and confirm it’s small relative to the input scale.
Diagonal dominance screen: run a quick pass over a, b, c before solving to confirm |b_i| meets or exceeds |a_i| + |c_i| for each row, which is a sufficient (not necessary) stability condition.
LAPACK’s DGTSV and DGTSVX routines provide the fallback when those checks fail. DGTSV factors the tridiagonal system with partial pivoting and solves for one or more right-hand sides, while DGTSVX adds a condition number estimate and iterative refinement with computed error bounds. Numerical linear algebra practice treats speed and reliability as separate concerns: fast no-pivot TDMA when stability is known, pivoted LAPACK routines when it isn’t.
How do you handle cyclic or periodic boundary conditions?
Periodic boundary conditions couple the first and last unknowns directly, adding nonzero entries in the top-right and bottom-left corners of the matrix. That breaks the strict three-band structure TDMA assumes, so the algorithm can’t be applied directly to a cyclic system.
Split off the corner coupling. Rewrite the cyclic matrix as a standard tridiagonal matrix plus a rank-1 correction that captures the corner terms.
Solve two ordinary TDMA systems. One uses the original right-hand side; the second uses a unit vector that represents the corner coupling’s structure.
Combine with the Sherman-Morrison formula. The two partial solutions combine algebraically to recover the exact solution to the original cyclic system, a technique documented for periodic tridiagonal matrices.
Because both auxiliary systems share the same a, b, c coefficients, you can reuse most of the forward sweep work between them. The main stability caveat carries over unchanged: if the underlying tridiagonal part isn’t diagonally dominant, the cyclic correction inherits that instability.
What are the parallel alternatives to serial Thomas?
TDMA’s forward sweep is inherently serial: each row’s pivot depends on the previous row’s result, so a single long system doesn’t parallelize well on its own. The workaround engineers use is batching: when a problem produces many independent small tridiagonal systems, as ADI sweeps and some spline-fitting workflows do, those independent systems can run in parallel even though each one is solved serially.
For a single large system on parallel hardware, several alternative algorithms trade arithmetic work for parallel steps:
Cyclic reduction (CR) eliminates alternating unknowns recursively, trading more total arithmetic for O(log n) parallel steps.
Parallel cyclic reduction (PCR) removes CR’s sequential back-substitution phase at the cost of more work per step.
Recursive doubling (RD) reformulates the recurrence for direct parallel evaluation, useful for very regular problem structures.
Hybrid tiled PCR plus p-Thomas splits large systems into tiles, reduces with PCR or CR, then finishes each tile with a parallel Thomas pass. One GPU study reported speedups up to 28 times over sequential LAPACK for certain workloads, though results depend heavily on hardware and problem shape.
Pro Tip: Treat algorithm choice as a function of workload shape, not raw system size: one long system favors serial TDMA; many systems or a GPU target favor a hybrid CR plus PCR approach.
Where does TDMA fit into applied engineering workflows?
Implicit time-stepping for 1D diffusion or Poisson problems is the canonical use case: each time step’s finite-difference discretization produces a fresh tridiagonal system, solved with TDMA before advancing. In ADI schemes for two-dimensional problems, this repeats once per grid line in each sweep direction, so TDMA often runs thousands of times inside one simulation.
Natural cubic spline interpolation reduces the continuity conditions between segments to a tridiagonal system for the second derivatives at each knot.
Block-tridiagonal systems, common in coupled multi-physics problems, replace each scalar entry with a small matrix, which changes the arithmetic and the conditioning behavior enough that scalar TDMA logic should not be reused without modification.
What’s the quick checklist before shipping TDMA code?
Confirm sentinel values are set correctly: a[0] = 0 and c[n-1] = 0 for zero-based indexing.
Copy b and d before elimination whenever the caller might need the originals afterward.
Run a pivot threshold check and a post-solve residual check, falling back to DGTSV or DGTSVX when either fails.
Keep arithmetic in-place and loops sequential for cache-friendly memory access; batch independent systems rather than trying to parallelize one.
Pro Tip: Write one unit test that constructs a matrix from a known x, generates d by multiplication, and checks that your solver recovers x within floating-point tolerance. It catches more bugs than any amount of manual review.
Author and publisher background on this guide
This guide was written by Joel, who covers numerical methods and applied simulation topics for engineers, including implementation walkthroughs for heat transfer and fluid flow solvers. We build and maintain the Jewlz Engineering Toolkit, which includes thermal analysis, CFD, and pressure vessel simulation modules built on the same discretization methods discussed here, giving engineers a higher-level workflow once the underlying solver choices are settled.
A practical take on choosing between TDMA and its alternatives
Plain TDMA remains the right default for a single serial solve where diagonal dominance holds: it’s shorter to code and faster per system than any pivoted alternative. Reach for a hybrid CR plus PCR approach only when you have many systems or a GPU target that justifies the added complexity. Keep pivot and residual checks in your test suite regardless of which path you choose, and escalate to block-aware methods the moment your problem’s unknowns stop being scalars.
— Joel
An adjacent option when you need full simulation workflows
TDMA libraries solve one piece of a larger problem. When you need the discretization, boundary conditions, and solver working together, our Jewlz Engineering Toolkit combines thermal analysis, CFD, and pressure vessel simulation in one place, with free basic tools and a built-in material property database. It’s not a drop-in replacement for a low-level solver library, but a workflow for engineers who’d rather not assemble one from scratch.
Explore the engineering toolkit to see which modules fit your next simulation project.
FAQ
What problem does the Thomas algorithm actually solve?
It solves a system of linear equations Ax = d where A is tridiagonal, meaning nonzero entries exist only on the main diagonal and its two adjacent bands. This specialized structure lets it run in O(n) time and O(n) storage instead of the cubic cost of general Gaussian elimination.
How does TDMA differ from general Gaussian elimination?
TDMA is Gaussian elimination restricted to a matrix where fill-in outside the three bands never occurs, so the forward sweep and backward substitution only ever touch neighboring entries. That structural guarantee is what keeps the arithmetic linear rather than cubic, but it also means TDMA fails silently if the matrix isn’t actually tridiagonal.
Is the Thomas algorithm always numerically stable?
No. Stability is reliable under diagonal dominance or symmetric positive definiteness, but plain Thomas elimination has no pivoting and no guaranteed stability for an arbitrary tridiagonal matrix. When those conditions aren’t confirmed, a pivoted routine like LAPACK’s DGTSV is the safer choice.
How do I handle periodic boundary conditions with TDMA?
Periodic boundaries add corner couplings that break strict tridiagonal structure, so the standard fix is to solve two ordinary TDMA systems and combine them with the Sherman-Morrison formula. This recovers the exact cyclic solution while reusing most of the forward sweep work between the two systems.
When should I use a parallel method instead of serial TDMA?
A single long tridiagonal system is usually faster with serial TDMA since the method’s own sequential dependency limits parallel speedup. For many independent systems or a GPU target, hybrid approaches combining cyclic reduction with parallel Thomas solves have shown substantial speedups over sequential LAPACK in targeted GPU benchmarks, though results depend on hardware and problem shape.
Sources
Recommended
Comments