- Generated by
1.9.8
|
numerics
|
Banded matrices arise naturally in the discretization of differential equations, particularly in:
A banded matrix \(A \in \mathbb{R}^{n \times n}\) has bandwidth parameters \((k_l, k_u)\) where:
The element \(a_{ij} = 0\) whenever \(j < i - k_l\) or \(j > i + k_u\).
Structure visualization (pentadiagonal, \(k_l = k_u = 2\)):
\[A = \begin{pmatrix} \times & \times & \times & 0 & 0 & \cdots \\ \times & \times & \times & \times & 0 & \cdots \\ \times & \times & \times & \times & \times & \cdots \\ 0 & \times & \times & \times & \times & \cdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots \end{pmatrix}\]
Key advantage: Exploiting band structure reduces:
The standard format for band storage, used in LAPACK's DGBTRF/DGBTRS, stores diagonals in rows of a 2D array with leading dimension \(\text{ldab} = 2k_l + k_u + 1\).
Layout: For matrix element \(a_{ij}\) within the band:
\[\text{AB}(k_l + k_u + i - j, \, j) = a_{ij}\]
where AB is the band storage array with 0-based indexing.
Visual representation (for \(n=5\), \(k_l=2\), \(k_u=1\)):
**Why extra \(k_l\) rows?** During LU factorization with row pivoting, fill-in can occur in the upper triangle up to \(k_l\) positions above the original upper bandwidth.
Column-major storage (as in the layout above) ensures that:
**Memory access pattern during elimination of column \(j\)**:
For a banded matrix \(A\) with bandwidths \((k_l, k_u)\), we compute \(PA = LU\) where:
Theorem (Band Preservation): If \(A\) has lower bandwidth \(k_l\) and upper bandwidth \(k_u\), then:
Proof: During Gaussian elimination of column \(j\):
With pivoting, row \(j\) may come from row \(j + k_l\) at worst, creating fill-in up to column \((j + k_l) + k_u\) for original row \(j\). But this is stored in the extra \(k_l\) rows we allocated. \(\square\)
Standard Gaussian elimination adapted for band structure:
For column \(j = 0, 1, \ldots, n-1\):
Given \(PA = LU\), solve \(A\mathbf{x} = \mathbf{b}\) in three steps:
Step 1: Apply permutation ( \(\mathbf{b} \leftarrow P\mathbf{b}\))
Step 2: Forward substitution (solve \(L\mathbf{y} = P\mathbf{b}\))
Step 3: Back substitution (solve \(U\mathbf{x} = \mathbf{y}\))
LU Factorization:
For column \(j\), the number of operations is:
Summing over all columns:
\[\text{FLOPs} = \sum_{j=0}^{n-1} O(k_l \cdot k_u) = O(n \cdot k_l \cdot k_u)\]
More precisely, accounting for boundary effects:
\[\text{FLOPs} \approx n \cdot k_l \cdot (k_l + k_u) + O(k_l^2 \cdot k_u)\]
Forward/Back Substitution:
| Method | Factor | Solve | Storage |
|---|---|---|---|
| Dense LU | \(\frac{2}{3}n^3\) | \(2n^2\) | \(n^2\) |
| Banded LU | \(n \cdot k_l \cdot (k_l + k_u)\) | \(2n(k_l + k_u)\) | \((2k_l + k_u + 1) \cdot n\) |
| Tridiagonal (Thomas) | \(8n\) | \(5n\) | \(4n\) |
Example: \(n = 10000\), pentadiagonal ( \(k_l = k_u = 2\))
For large \(n\), the algorithm is typically memory-bound rather than compute-bound.
Bytes per FLOP (banded factorization):
For pentadiagonal ( \(k_l = k_u = 2\)): \(\frac{8 \times 7}{2 \times 4} = 7\) bytes/FLOP
Modern CPUs achieve ~50 GFLOPS but only ~50 GB/s memory bandwidth. At 7 bytes/FLOP, we're limited to ~7 GFLOPS — memory bound.
Theorem: LU factorization with partial pivoting satisfies:
\[\|L\|_\infty \leq 1 + k_l\]
for banded matrices, where the bound comes from at most \(k_l\) nonzeros per column of \(L\), each with magnitude \(\leq 1\).
Growth Factor: The growth factor \(\rho = \max_{ij}|u_{ij}| / \max_{ij}|a_{ij}|\) is bounded:
\[\rho \leq 2^{k_l}\]
This is much better than the \(2^{n-1}\) worst-case for dense matrices.
The computed solution \(\hat{\mathbf{x}}\) satisfies:
\[\frac{\|\hat{\mathbf{x}} - \mathbf{x}\|}{\|\mathbf{x}\|} \lesssim \kappa(A) \cdot \epsilon_{\text{mach}}\]
where \(\kappa(A) = \|A\| \cdot \|A^{-1}\|\) is the condition number.
Condition number estimation for banded matrices:
Definition: \(A\) is diagonally dominant if:
\[|a_{ii}| \geq \sum_{j \neq i} |a_{ij}|, \quad \forall i\]
with strict inequality for at least one row.
Theorem: For diagonally dominant banded matrices:
Common diagonally dominant systems:
Element access mapping:
The inner loops in elimination and solve are vectorizable:
The #pragma omp simd directive hints to the compiler that iterations are independent and can use SIMD instructions (AVX-512 on modern Intel/AMD).
Blocking for L2 cache: For very wide bands ( \(k_l, k_u > 100\)), block the elimination:
When solving with multiple RHS vectors (common in radiative transfer for multiple wavelengths):
Parallelization over RHS is embarrassingly parallel — no synchronization needed.
For GPU efficiency, we solve many independent banded systems in parallel:
| Scenario | CPU | GPU | Winner |
|---|---|---|---|
| Single small system ( \(n < 1000\)) | 10 μs | 50 μs (launch overhead) | CPU |
| Single large system ( \(n > 10000\)) | 1 ms | 0.5 ms | GPU |
| 1000 small systems | 10 ms | 0.1 ms | GPU |
| 1000 large systems | 1 s | 50 ms | GPU |
The key insight: GPU excels at batched solves, which is exactly the pattern in radiative transfer (one system per spectral band or grid column).
The radiative transfer equation for diffuse radiation:
\[\mu \frac{dI}{d\tau} = I - S\]
In the two-stream approximation (up/down fluxes \(F^+\), \(F^-\)), discretized over \(N\) atmospheric layers:
\[\begin{pmatrix} T_1 & R_1 & & \\ R_1 & T_1 & R_2 & \\ & \ddots & \ddots & \ddots \\ & & R_{N-1} & T_N \end{pmatrix} \begin{pmatrix} F^+_1 \\ F^-_1 \\ \vdots \\ F^-_N \end{pmatrix} = \begin{pmatrix} S_1 \\ S_2 \\ \vdots \\ S_N \end{pmatrix}\]
where:
Matrix structure: Block tridiagonal with \(2 \times 2\) blocks → banded with \(k_l = k_u = 2\).
TUVX (Tropospheric Ultraviolet-Visible) computes photolysis rate coefficients:
\[J_i = \int_\lambda \sigma_i(\lambda) \phi_i(\lambda) F(\lambda) \, d\lambda\]
where \(F(\lambda)\) is the actinic flux from the radiative transfer solve.
Computational pattern:
This is ideal for batched banded solvers!
NCAR's Derecho supercomputer specifications:
Optimization for Derecho:
Key flags:
-march=znver3: Zen 3 (Milan) specific instructions including AVX2-fopenmp: Enable OpenMP for multi-RHS parallelization-ffast-math: Allow reordering for SIMD (safe for banded solvers)Expected performance:
| System Type | Recommended Solver |
|---|---|
| Tridiagonal, single | Thomas algorithm |
| Tridiagonal, batched | GPU batched Thomas |
| General banded, single | banded_solve() |
| General banded, reuse matrix | banded_lu() + banded_lu_solve() |
| General banded, multi-RHS | banded_lu_solve_multi() |
| Very large single system | Consider iterative (CG with banded preconditioner) |
banded_rcond() after factorization| Function | Purpose | Complexity |
|---|---|---|
BandedMatrix(n, kl, ku) | Construct band storage | \(O(n \cdot \text{ldab})\) |
banded_lu(A, ipiv) | LU factorization | \(O(n \cdot k_l \cdot (k_l + k_u))\) |
banded_lu_solve(A, ipiv, b) | Solve with factored matrix | \(O(n \cdot (k_l + k_u))\) |
banded_lu_solve_multi(A, ipiv, B, nrhs) | Multiple RHS solve | \(O(n \cdot (k_l + k_u) \cdot \text{nrhs})\) |
banded_solve(A, b, x) | One-shot solve | Factor + Solve |
banded_matvec(A, x, y) | Matrix-vector product | \(O(n \cdot (k_l + k_u))\) |
banded_norm1(A) | 1-norm | \(O(n \cdot (k_l + k_u))\) |
banded_rcond(A, ipiv, anorm) | Condition estimate | \(O(n \cdot (k_l + k_u))\) |
Claim: With partial pivoting, \(U\) has at most \(k_l + k_u\) super-diagonals.
Proof:
Consider eliminating column \(j\). The pivot row \(p\) satisfies \(j \leq p \leq j + k_l\).
After swapping row \(p\) into row \(j\):
Since we're in row \(j\) (the pivot row), this creates fill-in up to column \(j + k_l + k_u\).
Thus \(U\) has upper bandwidth \(k_l + k_u\). \(\square\)
Claim: If \(A\) is strictly diagonally dominant, so is the Schur complement during elimination.
Proof:
After eliminating column 0, the \((1,1)\) element of the Schur complement is:
\[\tilde{a}_{11} = a_{11} - \frac{a_{10}}{a_{00}} a_{01}\]
The diagonal dominance of row 1 gives:
\[|a_{11}| > |a_{10}| + |a_{12}| + \cdots\]
We need to show \(|\tilde{a}_{11}| > |\tilde{a}_{12}| + \cdots\) where \(\tilde{a}_{1k} = a_{1k} - \frac{a_{10}}{a_{00}} a_{0k}\).
Using triangle inequality and the fact that \(|a_{10}/a_{00}| < 1\) (from diagonal dominance of row 0):
\[|\tilde{a}_{11}| \geq |a_{11}| - \left|\frac{a_{10}}{a_{00}}\right| |a_{01}|\]
and
\[|\tilde{a}_{1k}| \leq |a_{1k}| + \left|\frac{a_{10}}{a_{00}}\right| |a_{0k}|\]
Summing over \(k \neq 1\) and using diagonal dominance of rows 0 and 1:
\[\sum_{k \neq 1} |\tilde{a}_{1k}| < \sum_{k \neq 1} |a_{1k}| + \left|\frac{a_{10}}{a_{00}}\right| \sum_{k \neq 0} |a_{0k}|\]
\[< |a_{11}| - |a_{10}| + \left|\frac{a_{10}}{a_{00}}\right| \cdot |a_{00}|\]
\[= |a_{11}|\]
But we also have \(|\tilde{a}_{11}| > |a_{11}| - |a_{10}||a_{01}|/|a_{00}| > |a_{11}| - |a_{10}|\) since \(|a_{01}| < |a_{00}|\).
Combining these bounds shows \(|\tilde{a}_{11}| > \sum_{k \neq 1} |\tilde{a}_{1k}|\). \(\square\)
This explains why Thomas algorithm (no pivoting) is stable for diagonally dominant tridiagonal systems.
1.9.8