The arithmetic is only half the story算得一样多，为什么速度不同？

Think of a cook who needs ingredients from a pantry. Two recipes may use the same number of cuts and stirs, yet one walks to the pantry for every pinch while the other brings a tray. Matrix performance has the same tension: arithmetic matters, but so does how often data must travel.想象厨师从储物间取食材。两份菜谱切、炒的次数一样，一份却每取一撮盐都要走一趟，另一份先把食材装成一盘。矩阵性能也有这样的差别：计算次数重要，数据来回搬运多少次同样重要。

This note follows one operation: update C by AB. A has m rows and k columns, B has k rows and n columns, and C has m rows and n columns. We use row-major C++ arrays throughout; the final section explains how the conclusion changes for column-major storage.本文围绕一个操作：把 AB 加到 C 上。A 是 m 行 k 列，B 是 k 行 n 列，C 是 m 行 n 列。C++ 示例统一采用行优先数组，后面会说明列优先存储会怎样改变结论。

1.1

BLAS levels: opportunities to reuse dataBLAS 层级：数据复用的机会

A level describes the shapes involved, not a speed ranking. Level 1 works mainly with vectors, level 2 with a matrix and vectors, and level 3 with matrices. For square dense problems, the following rough counts explain why matrix multiplication can reuse each loaded value more often.层级描述参与运算的数据形状，不是速度排行榜。Level 1 主要处理向量，Level 2 处理矩阵和向量，Level 3 处理矩阵。对稠密方阵问题，下面的量级说明矩阵乘法为何有机会更充分地复用每次载入的数据。

Level / example层级 / 示例Stored data存储数据量Arithmetic算术量
1 · xᵀyO(n)O(n)
2 · y ← y + AxO(n²)O(n²)
3 · C ← C + ABO(n²)O(n³)

These counts are not memory-traffic measurements. A poorly ordered matrix product can load the same values repeatedly. Structure also changes the problem: diagonal matrices skip most entries, banded matrices skip distant interactions, and symmetric matrices avoid storing both halves. Exploit a valid structural property before optimizing a dense loop.这些量级不是实际内存流量。循环安排不当的矩阵乘法可能反复载入同一批数。矩阵结构还会改变问题本身：对角矩阵跳过大多数元素，带状矩阵跳过远距离作用，对称矩阵不必保存两份镜像。先利用确实成立的结构，再优化稠密循环。

Read addresses before reading timings先看访问地址，再看计时

1. Stride: follow the inner loop1. 步长：跟着最内层循环走

In row-major storage, moving across a row advances one element; moving down a column jumps by the row length. With i-j-t order, the innermost loop reads B(t,j) down a column. With i-t-j order, it reads B(t,j) across a row while reusing one value A(i,t). C also becomes a contiguous row update.行优先存储中，横着走一列只前进一个元素，竖着走一行却要跳过整行长度。按 i-j-t 顺序循环，最内层沿 B(t,j) 的一列读取；按 i-t-j 顺序循环，则沿一行读取，并反复使用同一个 A(i,t)。C 也变成连续的行更新。

3.1
B · numbers are flat element offsets, not matrix valuesB · 数字是线性元素偏移，不是矩阵数值

The display assumes eight doubles per aligned 64-byte cache line and an 8 × 8 row-major B. Eight row reads fall in one line; eight column reads touch eight different lines. This counts distinct lines in a toy access trace, not actual cache misses or a measured speedup. Cache capacity, prior accesses, prefetching, and alignment all matter.图中假设一个对齐的 64 字节缓存行容纳八个 double，B 为 8 × 8 行优先矩阵。横向读八次落在同一缓存行，纵向读八次涉及八个不同缓存行。这只是简化访问序列中涉及的缓存行计数，不是真实缓存未命中次数，也不是实测加速比；容量、历史访问、预取与对齐都会影响结果。

2. Chunking is not a promise of SIMD2. 分段循环不等于保证 SIMD

Split a long vector into fixed-size chunks and handle the final short chunk explicitly. That explains tail processing without depending on a CPU-specific intrinsic. A compiler may vectorize the inner loop, but source-level chunking alone proves nothing about the generated instructions; inspect the assembly or a vectorization report.把长向量切成固定大小的小段，并显式处理最后不足一段的尾部，就能理解尾部处理，而不依赖某种 CPU 的专用指令。编译器可能向量化内层循环，但仅凭源码分段不能证明实际生成了 SIMD 指令，应查看汇编或向量化报告。

An older vector-pipeline model assigns a startup cost τ_add to a vector add and τ_data to each load or store, followed by one element per cycle. If n is divisible by vector length v, its arithmetic and data-motion counts give the ratio below. The denominator must use τ_data. On modern processors, overlap and cache behavior make this a teaching model rather than a timing predictor.旧式向量流水线模型给向量加法分配启动代价 τ_add，给每次载入或存储分配 τ_data，随后每周期处理一个元素。若 n 能被向量长度 v 整除，就能得到下面的算术与数据搬运周期比，分母必须使用 τ_data。现代处理器存在执行重叠和复杂缓存行为，因此这只是教学模型，不能直接预测耗时。

3.2

3. GAXPY and outer products move different data3. GAXPY 与外积，搬的是不同的数据

3.3

For square n × n data, both do about 2n² FLOPs. GAXPY reads A and updates a vector, whereas an outer-product update reads and writes A itself. In a simplified model where whole vectors fit in registers and x and y remain resident, GAXPY needs n+3 vector transfers and an outer update 2n+2. The ratio approaches two; it is not exactly two for finite n and says nothing by itself about elapsed time.对 n × n 方阵，两者都约做 2n² 次 FLOP。GAXPY 读取 A、更新向量；外积更新却要读写矩阵 A 本身。在整条向量能装入寄存器、且 x、y 驻留其中的简化模型里，GAXPY 需要 n+3 次整向量传输，外积更新需要 2n+2 次。两者比值趋近二，但有限 n 时并不恰好为二，也不能单凭这个比值判断耗时。

The small program below makes both operations and the vector tail concrete. The vector has five elements and a chunk size of four. Output checks prevent a missing tail or an accidental overwrite from looking like a successful optimization.下面的小程序把这两种运算和向量尾部都写成可运行的例子。向量长度为五，分段长度为四。输出校验会防止漏掉尾部或误把累加写成覆盖，却仍被当作成功的优化。

#include <algorithm>
#include <iostream>
#include <vector>
using Vec=std::vector<double>;
int main() {
    const Vec x={1,2,3,4,5},y={5,4,3,2,1};
    Vec z(x.size());
    constexpr std::size_t chunk=4;
    for(std::size_t first=0;first<x.size();first+=chunk)
        for(std::size_t i=first;i<std::min(first+chunk,x.size());++i)
            z[i]=x[i]+y[i]; // The final chunk contains just one element.
    std::cout << "chunked sum:";
    for(double v:z) std::cout << ' ' << v;
    // Row-major 2x2 examples with nonzero initial outputs.
    Vec a={1,2,3,4},v={2,3},out={10,20},outer=a;
    for(int i=0;i<2;++i) for(int j=0;j<2;++j) out[i]+=a[i*2+j]*v[j];
    for(int i=0;i<2;++i) for(int j=0;j<2;++j) outer[i*2+j]+=v[i]*v[j];
    std::cout << "\ngaxpy: " << out[0] << ' ' << out[1] << "\nouter:";
    for(double value:outer) std::cout << ' ' << value;
    std::cout << '\n';
    return z==Vec(5,6) && out==Vec{18,38} && outer==Vec{5,8,9,13} ? 0 : 1;
}

4. Blocking: keep a working set nearby4. 分块：把一盘食材留在手边

Divide A, B, and C into tiles. For one C tile, walk through the matching A and B tiles and accumulate partial products. For three square b × b tiles of doubles, the nominal footprint is 24b² bytes. Choose b with room for other working data; nominal capacity alone does not account for associativity, registers, packing, or multiple threads.把 A、B、C 划成小块。固定一块 C，依次取对应的 A、B 小块累加部分乘积。三个 b × b 的 double 方块名义上占 24b² 字节。选择 b 时还要给其他工作数据留空间；只看缓存标称容量，并没有考虑组相联、寄存器、打包和多线程的影响。

3.4

Our blocked example still uses row-major storage; it does not magically change the physical layout into tiles. Every tile-end loop uses min to handle rectangular matrices and dimensions that are not multiples of the tile size. Production matrix multiplication adds microkernels, packing, vector instructions, and often threading.下方分块示例仍然采用行优先存储，并没有自动把物理布局改成块布局。每个块末尾都用 min 裁剪，以支持长方形矩阵和不能整除块长的尺寸。生产级矩阵乘法还会使用微内核、打包、向量指令，并经常配合多线程。

Compile, check, then measure先编译校验，再做计时

The complete C++17 program compares i-j-t, i-t-j, and a 16 × 16 blocked traversal. First it checks small rectangular cases with nonzero initial C. Then it warms up each kernel twice and records nine samples in rotating order, reporting median and quartiles. Initialization, copying, and result checks stay outside each timed interval.这个完整 C++17 程序对比 i-j-t、i-t-j 和 16 × 16 分块遍历。先用非零初始 C 检查小型长方形案例，再对每个内核预热两次，按轮换顺序收集九次采样，并报告中位数和四分位数。初始化、复制和结果校验都在计时区间之外。

#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <vector>
using Vec=std::vector<double>;
// Row-major A(m,k), B(k,n), C(m,n). C is UPDATED, not overwritten.
void ijk(const Vec& a,const Vec& b,Vec& c,int m,int n,int k) {
    for(int i=0;i<m;++i) for(int j=0;j<n;++j) for(int t=0;t<k;++t)
        c[i*n+j]+=a[i*k+t]*b[t*n+j];
}
void ikj(const Vec& a,const Vec& b,Vec& c,int m,int n,int k) {
    for(int i=0;i<m;++i) for(int t=0;t<k;++t) {
        const double value=a[i*k+t];
        for(int j=0;j<n;++j) c[i*n+j]+=value*b[t*n+j];
    }
}
void blocked(const Vec& a,const Vec& b,Vec& c,int m,int n,int k) {
    constexpr int tile=16;
    for(int ii=0;ii<m;ii+=tile) for(int jj=0;jj<n;jj+=tile)
        for(int tt=0;tt<k;tt+=tile)
            for(int i=ii;i<std::min(ii+tile,m);++i)
                for(int t=tt;t<std::min(tt+tile,k);++t) {
                    const double value=a[i*k+t];
                    for(int j=jj;j<std::min(jj+tile,n);++j)
                        c[i*n+j]+=value*b[t*n+j];
                }
}
double error(const Vec& a,const Vec& b) {
    double result=0;
    for(std::size_t i=0;i<a.size();++i) result=std::max(result,std::abs(a[i]-b[i]));
    return result;
}
int main() {
    using Kernel=void(*)(const Vec&,const Vec&,Vec&,int,int,int);
    const std::array<Kernel,3> kernels={ijk,ikj,blocked};
    const std::array<const char*,3> names={"ijk","ikj","blocked"};
    double max_error=0;
    // Rectangular sizes and partial tiles catch loop-bound mistakes.
    for(auto shape : {std::array<int,3>{1,1,1},{3,5,2},{17,19,13}}) {
        const auto [m,n,k]=shape;
        Vec a(m*k),b(k*n),reference(m*n,1.0);
        for(int i=0;i<m*k;++i) a[i]=(i%7)-3;
        for(int i=0;i<k*n;++i) b[i]=(i%5)-2;
        ijk(a,b,reference,m,n,k);
        for(auto kernel:kernels) {
            Vec actual(m*n,1.0);kernel(a,b,actual,m,n,k);
            max_error=std::max(max_error,error(actual,reference));
        }
    }
    // Teaching microbenchmark. Not a claim about all matrices or computers.
    const int m=127,n=113,k=95;
    Vec a(m*k),b(k*n),initial(m*n,1.0),reference=initial,c=initial;
    for(int i=0;i<m*k;++i) a[i]=std::sin(0.01*i);
    for(int i=0;i<k*n;++i) b[i]=std::cos(0.02*i);
    ijk(a,b,reference,m,n,k);
    std::array<std::vector<double>,3> samples;
    for(int warmup=0;warmup<2;++warmup) for(auto kernel:kernels) {
        c=initial;kernel(a,b,c,m,n,k);
        max_error=std::max(max_error,error(c,reference));
    }
    double checksum=0;
    for(int round=0;round<9;++round) for(int offset=0;offset<3;++offset) {
        const int which=(round+offset)%3;
        c=initial; // Reset and allocation are outside the timed interval.
        const auto start=std::chrono::steady_clock::now();
        kernels[which](a,b,c,m,n,k);
        const auto stop=std::chrono::steady_clock::now();
        samples[which].push_back(std::chrono::duration<double,std::milli>(stop-start).count());
        max_error=std::max(max_error,error(c,reference));
        checksum+=std::accumulate(c.begin(),c.end(),0.0);
    }
    std::cout << std::setprecision(8) << "shape=" << m << ',' << n << ',' << k
              << " max_error=" << max_error << " checksum=" << checksum << '\n';
    for(int i=0;i<3;++i) {
        auto& s=samples[i];std::sort(s.begin(),s.end());
        std::cout << names[i] << " median_ms=" << s[4]
                  << " Q1=" << s[2] << " Q3=" << s[6] << '\n';
    }
    return max_error<1e-9 && std::isfinite(checksum) ? 0 : 1;
}

Read max_error before comparing the timings. It should be below 10⁻⁹ for these inputs. The checksum makes computed values observable. No fixed speedup is promised: try tile sizes 8, 16, and 32, or change the matrix shape. A slower blocked result is a valid observation, not a broken experiment.比较计时前先看 max_error：这组输入应小于 10⁻⁹。checksum 让计算结果可被外部观察。这里不承诺固定加速比：可以把块长改为 8、16、32，或改变矩阵形状。分块版本更慢也是有效观测，不代表实验坏了。

Take the rule back to the layout把规则还给具体的存储布局

For column-major arrays, the first index is contiguous instead, and j-t-i is a natural update order. The transferable rule is to follow contiguous data and preserve useful reuse, not to memorize one universally best loop order. Faster arithmetic cannot rescue an invalid index, a missing tail, or a changed mathematical operation.若数组采用列优先，连续的是第一个索引，j-t-i 就是自然的更新顺序。真正可以迁移的规则是顺着连续数据访问、保留有用的复用，而不是背下某一种永远最快的循环顺序。再快的算术也救不了错误索引、遗漏尾部或被改变的数学运算。

Read further继续阅读

Netlib BLAS: operation levels and reference routinesNetlib BLAS：运算层级与参考例程

Benchmark with evidence: repeat, verify, report conditions可信基准测试：重复、校验、记录条件

Exploit a narrow band before optimizing dense work优化稠密运算前，先利用窄带结构