Store the road, not the empty landscape只存道路,不存两旁的空地
Imagine a row of rooms exchanging heat only with their immediate neighbors. The coupling matrix has entries near the diagonal and zeros far away. Keeping every zero is like reserving a parking space for every place a car cannot go. A banded layout keeps the useful diagonal strip.想象一排房间只和左右邻居交换热量:对应矩阵的非零元素集中在对角线附近,远处都是零。把所有零也存下来,就像给车根本开不到的地方分配停车位。带状存储只保留有用的那条“道路”。
Our operation is GAXPY: add a matrix–vector product to an existing vector. The initial y matters: replacing it with zero silently changes the problem. All indices below start at zero; the compact array is stored column by column.这里计算 GAXPY:把矩阵向量积加到已有向量上。初始 y 不能忽略,把它清零就改变了问题。下文全部采用从 0 开始的索引,紧凑数组按列存放。
See where each entry goes亲手看一个元素搬到哪里
Choose the lower and upper bandwidths, then select a row and column. The highlighted dense entry and compact slot represent the same value. A dash in the compact array is unused padding, not another matrix entry.调整上下带宽,再选择行和列。高亮的稠密矩阵元素和紧凑槽位代表同一个数。紧凑数组里的短横线是未使用的填充位置,不是额外的矩阵元素。
Turn a diagonal into a row把斜着的对角线摆成一行
The main diagonal has i − j = 0, so it belongs to compact row q. Moving one step below the diagonal increases the compact row by one. Moving above decreases it. A compact column has h = p + q + 1 slots, so the flat offset combines the column start with that row.主对角线上 i − j = 0,因此位于紧凑数组的第 q 行。向下偏一条对角线,紧凑行号加一;向上偏则减一。每个紧凑列占 h = p + q + 1 个槽位,线性地址就是列起点再加行号。
For n = 6, p = 2, q = 1, entry A(3,2) goes to row 2 of compact column 2, hence offset 10. The dense row-major offset would be 20. Never reuse a row-major offset formula for this column-packed array.例如 n = 6、p = 2、q = 1 时,A(3,2) 被放到紧凑数组第 2 列、第 2 行,线性偏移是 10;它在行优先稠密数组中的偏移却是 20。不要把行优先的地址公式直接套到按列压缩的数组上。
These limits clip the band at the matrix edges. They also prevent reading the padding slots. For a tridiagonal matrix, an interior column touches three rows, while the first and last columns touch only two.这组边界把带状区域裁剪到矩阵内部,也避免读取填充槽位。三对角矩阵的内部列涉及三行,第一列和最后一列则只涉及两行。
Run a checked C++ implementation运行带校验的 C++ 实现
This complete C++17 program tests diagonal, asymmetric-band, and full-band cases, including n = 1. It compares against a separately indexed dense product with a nonzero initial y. The packing step deliberately inspects the whole dense matrix; only the subsequent band update has band-linear cost.这个完整 C++17 程序覆盖对角、上下带宽不等和满带等情况,也包含 n = 1。它用不同索引方式的稠密乘法作对照,并保留非零初始 y。转换函数会特意检查整个稠密矩阵;只有转换完成后的带状更新才具有按带宽线性增长的成本。
#include <algorithm>
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <vector>
using Vec = std::vector<double>;
// Zero-based indices. Each packed column occupies p+q+1 slots.
Vec pack_band(const Vec& a, int n, int p, int q) {
if (n <= 0 || p < 0 || q < 0 || p >= n || q >= n ||
a.size() != std::size_t(n) * n)
throw std::invalid_argument("invalid band dimensions");
const int h = p + q + 1;
Vec band(std::size_t(h) * n, 0.0);
for (int j = 0; j < n; ++j)
for (int i = 0; i < n; ++i) {
const bool inside = i >= j-q && i <= j+p;
if (!inside && a[i*n+j] != 0.0)
throw std::invalid_argument("nonzero outside band");
if (inside) band[j*h + q+i-j] = a[i*n+j];
}
return band;
}
void band_gaxpy(const Vec& band, const Vec& x, Vec& y,
int n, int p, int q) {
if (n <= 0 || p < 0 || q < 0 || p >= n || q >= n ||
x.size() != std::size_t(n) || y.size() != std::size_t(n) ||
band.size() != std::size_t(n) * (p+q+1))
throw std::invalid_argument("incompatible dimensions");
const int h = p+q+1;
for (int j = 0; j < n; ++j)
for (int i = std::max(0,j-q); i <= std::min(n-1,j+p); ++i)
y[i] += band[j*h + q+i-j] * x[j];
}
int main() {
double max_error = 0;
int cases = 0;
for (int n : {1, 2, 5, 9})
for (int p = 0; p < n; ++p)
for (int q = 0; q < n; ++q) {
Vec a(n*n,0.0), x(n), y(n,1.0), expected=y;
for (int i=0; i<n; ++i) {
x[i]=i+1;
for (int j=0; j<n; ++j)
if (i>=j-q && i<=j+p) a[i*n+j]=(i==j?2.0:-1.0);
}
for (int i=0; i<n; ++i)
for (int j=0; j<n; ++j) expected[i]+=a[i*n+j]*x[j];
auto band=pack_band(a,n,p,q);
band_gaxpy(band,x,y,n,p,q);
for (int i=0; i<n; ++i)
max_error=std::max(max_error,std::abs(y[i]-expected[i]));
++cases;
if (n==5 && p==1 && q==1) {
std::cout << "tridiagonal y:";
for (double v:y) std::cout << ' ' << v;
std::cout << '\n';
}
}
std::cout << "cases=" << cases << " max_error=" << max_error << '\n';
return max_error < 1e-12 ? 0 : 1;
}
Expected checks: the tridiagonal example prints y = 1, 1, 1, 1, 7; all 111 cases should have zero error for these small integer-valued inputs. Try swapping p and q or removing an edge bound and observe which checks fail.预期结果:三对角示例输出 y = 1、1、1、1、7;这组小整数输入的 111 个测试应全部误差为零。可以尝试交换 p、q,或删去边界裁剪,看看哪些测试会失败。
What is actually saved?到底省下了什么?
For 0 ≤ p,q < n, this is the number of positions inside the band, including any numerical zeros there. The update performs one multiply and one add per position: 2N_band FLOPs. Allocated storage is n(p+q+1) values, slightly larger because of padding. With fixed bandwidth, both grow linearly with n.当 0 ≤ p,q < n 时,这是带内位置的准确个数,带内恰好为零的数也算在内。每个位置做一次乘法和一次加法,合计 2N_band 次浮点运算。分配的存储为 n(p+q+1) 个数,因为填充而稍大。带宽固定时,两者都随 n 线性增长。
For n = 1000 and p = q = 1, dense storage needs 1,000,000 doubles, compact storage 3,000, and the update visits 2,998 positions. These are element counts, not measured speedups. A wide band can erase the benefit; irregular sparsity often needs a different format such as CSR.当 n = 1000、p = q = 1 时,稠密存储需要 1,000,000 个 double,带状存储分配 3,000 个,更新只遍历 2,998 个位置。这些是元素数量,不是实测加速比。带宽很大时优势会消失;不规则稀疏结构通常更适合 CSR 等格式。
Read further继续阅读
LAPACK Users’ Guide: band storage and factorization paddingLAPACK 用户指南:带状存储与分解时的填充
Next: reuse symmetry instead of storing both halves下一篇:利用对称性,不存两份镜像