Thomas algorithm: O(N) solver for tridiagonal systems. More...
#include "config.h"#include "utilities/expected.h"#include "matrix_traits.h"#include "fused/fused_vector.h"#include "math/math_utils.h"
Include dependency graph for tridiagonal.h:
Go to the source code of this file.
Namespaces | |
| namespace | matrix_algorithms |
Functions | |
| template<typename T , my_size_t N> | |
| Expected< FusedVector< T, N >, MatrixStatus > | matrix_algorithms::thomas_solve (const FusedVector< T, N > &a, const FusedVector< T, N > &d, const FusedVector< T, N > &c, const FusedVector< T, N > &b) |
| Solve a tridiagonal system Ax = b using the Thomas algorithm. | |
Detailed Description
Thomas algorithm: O(N) solver for tridiagonal systems.
Solves Ax = b where A is tridiagonal, represented by three vectors:
- a: sub-diagonal (N−1 elements, a(0) unused)
- d: main diagonal (N elements)
- c: super-diagonal (N−1 elements, c(N−1) unused)
[ d(0) c(0) ] [x(0)] [b(0)] [ a(1) d(1) c(1) ] [x(1)] [b(1)] [ a(2) d(2) c(2) ] [x(2)] = [b(2)] [ ... ... ... ] [ ... ] [ ... ] [ a(N-1) d(N-1) ] [x(N-1)] [b(N-1)]
ALGORITHM
Forward sweep (eliminate sub-diagonal): For i = 1 … N−1: w = a(i) / d'(i−1) d'(i) = d(i) − w · c(i−1) b'(i) = b(i) − w · b'(i−1)
Back substitution: x(N−1) = b'(N−1) / d'(N−1) For i = N−2 … 0: x(i) = (b'(i) − c(i) · x(i+1)) / d'(i)
Complexity: O(N) — optimal for tridiagonal systems.
USE CASES
- Cubic spline interpolation
- Trajectory smoothing
- Heat equation / diffusion (implicit methods)
- Boundary value problems (finite differences)
FAILURE MODES
- MatrixStatus::Singular — zero pivot encountered during forward sweep
