One value, two positions一个数，对应两个位置
A symmetric matrix is a mirror about its main diagonal. If a spring links objects i and j with the same coupling in both directions, the two matrix entries repeat one piece of information. We can store it once and use it twice. This saves matrix storage; it does not remove either contribution to the product.对称矩阵就像沿主对角线放了一面镜子。如果一根弹簧连接物体 i 和 j，两个方向具有相同的耦合系数，那么矩阵里的两个位置重复记录了同一条信息。我们可以只存一次、使用两次。这省下的是矩阵存储，不是把乘积里的一份贡献删掉。
The discussion and C++ implementation below concern real symmetric matrices. For Hermitian matrices the reflected entry must be conjugated; for skew-symmetric matrices it changes sign and the diagonal is zero. Those cases cannot all use the same real-symmetric update unchanged.下文讲解和 C++ 实现针对实对称矩阵。Hermitian 矩阵的镜像元素需要取共轭；反对称矩阵的镜像元素需要变号，且对角线为零。不能把这些情况都直接套进同一份实对称更新代码。
Follow an entry into the packed array跟着一个元素走进压缩数组
Select any position in the 4 × 4 matrix. Its mirror and its unique packed slot light up together. We store the lower triangle by columns: first four values, then three, then two, then one. Equal colors mean shared storage, not an extra copy.选择 4 × 4 矩阵中的任意位置，它的镜像位置和唯一的压缩槽位会一起亮起。这里按列存下三角：先存四个数，再存三个、两个、一个。同色表示共用存储，不是多存一份。
Count the columns before computing the address先数前面有多少列，再算地址
For zero-based column j, the previous columns contain n + (n−1) + … + (n−j+1) values. Within column j, row i is i−j steps from the diagonal. Adding these two counts gives an integer offset; no floating-point arithmetic belongs in this index calculation.对从 0 开始的列号 j，前面的列一共存了 n + (n−1) + … + (n−j+1) 个数。在第 j 列内，第 i 行距离对角线 i−j 步。两者相加就得到整数偏移，索引计算不需要浮点运算。
If i < j, swap the indices before using this formula. For n = 4, A(1,2) is read as A(2,1): base(1) = 4 and offset = 5. Ten slots replace sixteen, and the fraction stored approaches one half as n grows.若 i < j，先交换行列再套公式。例如 n = 4 时，A(1,2) 转为读取 A(2,1)：base(1) = 4，最终偏移是 5。十个槽位代替十六个；随着 n 增大，所需存储比例趋近一半。
Read once, update both outputs读一次，同时更新两个输出
For an off-diagonal entry, the mirrored matrix position multiplies a different input component, so both updates are necessary. On the diagonal, there is only one contribution. Accidentally applying the two-update rule there doubles the diagonal term.对于非对角元素，镜像位置乘的是另一个输入分量，因此两次更新都不能省。对角线上只有一份贡献，如果把两次更新的规则也用于对角线，就会把该项翻倍。
With A = [[2,1,0],[1,3,2],[0,2,4]], x = [1,2,3], and initial y = [10,20,30], the contribution A(2,1) = 2 adds 4 to y₂ and 6 to y₁. The full result is [14,33,46]. Use this small case to reason about the code before running larger checks.取 A = [[2,1,0],[1,3,2],[0,2,4]]、x = [1,2,3]、初始 y = [10,20,30]。其中 A(2,1) = 2 会给 y₂ 加 4、给 y₁ 加 6。完整结果为 [14,33,46]。先用这个小例子想清楚代码，再做更大的校验。
Run and modify the C++ implementation运行并修改 C++ 实现
#include <algorithm>
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <vector>
using Vec = std::vector<double>;
// Lower triangle, packed by columns. Requires 0 <= j <= i < n.
std::size_t offset(std::size_t i, std::size_t j, std::size_t n) {
return j*(2*n-j+1)/2 + (i-j);
}
Vec pack_lower(const Vec& a, std::size_t n) {
if (a.size()!=n*n) throw std::invalid_argument("matrix size");
Vec ap(n*(n+1)/2);
for (std::size_t j=0; j<n; ++j)
for (std::size_t i=j; i<n; ++i) {
if (a[i*n+j]!=a[j*n+i]) throw std::invalid_argument("not symmetric");
ap[offset(i,j,n)]=a[i*n+j];
}
return ap;
}
void symmetric_gaxpy(const Vec& ap, const Vec& x, Vec& y) {
const std::size_t n=x.size();
if (y.size()!=n || ap.size()!=n*(n+1)/2)
throw std::invalid_argument("incompatible dimensions");
for (std::size_t j=0; j<n; ++j) {
y[j]+=ap[offset(j,j,n)]*x[j]; // The diagonal contributes once.
for (std::size_t i=j+1; i<n; ++i) {
const double a=ap[offset(i,j,n)];
y[i]+=a*x[j];
y[j]+=a*x[i]; // Mirror contribution; no second matrix load.
}
}
}
int main() {
const Vec a={2,1,0, 1,3,2, 0,2,4}, x={1,2,3};
Vec y={10,20,30};
symmetric_gaxpy(pack_lower(a,3),x,y);
std::cout << "y:";
for (double v:y) std::cout << ' ' << v;
std::cout << '\n';
double max_error=0;
for (std::size_t n : {1u,4u,7u}) {
Vec dense(n*n), input(n), actual(n,2.0), expected=actual;
for (std::size_t i=0; i<n; ++i) {
input[i]=double(i)-2.0;
for (std::size_t j=0; j<n; ++j)
dense[i*n+j]=double(i+j+1)+(i==j?2.0:0.0);
}
for (std::size_t i=0; i<n; ++i)
for (std::size_t j=0; j<n; ++j) expected[i]+=dense[i*n+j]*input[j];
symmetric_gaxpy(pack_lower(dense,n),input,actual);
for (std::size_t i=0; i<n; ++i)
max_error=std::max(max_error,std::abs(actual[i]-expected[i]));
}
std::cout << "max_error=" << max_error << '\n';
return y==Vec{14,33,46} && max_error<1e-12 ? 0 : 1;
}
The program checks symmetry when packing, validates buffer sizes, and compares packed multiplication against a dense reference for n = 1, 4, and 7. It should print y: 14 33 46 and max_error=0. Try deleting the mirror update; the result will no longer match.程序打包时检查对称性，验证缓冲区大小，并在 n = 1、4、7 时比较压缩乘法与稠密参考结果。应输出 y: 14 33 46 和 max_error=0。试着删掉镜像更新，结果就无法通过对照校验。
Half the storage is not half the work存储近乎减半，算术量并未减半
Each diagonal entry costs one multiply-add pair; every stored off-diagonal entry costs two. The arithmetic count is therefore the same as dense GAXPY under the usual two-FLOP counting convention. Reading fewer matrix values can still help when memory traffic dominates.每个对角元素做一对乘加，每个已存的非对角元素做两对。按一次乘法加一次加法算两次 FLOP 的约定，算术量与稠密 GAXPY 相同。当访存占主导时，少读取一份矩阵值仍可能有帮助。
Read further继续阅读
LAPACK Users’ Guide: packed storage conventionsLAPACK 用户指南：压缩存储约定