高性能数学库优化指南：矩阵乘法从310ms到9ms的工程实践

1. 为什么“同样一份数学代码”别人的库能快十倍做数学库这件事最早我是很不以为然的。上学那会儿写矩阵乘法三层循环套上去能跑出结果就算完事。直到第一次把一份 N2048 的矩阵乘法丢进项目里发现耗时整整两秒多而同门的程序用同样的算法只花了不到一百毫秒。那一刻我才意识到一个问题一个数学库的性能和你“会用公式”之间隔着一整座冰山。1.1 先厘清“高性能数学库”到底在解决什么问题所谓高性能数学库它不是指把 sin、cos、sqrt 这些函数重新发明一遍而是指在一个特定硬件平台上把数学运算的吞吐量、延迟、精度和资源占用压到极致。常见的形态包括 BLAS 级别的向量与矩阵运算、LAPACK 级别的线性方程组求解、FFT、稀疏矩阵运算以及统计与随机数生成。它们共同的特点是计算模式高度规则、数据规模大、可并行性强因此非常吃“工程优化”而不只是算法复杂度。拿矩阵乘法举例经典的 i-k-j 三重循环时间复杂度都是 O(n^3)这是算法层面的天花板。但同样的 O(n^3)不同写法在真实机器上的耗时差距能拉到 20 倍以上。原因是现代 CPU 的算力远远超过内存系统能供给的数据速度你的代码能不能让计算单元一直有数据可算决定了最终效率。高性能数学库的本质工作就是“喂养”计算单元而不是炫数学技巧。1.2 性能瓶颈的第一性原理访存 计算很多刚接触优化的同学会把注意力放在“减少乘法次数”上比如用斯特拉森算法。但真实项目中斯特拉森很少有人直接用原因很简单它的常数因子大、数值稳定性差更重要的是对现代 CPU 来说乘法本身已经是极廉价的指令了真正贵的是把数据从内存搬进寄存器的那一下。一次 L1 cache miss 大概是 4 个时钟周期L2 miss 约 12 个L3 miss 约 40 个而一次主存访问要 200 个时钟周期以上。相比之下一次浮点乘法在 Haswell 以后的 CPU 上只需要 0.5 个周期左右的吞吐。算一下就知道省一个乘法指令远不如省一次 cache miss 来得实在。所以我在优化任何数学函数之前第一件事就是计算访存强度也就是“每访问一次内存能支撑多少次浮点运算”。对矩阵乘法这种访存密集型任务如果能把数据块塞进 L2/L1 cache让同一块数据被反复复用性能自然就上来了。这个思路贯穿了后面所有优化步骤。1.3 一次真实对比三种矩阵乘法实现的时间差异为了让你快速感知差距我放一张我在自己机器上测过的数据Intel i5-12400单线程N1024实现方式耗时相对速度朴素 i-j-k 三重循环约 310ms1x矩阵分块 转置优化约 76ms约 4x分块 AVX2 向量化 多线程4线程约 9ms约 34x同一个算法同一个精度float只是让数据的流动方式和计算指令的组织方式不同结果差了 34 倍。这不是魔法而是下面每一层优化叠加起来的效果。接下来我会把每一步拆开讲并且把为什么这样做、踩过什么坑都交代清楚。2. 数据布局与内存子系统数学库最快的那条腿如果你只记住一句话那就是数学库性能优化的第一步是数据布局不是 CPU 指令。很多人在这一层还没做好就急着上 AVX、上多线程结果收益很有限因为瓶颈在内存搬运不在计算。2.1 行主序与列主序一次cache miss就足以改变命运在 C/C 里二维数组默认是行主序存储也就是 a[i][j] 在内存里是连续的一行一行排下来的。而 Fortran 和 MATLAB 的默认存储是列主序a[i][j] 的 i 变化最快j 变化最慢。这带来一个直接影响访问 a[i][k] 时按行读是连续的、高速的但如果你在三重循环里对 b[k][j] 按列去读那每一次读取都大概率是 cache miss。我最早写朴素矩阵乘法的时候就是这样内层循环对 b 做列访问每一列都跨越整个矩阵的内存范围几乎全 miss。这 310ms 就是这么来的。解决办法有两个方向一是把 b 矩阵转置成行主序的临时矩阵二是交换循环顺序让内层循环对 a、b 都是行访问。我实际项目中两个都用优先用转置临时矩阵因为它不需要改动计算逻辑只是多一次 O(n^2) 的拷贝对 n1024 来说这个拷贝成本远低于列访问多出来的 cache miss。2.2 AoS/SoA选择线性代数库如何为SIMD铺路矩阵是天生适合行主序的但很多数学库不只是处理矩阵还需要处理“一批点”的坐标变换、批量向量的运算。这时候就涉及到 AoSArray of Structures结构体数组和 SoAStructure of Arrays数组结构体的选择。举个例子处理 10000 个三维点的旋转如果定义成struct Point { float x, y, z; }; Point points[10000];这样声明points[0].x、points[0].y、points[0].z 在内存里是紧挨着的。你要对每个坐标做乘法访问模式是“隔着 12 字节取一个 float”这对 SIMD 非常不友好想要一次加载 4 个 float必须做 gather 或者手动拆包。改成 SoA 就清爽了struct Points { float x[10000]; float y[10000]; float z[10000]; };这样 x 数组是一整块连续内存一次可以 load 4 个甚至 8 个 float直接进入向量寄存器做运算效率提升非常明显。代价是代码可读性略差抽象层级要设计好。我在自己的小库里对外暴露面向对象的接口内部存储全部用 SoA。2.3 对齐分配与跨平台陷阱数据布局定了还有一个容易忽略的细节内存对齐。AVX2 一次可以加载 32 字节AVX-512 是 64 字节。如果你用 malloc 动态分配一个 float 数组它默认只保证 16 字节对齐运气好落在 32 字节边界运气不好就落到一个尴尬的位置。当未对齐地址传给对齐 load 指令时轻则性能回退重则直接段错误。我吃过一次亏在调试版里一切正常开了 -O2 用了对齐 load intrinsic 后程序随机崩。后来用 addr2line 才发现是数组首地址没对齐。解决方式有三种posix_memalign或_aligned_malloc明确分配对齐内存。C17 的std::aligned_alloc注意它要求 size 是对齐值的整数倍容易踩坑。或者更简单声明一个 64 字节对齐的结构体数组内部放一个足够大的alignas(64) float data[]。对齐这件事同一个代码在不同编译器下的行为也会不一样尤其是 MSVC 的__declspec(align)和 GCC 的__attribute__((aligned))写法不同封装一层宏是值得的。2.4 Bandwidth vs Compute Intensity理解Roofline模型给优化划边界“Roofline”听起来很高深说白了就是画两条线一条叫内存带宽上限一条叫计算峰值上限。横轴是计算强度也就是每读入一个字节数据能做多少次浮点运算纵轴是实际能达到的性能。如果你的计算强度很低性能就会被带宽压着走再怎么优化计算指令都没用因为数据喂不进来。只有把计算强度推进到超过某个临界点性能才会跑到计算上限那一段。做数学库优化我习惯在动手前先算这个临界点。例如我的机器内存带宽实测约 30GB/sAVX2 单线程浮点峰值约 50 GFLOPs那么临界计算强度就是 50/30约等于 1.7 FLOP/Byte。如果某个运算的计算强度不足 1.7优化的重心应该放在缓存友好和数据复用上如果远超 1.7那就放心去压榨计算指令和 ILP。这解释了一个常见困惑为什么同样一台机器矩阵乘法能做得很高效而一个简单的逐元素加和却永远跑不上高 GFLOPs。逐元素加和每个元素只读一次写一次几乎没有复用计算强度极低上限就是内存带宽。这类运算再牛的库也只是把带宽用到极限而已。3. 让编译器把你的循环变成向量微架构级优化实操数据布局搞定之后下一步就是让你的循环尽量贴近 CPU 的向量执行单元。这一步是“榨干硬件”的关键也是最容易写出“看着优化了实际上更慢”代码的地方。3.1 自动矢量化编译器不是万能的很多编译器确实有自动向量化能力GCC 的 -O3、Clang 的 -O3 以及 Intel 编译器的 -xHost 都能把简单循环向量化。但这里有个前提循环结构必须足够简单、内存访问必须连续、循环次数最好能在编译期或运行期确定且不被复杂的别名分析拦下。我遇到最典型的自动向量化失败场景是这样的void vector_add(const float* a, const float* b, float* c, int n) { for (int i 0; i n; i) { c[i] a[i] b[i]; } }逻辑再简单不过了可如果调用方传入了c和a指向同一块内存的别名编译器为了安全就没法放心重排指令向量化就直接放弃。解决办法是加__restrict__告诉编译器这些指针没有重叠或者提前拷贝到临时缓冲区。GCC 的-fopt-info-vec能输出向量化报告建议读者开起来看一眼很多时候你会发现编译器“卡”在哪儿。3.2 依赖链与延迟为什么你的累加循环跑不满还有一个常被忽略的坑循环携带依赖。看这段代码float sum 0.0f; for (int i 0; i n; i) { sum a[i]; }每次sum a[i]都依赖上一次sum的值下一次加法必须等上一次的 FADD 指令完成。FADD 的延迟通常是 4 个时钟周期左右这意味着这个循环的吞吐被限制在每 4 个周期加一个元素和 CPU 实际上每周期能做 2 次浮点加法的能力相去甚远。解决办法是“多路累加”把累加器拆成 4 个独立的变量每次循环交替加float sum0 0, sum1 0, sum2 0, sum3 0; for (int i 0; i n; i 4) { sum0 a[i]; sum1 a[i1]; sum2 a[i2]; sum3 a[i3]; } float sum (sum0 sum1) (sum2 sum3);这样四条依赖链并行执行吞吐立刻翻 3-4 倍。更优雅的写法是让编译器在自动向量化时自己展开但显式写多路累加在大多数编译器上表现更稳定。拆成 4 路还是 8 路取决于目标微架构拥有几个加法执行单元一般 4 到 8 路比较合适再多也不会更快反而容易寄存器溢出。3.3 FMA、reduce与分支消除一个累加器的完整改造真正到库级别累加还会用到 FMAFused Multiply-Add也就是一条指令同时完成乘法和加法既省时间又减少中间舍入误差。对多项式求值、点积这类操作FMA 几乎是性能神器。看一个多项式求值的例子。普通写法是霍纳法逐项累乘加依赖链很长。优化后我可以把它拆成偶次项和奇次项两条独立链分别用 FMA 计算最后合并。这样依赖链深度减半吞吐自然上去。不过要注意FMA 改变了浮点运算的舍入顺序如果你的库对结果精度有逐位匹配的硬性要求这里需要做版本控制。另外能不用分支就不用分支。现代 CPU 的分支预测在简单规则下很准但面对数据相关的分支比如if (x 0) y 1; else y 0;预测失败代价大概是 15-20 个周期。这种场景可以直接用比较指令生成掩码再和常数做位运算做到无分支。数学库里的符号函数、绝对值、clamp 都可以这样优化。3.4 编译优化选项的组合拳单文件里写了很多优化代码但如果编译选项不对前面的辛苦全部白费。我常用的组合是GCC/Clang-O3 -marchnative -ffast-math -fno-math-errno -fomit-frame-pointer -funroll-loops用-ffast-math之前要想清楚它会假设 NaNs/Infs 不会出现、允许重排浮点运算、关闭 errno 设置对安全性要求高的代码慎用。我自己一般会单独为内部性能敏感模块开对外 API 保持保守。大项目里开-flto链接时优化收益也很大它能让跨编译单元的常量传播和内联共享变成可能但会显著增加编译时间和内存占用。我在实际项目中还发现一个细节-marchnative会启用编译器对当前 CPU 的自动检测如果你把编译好的库放到另一台老机器上可能会出现 illegal instruction 崩溃。分发库的时候应该用-marchx86-64-v3这类更保守的指令集级别而不是图省事直接 native。4. 多线程并行线程数不等于性能单核优化做到位之后接下来必然是并行。但并行不是“开几个线程然后等奇迹”这一层踩坑最多也最容易做出负优化。4.1 两种并行切分策略数据分片 vs 管道分片数学库的并行通常用数据分片因为每个数据元素的计算是独立的天然适合把矩阵按行块或列块切给不同线程。管道分片在数学库里用得少因为每一步的输出要被下一步消费线程间经常需要同步同步成本很容易吃掉并行收益。举个例子矩阵乘法把输出矩阵的 i 维按行切成 4 段每段交给一个线程去算对应的行块这是一种典型的数据分片。注意切分要尽量让每段运算量相等避免一个线程算得飞快另一个还在跑长尾。如果矩阵规模不确定最简单的办法是动态调度比如用 OpenMP 的 schedule(dynamic, chunk_size)。4.2 伪共享你不敢信的“负优化”案例伪共享是我见过最隐蔽的并行性能杀手。它的原理是CPU 的缓存一致性是以 cache line通常 64 字节为单位的。如果两个线程分别修改不同变量但这两个变量凑巧落在同一个 cache line 里缓存一致性协议会让两个核心不断互相通知“这个缓存行被改了”导致明明互不相干的写入变成串行化性能暴跌。我做过一个实验开 4 个线程每个线程对数组的某一段独立求和然后把结果写进四个不同的全局变量。仅仅因为这四个 float 变量紧挨着存放在一起性能比把它们隔开 64 字节对齐慢了三倍多。解决办法是给不同线程的结果变量加 padding让它们落在不同 cache linestruct alignas(64) PerThreadResult { float value; char padding[60]; };这样就把每个线程的写入隔离到独立缓存行。类似的还有循环切块大小最好也按 cache line 的整数倍去切避免两个线程把同一行数据拆开处理。4.3 时间分片与调度开销OpenMP schedule的选择用 OpenMP 时schedule 的选择直接影响负载均衡。static调度是编译期均分开销最小但如果数据分布导致不同区域计算量不等线程会空等。dynamic调度是运行时按块领取能自动均衡但每领一块都有锁开销。guided是介于两者之间块大小递减适合循环内部计算量差异大的情况。我自己的经验是如果循环内计算量均匀优先static如果计算量难以预估选guided而非dynamic因为它的调度次数更少。有一回矩阵维度正好让 8 线程的 static 切分把一个线程分到几乎为空的行区间另一个线程却要算整个最重的列块性能直接掉了一半。换成 guided 之后才恢复正常。4.4 同步原语与原子操作在数学库里合理边界数学库内部的同步越少越好。多线程累加场景里最常见的错误是直接把sum local_value放进临界区或者频繁使用原子加。这类操作在低竞争下勉强能用高竞争时锁和原子操作的开销甚至比计算本身还高。正确的做法是每线程维护局部结果最后用一个 O(n_threads) 的归约合并结果。这个归约阶段即使加锁竞争也只发生一次开销忽略不计。如果并发级别不高也可以直接用无锁的原子加做最终归约但要注意原子操作对 cache line 的影响和上面伪共享是同一个问题。5. GEMM三层优化实战从朴素实现到接近硬件峰值泛泛谈优化容易陷入空洞下面我用最经典的 GEMM通用矩阵乘法C A * B完整走一遍优化流程。GEMM 是整个 BLAS 体系的基石也是几乎所有线性代数库评测里的重头戏。5.1 版本0朴素三重循环性能基准void matmul_naive(const float* a, const float* b, float* c, int n) { for (int i 0; i n; i) { for (int j 0; j n; j) { float acc 0.0f; for (int k 0; k n; k) { acc a[i * n k] * b[k * n j]; } c[i * n j] acc; } } }这个版本的问题前面已经讲过b 是按列访问的内层循环每次b[k][j]都在跳跃地址cache miss 率极高。在 n1024 时实测单线程耗时 310ms 左右。作为基准它的价值是告诉我们在没有任何优化时的代价。5.2 版本1转置缓存分块向量化第一步先把 b 转置一遍让内层循环里两个矩阵都变成按行访问void matmul_t1(const float* a, const float* b, float* c, int n) { float* bt new float[n * n]; for (int i 0; i n; i) for (int j 0; j n; j) bt[j * n i] b[i * n j]; for (int i 0; i n; i) { for (int j 0; j n; j) { float acc 0.0f; for (int k 0; k n; k) { acc a[i * n k] * bt[j * n k]; } c[i * n j] acc; } } delete[] bt; }这一步已经能降到 76ms 左右。但还不够因为每次内层循环取 a 的一行和 bt 的一行这两行都来自主存没有充分利用 cache。下一步是分块让计算的最小单元在一个块内反复复用 L1/L2 缓存void matmul_t2(const float* a, const float* b, float* c, int n, int block) { for (int i0 0; i0 n; i0 block) { for (int j0 0; j0 n; j0 block) { for (int k0 0; k0 n; k0 block) { for (int i i0; i i0 block; i) { for (int k k0; k k0 block; k) { float a_tmp a[i * n k]; for (int j j0; j j0 block; j) { c[i * n j] a_tmp * b[k * n j]; } } } } } } }block 大小一般选 32 或 64。选小了块太小每次加载的 a、b 片段复用次数不够选大了块内数据超过 L2 容量又会频繁换出。这个版本在 32 块大小时耗时能再降到 30ms 上下。最后叠加 AVX2 向量化内层 j 循环一次算 8 个 float耗时进一步降到 15ms 左右。5.3 版本2多线程进一步调优把分块后的 i0 层直接抛给 OpenMP#pragma omp parallel for schedule(guided) for (int i0 0; i0 n; i0 block) { ... }注意我选择对 i0 并行而不是内层 j0因为并行粒度越粗线程切分和合并的开销就越小。在 4 线程实测下这个版本能跑到 9ms 左右。进一步调优还有两个方向一是把 block 的尺寸按矩阵维度动态调整二是让内层循环用restrict修饰指针、避免别名误判。还有一个容易忽视的点c[i * n j]在内层循环里被反复读写如果 block 不太大c 的对应片段会留在 L2 里但为了保险可以为每个输出块开一个 block×block 的临时累加数组算完再写回。这样 c 的读改写从 main memory 变成了 cache 内的局部操作效率更高。5.4 性能数据对比与提升路径检验把这几个版本放一起就能看到每一层优化的贡献版本优化手段耗时N1024, 单线程提升倍数v0朴素三重循环310ms1xv1转置 b76ms4.1xv2转置 分块(block32)30ms10.3xv3v2 AVX2 向量化15ms20.7xv4v3 4线程9ms34.4x从数据看转置和分块带来的收益甚至超过向量化。这验证了之前说的“访存优先”原则先解决数据怎么流动再考虑计算指令怎么排。向量化的收益是在缓存优化做完之后才显出来如果一开始就上 AVX 而没有分块效果会大打折扣。6. 验证与基准测试跑得快的前提是算得对优化做得越多正确性验证就越不能省。尤其是做数学库服务对象可能是别的工程师在使用一个错误的结果比慢十倍还要糟糕。6.1 测试数据设计三组必测用例我一般会设计三组测试第一组是“已知答案”用例。比如让矩阵元素全是 1或者单位矩阵乘任意矩阵结果必定一眼能看出对错。这类用例跑一遍基本功能错误立刻暴露。第二组是“随机数据 高精度参考”。用 double 甚至 long double 实现一个相同运算的朴素版本作为参考然后对 float 的优化版本做相对误差对比。随机数据能覆盖到各种数值区间避免某个特定输入下精度刚好够好。第三组是“边界条件”。包括 n0、n1、n非 4 的倍数、NaN/Inf 输入、非常大和非常小的数值。很多优化代码在 n 不是 SIMD 宽度倍数时会写出错误的后处理逻辑这组用例就是抓这个的。6.2 ULP误差控制浮点数学库的“误差牌照”数学库和普通业务代码不同它不能只说“结果大概差不多”而是要以 ULPUnit in the Last Place为单位定义误差上限。ULP 是指两个相邻浮点数之间的间隔1 ULP 误差就是你结果的最后一位可能有偏差。对不同函数业界对不同函数有不同要求加减乘除这类四则运算IEEE 754 要求正确舍入误差不超过 0.5 ULPexp、log、sin、cos 这类超越函数通常允许 1-4 ULP。如果你的库要跟 BLAS/LAPACK 对标误差也得定在相近水平。所以我在优化每个函数时都会先写一个“慢但正确”的参考实现再用随机输入批量对比要求误差在预设 ULP 范围内。这里有个经验浮点优化里常见的错误不是算错而是重排运算顺序后误差被放大。比如前面提到的多路累加和 FMA改变累加顺序会改变舍入误差虽然通常在一个 ULP 以内但如果库的接口契约里写死了“逐位匹配”就需要专门跑全量黄金样本。6.3 基准测试方法论高频抖动、预热与统计口径基准测试也不能随便跑一遍就看数字。现代 CPU 的频率是动态的同一个循环冷启动和热启动能差 30% 以上。我做基准的标准流程是用taskset或numactl把进程绑到固定核避免操作系统调度导致核心漂移。先在循环外跑一轮预热把缓存填起来触发频率拉升。正式测量时每个尺寸跑 5-10 次取中位数而不是平均值。中位数能有效剔除系统中断和偶发抖动的影响。同时记录功耗和温度防止机器过热降频影响结论。另外我刚入行时犯过一个错优化完一个函数后随手用系统自带的 chrono 计时结果把打印、内存分配也一起算进去了导致“优化效果”被其他开销稀释。正确的做法是只对被测试的核心计算循环计时内存分配、初始化、打印都放到计时区间之外。7. 真实项目的额外复杂度不是只有矩阵乘法GEMM 是最典型的例子但一个真正可用的数学库还要面对大量“不那么规整”的问题。这些才是拉开库之间差距的地方。7.1 特殊函数、向量运算与数值稳定性比如指数、对数、三角函数这些特殊函数。它们的优化不完全是“算得快”而是“在快的同时保持数值稳定”。比方说计算log(1x)如果直接实现log(1.0 x)当 x 很小时1.0 x 会先被舍入到 1.0导致结果直接变成 0。正确做法是用专门的log1p展开式对很小和很大的 x 分别处理。同样地exp(x)在 x 接近边界时要防止溢出到 Inftan(x)在接近π/2时会爆炸这些边界行为必须仔细定义。数学库的每一个函数背后都有一本“数值食谱”有的来自经典参考书有的来自平台原厂库的对照。7.2 稀疏矩阵与数据结构取舍真实场景里很多矩阵是稀疏的比如图计算、有限元模拟中的矩阵非零元素占比可能只有 1%。这时如果你还用稠密矩阵的存储方式去优化即使 GEMM 再快也是巨大的浪费。稀疏矩阵的优化重点变成“只访问非零元素”。常用的存储格式有 CSRCompressed Sparse Row和 CSCCompressed Sparse Column。CSR 用一个数组存所有非零值再用两个数组记录每一行的起始位置和每行非零元素的列号。这种方式省内存但访问模式变得不规则很难像稠密矩阵那样用 SIMD 吃到收益。我做过一个稀疏矩阵向量乘普通 CSR 实现遇到 GPU 上的分支发散严重换了 ELLPACK 格式把最多非零元素的那行定成宽度其他行补零之后SIMD 友好度大幅提升速度翻了 3 倍。代价是内存略多一点非零分布极不均匀时会浪费严重。7.3 精度与性能的矛盾什么时候该用float什么时候必须double很多人以为 float 一定比 double 快这是误解。在支持双精度 FMA 的芯片上double 的吞吐和 float 差不多比如某些服务器 CPU 的 FP64 单元是完整宽度的只是内存带宽消耗减半。而有些消费级 CPU 的 FP64 单元是被阉割的跑 double 就明显慢。所以“用 float 还是 double”取决于目标硬件和数值需求。我的经验是如果业务数值范围很大、需要高精度累加先用 double 做原型确认误差水平后再决定是否下沉到 float。如果确实需要 float 性能但担心累加误差可以考虑用 Kahan 补偿求和——它只增加少量指令却能把累加误差降到一个很低的水平。8. 我在多次重写数学库后最想分享的几点经验最后说点个人感受可能比前面的技术细节更有价值。第一优化是迭代式的不是一蹴而就。不要试图“一次写出最优代码”而是先写好正确版本再逐步按访存、指令级并行、多线程的顺序一层层加优化。每做一步就测一次既能定位瓶颈也能防止后面的改动破坏前面的成果。第二用性能剖析器验证你的直觉不要靠猜。我经常打开perf stat看 cache-misses、branch-misses、IPC 这些指标。有一次我以为瓶颈是分支预测失败数据出来却是 L1 cache miss 占了大头方向完全反了。第三版本管理里记录每个 commit 的性能增量。数学库的重构风险极高没有性能回归记录的优化很容易在几个月后被一个看似无害的改动悄悄毁掉。我会在 CI 里挂一个简单的性能基准作业每次 merge 前自动跑一遍关键 case超过阈值就报警。第四不要迷信“某个库就是最快”的传说。每个库都有自己的适用场景。MKL 在 Intel 平台上是标杆但换到 AMD 或者 ARM 上未必最优。真正可靠的做法是拿你自己的业务矩阵尺寸和访问模式去实测跑出一个可复现的对比数据再下结论。数学库这条路做久了你会发现真正难的不是数学公式而是让公式在真实的硬件上跑得又快又稳。每一次优化都是和硬件的一次对话。希望这篇文章能帮你把第一句“对话”说清楚。