Jack Dongarra 的图灵奖与并行数值线性代数库之路
一点从业人员的吐槽。
1. Dongarra 的图灵奖，和它的反面
前几天，Jack Dongarra 教授荣获 2021 年图灵奖。当天，笔者的微信朋友圈里充满了各种和 Dongarra 教授的合影，笔者身边几乎所有从事（过）与高性能计算相关工作的人都为这个消息感到高兴。当然，也有反面的声音，比如 这个知乎回答。在这个知乎问题下面，这个回答算是比较难得地提到了一些值得思考的问题，其他的回答除了一些同行的肯定，大部分都没什么价值（师老党立也来凑热闹了，他自称是做计算的，但是他的论文在哪里呢，笔者暂且蒙在古里）。这个唱反调的回答主要强调了两点：
- 现在分布式稠密线性代数的并行算法的主要组成部分是 Communication-avoiding 算法，主要贡献人是James Demmel 教授等。Jack Dongarra 现在最新的 SLATE 这种分布式稠密线性代数库的算法就是 CA 那套。
- Top500 的 LINPACK 更加是一个没什么太大用处的排名。实际问题大多都是以稀疏矩阵为基础的。
这两个论点都很有意义。不幸的是，具体的回答内容显示出答主似乎对自己的两个论点没有很好的理解。下面我们深入探讨一下这两个论点。
2. 新一代并行数值线性代数库 SLATE
经过了几十年的发展，稠密数值线性代数库在共享内存机器上有了两套标准接口 BLAS 和 LAPACK，并且有 netlib 的未调优参考实现。而在分布式内存机器上，PBLAS (Parallel BLAS) 和 ScaLAPACK (Scalable LAPACK) 虽然不是名义上的标准接口，但是事实上的标准接口。由于很少有人单独使用 PBLAS, 而 ScaLAPACK 自某个版本以后就自带了 PBLAS，故我们一般不加区分地只称 ScaLAPACK。ScaLAPACK 的参考实现是用 Fortran 写的，已经有超过 20 年的开发历史了。一方面，ScaLAPACK 代码的可读性已经大不如前；另一方面，ScaLAPACK 设计之初仅仅面向 MPI 并行，无法很好驾驭现代的多核和 GPU 计算节点。
自从 Exascale Computing 被摆上桌面，Dongarra 领导的团队就开始为之准备和开发新一代的并行数值线性代数库 SLATE, Software for Linear Algebra Targeting Exascale. SLATE 采用了现代的 C++ 和 OpenMP 来实现各种线性代数算法（因此代码量比起 ScaLAPACK 要小了一个量级），并且广泛吸收了近 20 年来数值线性代数算法的新成果。反面回答认为 SLATE 现在采用的算法都是 Demmel 的那一套 CA 算法，因此 Dongarra 在并行数值代数算法上贡献并没有那么大。我们姑且不讨论前半句这种文无第一的问题，前半句本身也不太对。让我们来看看 SLATE 自己怎么说。
数值线性代数库主要包含几类算法：矩阵乘法、LU 与 Cholesky 分解、QR 分解、特征值与奇异值分解。我们逐个讨论。
首先是矩阵乘法。笔者在 前一篇文章 中分析和介绍了主要的几种矩阵乘法算法，其中包括了 Demmel 的 Communication-Optimal 2.5D 算法以及摸到了通信量理论下界 CARMA 算法。然而 SLATE 采用的是 SUMMA 算法 [1]，和 ScaLAPACK 一样。对此，SLATE 的解释是：
The current emphasis in SLATE is on simplicity and robustness. More advanced algorithms are up for consideration in the future.
然后是 LU 和 Cholesky 分解。这两个分解的计算模式高度相似，所以我们只讨论 LU 分解。为了偷懒，我从 [2] 中截取了一段内容来说明 LU 分解在做什么：
SLATE 在做 LU 分解的时候，对 panel factorization 采用了 Communication-Avoding LU 算法 [3]。在 panel factorization 完成以后，SLATE 调用自己的矩阵乘法算法来计算 rank-k update 以更新 trailing matrix。实际上，这个算法并不是通信最小化的算法。通信最小化的算法仍然是 Demmel 提出的 2.5D LU [4]。与 2D LU 相比，2.5D LU 其实非常好理解：用 2.5D 矩阵乘法来计算 trailing matrix rank-k update。
接着是 QR 分解。关于 QR 分解，Demmel 有两个 CA 算法：TSQR 和 CAQR [5], SLATE 采用了 CAQR [6]。TSQR 针对 tall-and-skinny matrix, 假设这样的矩阵只进行行划分，每个进程持有若干行，通过二叉树归并的方式计算出矩阵的 Householder QR 分解。CAQR 是为通用矩阵设计的，其主体框架是分块化的 Householder QR。在算 panel factorization 的时候，CAQR 采用 TSQR 来算出当前 panel 的 householder vectors，随后对 trailing matrix 的更新和 ScaLAPACK 中的做法完全一致。[5] 中的 Table 16 也指出，Parallel CAQR 和 ScaLAPACK 中的 Householder QR 有着同阶的通信复杂度，前者比后者好的地方是复杂度中的某些部分少了一些常系数。
最后是特征值和奇异值分解。这两种分解的算法比较复杂，当前并没有什么 communication-avoiding 的算法，SLATE 的实现也是使用现有的 1-stage 和 2-stage 方法 [7]。
总结一下，SLATE 采用了很多 CA 算法，但并没有全部采用 CA 算法；同时，其采用的 CA 算法不全部是最小化通信的 CA 算法。在渐进通信复杂度上，SLATE 和 ScaLAPACK 并没有区别。总体而言，SLATE 依靠常系数优化和对 MPI + OpenMP 优化来提升了并行数值线性代数操作的可扩展性。
3. 『历史的终结及最后之人』
上一节对 SLATE 里面矩阵乘法和 LU 分解算法的分析很容易导向一个问题：为什么不能采用达到理论通信复杂度下界的 3D 矩阵乘法算法？这样的话起码矩阵乘法和 LU 分解都可以达到最优的渐进通信复杂度了。
这个问题的答案其实笔者在 前一篇文章 中已经指出来了：CARMA 和 COSMA 这些达到理论通信复杂度下界的矩阵乘法算法，都对矩阵分布有特殊的要求。笔者认为，ScaLAPACK 设计中的一个要点是它采用了 2D block-cyclic (2DBC) 的矩阵分布 [8]，并且所有的算法都是基于这个矩阵分布的。笔者非常讨厌这个反直觉的矩阵分布方式，但是在经历过相关的工作以后，笔者认为这是最不坏的方式。之所以这么说，是因为线性代数的主要算法是各种矩阵分解，随着这些分解算法的进行，trailing matrix 都在变小，大部分是行和列都在变小，有一些是行或者列变小。因此，所有静态的、连续的 1D 或者 2D 矩阵划分，都将随着算法的运行出现有些进程不再持有任何 trailing matrix 的情况。此时如果不重新分布 trailing matrix，则会负载不均浪费计算资源；如果重新分布 trailing matrix，则通信模式变化巨大，且通信量也巨大。2DBC 的使用保证了在矩阵分解的任意阶段，各个进程持有的 trailing matrix 大小几乎相同，差值是一个较小的常数，因此计算和通信的负载大体上是均衡的。如果矩阵乘法使用 COSMA，则必然要在矩阵乘法和其他计算之间来回转换矩阵分布，这个代价同样是不可接受的。
由此看来，并行稠密线性代数算法已经摸到了它的天花板。矩阵乘法、LU 分解、QR 分解都已经有达到理论通信复杂度下界的算法，但这些算法又难以统一在相同的矩阵分布和并行框架下面。SLATE 已经充分利用了多层次的硬件并行和数据访问模式来提高性能，10年以后如果人类的超级计算机还是像现在一样是十年前的放大版，那么接替 SLATE 的新软件还能如何提高性能？
让我们回到反对答案的第二点，即稠密矩阵线性代数已经不实用了，HPL 也只是个用处不大的排名。诚然，现在大规模的科学计算与仿真，其算法主体基本都不再使用稠密线性代数了，唯一的例外是深度神经网络。但是有些基于稀疏矩阵的算法依然需要稠密线性代数运算 [9]。因此，稠密矩阵线性代数并行算法的开发仍然有其必要性。至于 HPL 的价值，笔者非常同意上海交通大学廖秋承研究员的观点，大意为 HPL 仍然是最好的超级计算机稳定性测试工具，其性能波动、能耗波动仍然具有许多值得研究的地方。至于说 Top500 只是给了各个超算中心虚名去做广告要经费，那么哪个 benchmark 不是这样呢？HPL 过度强调浮点性能，HPCG 过度强调内存带宽（甚至可以被模型很好地预测出性能，[10]），Graph500 则毫无反应浮点计算能力。只要怪物亮血条不锁血，总有刮死它的办法。只要有 benchmark 允许改代码，总有针对性优化来刷榜的方法。求全责备，大可不必。
致谢
笔者衷心感谢 SYSU-SCC 水群群友的全天候高强度水群行为，群友的正论激发了本文的写作。
参考文献
- Parallel BLAS Performance Report, link
- 面向国产异构处理器的 HPL-AI 基准实现及优化, link
- Communication Avoiding LU with Tournament Pivoting in SLATE, link
- Communication-Optimal Parallel 2.5D Matrix Multiplication and LU Factorization Algorithms, link
- Communication-Optimal Parallel and Sequential QR and LU Factorizations, link
- Least Squares Performance Report, link
- Implementing Singular Value and Symmetric/Hermitian Eigenvalue Solvers, link
- The Two-dimensional Block-Cyclic Distribution link
- Numerical Methods for Large Eigenvalue Problems, link
- Performance Modeling of the HPCG Benchmark, link