PAYATHON 2026

为什么 ScaLAPACK 矩阵乘法从单节点扩展到双节点后反而变慢？

数据阿明

结论

现有信息不足以确定性能下降的唯一原因。代码中需要先检查两个问题：大矩阵的内存分配和下标计算可能发生 int 溢出；输出公式使用的单位是 GFLOP/s。

两个 MPI 进程会组成 1 × 2 的进程网格，每个进程仍承担节点内全部核心的线程计算。性能排查应关注实际 BLAS 线程数、CPU 绑定、NUMA 分布和 MPI 通信时间。

“互连带宽为 12.5 GB/s”不足以证明通信只需要几秒。还需要测量实际使用的通信路径、传输量，以及通信是否阻塞计算。

先确认单节点基准是否有效

内存分配与下标计算可能发生整数溢出

原代码使用：

double *A = malloc(m_A * k_A * sizeof(double));

m_A 和 k_A 都是 int，所以 m_A * k_A 会先进行 int 乘法，再与 sizeof(double) 相乘。后面的 sizeof(double) 无法避免前一步溢出。

如果实际尺寸为 50000 × 50000，单进程需要保存的元素数是：

50000 × 50000 = 2,500,000,000

这个数超过了常见的 32 位有符号 int 上限 2,147,483,647。有符号整数溢出属于未定义行为，初始化中的 i + j * m_A 也存在同样的风险。

原问题只说尺寸“约为 50K”，是否触发溢出仍取决于实际参数、整数宽度和本地矩阵尺寸。如果确实使用上述尺寸和 32 位 int，应先修复这些表达式，再重新建立性能基准。

下面是 A 的分配与初始化示例，适用于尺寸已验证为正、且 size_t 足以表示所需容量的环境。B、C 也要作同样处理：

#include <stdint.h>

size_t rows_A = (size_t)m_A;
size_t cols_A = (size_t)k_A;

if (cols_A != 0 && rows_A > SIZE_MAX / cols_A) {
    MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
}

size_t count_A = rows_A * cols_A;
if (count_A > SIZE_MAX / sizeof(double)) {
    MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
}

double *A = malloc(count_A * sizeof(*A));
if (A == NULL && count_A != 0) {
    MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
}

#pragma omp parallel for
for (int j = 0; j < k_A; j++) {
    for (int i = 0; i < m_A; i++) {
        A[(size_t)i + (size_t)j * rows_A] = 1.0;
    }
}

这里修复的是应用代码中的容量和下标计算。所链接的 ScaLAPACK、BLAS 是否支持该问题规模，以及 C 声明与库的整数 ABI 是否匹配，还需要结合实际构建配置确认。

检查内存容量、描述符和计算结果

当 M = K = N = 50000 时，三个双精度矩阵的数据合计约为 60 GB，尚未计入库内部缓冲区和其他内存开销。需要确认单节点有足够的可用内存，运行期间没有交换或明显的内存压力。

每次调用 descinit_ 后都应立即检查 info。原代码连续覆盖同一个 info，也没有检查分配结果，初始化失败后可能继续运行。

当前输入全部为 1.0，且 alpha = 1.0、beta = 0.0，因此正确结果中的每个元素都应为 K。可以在计时区间之外检查各进程持有的 C，并通过 MPI 汇总最大误差。确认结果正确后再比较性能。

明确吞吐量单位

原公式为：

2.0 * 1e-9 * M * N * K / elapsed

输出单位是 GFLOP/s。若要直接输出 TFLOP/s，应使用 1e-12：

double tflops = 2.0 * (double)M * (double)N * (double)K
              / elapsed / 1e12;

如果报告的 3.8 和 1.5 已经由原输出换算得到，就没有单位问题；如果直接读取程序输出，则需要重新核对。

为什么两个节点值得单独排查

MPI 进程网格发生了变化

这里每节点只有一个 MPI task，因此节点数等于 MPI 进程数。根据代码中的网格选择逻辑：

MPI 进程数BLACS 网格
11 × 1
21 × 2
31 × 3
42 × 2

从一个进程增加到两个进程后，计算开始需要跨进程交换数据；增加到四个进程时，网格又变成了二维布局。各配置的通信模式和本地计算形状不同，扩展曲线因此不一定平滑。

1 × 2 网格本身不能解释性能为什么恰好从 3.8 降到 1.5 TFLOP/s，但它值得对照测试。可以在双进程时比较 1 × 2 和 2 × 1；增加每节点的 MPI 进程数后，再比较接近方形的网格。

峰值带宽无法直接换算成 pdgemm 的通信耗时

将一个矩阵的大小除以链路带宽，只能得到一次理想化传输的时间。实际 pdgemm_ 涉及分块数据交换、集合通信，以及同步和本地计算之间的依赖，不能直接等同于一次完整矩阵传输。

需要确认：

  • MPI 实际使用了哪种传输通道，是否走预期的高速互连。
  • 消息规模接近实际情况时，跨节点通信能达到多大带宽。
  • pdgemm_ 执行期间，多少时间用于 MPI，多少时间用于本地 BLAS。
  • 是否有进程计算较慢，导致其他进程等待。

原程序末尾的 MPI_Barrier 会把进程间的等待计入总时间。这适合衡量整次分布式运算的完成时间，但无法直接区分通信开销和负载不均衡。

检查线程数、绑定和 NUMA

OMP_PLACES=cores 指定 OpenMP 线程可使用的位置，但不能单独确定实际线程数。Slurm 为 task 分配 128 个 CPU，也不能证明所链接的 BLAS 确实启动了 128 个计算线程。

先记录实际使用的 BLAS 实现、线程后端和线程配置。应用初始化循环的 OpenMP 线程与 BLAS 计算线程可能采用不同设置。

对于使用 OpenMP 线程后端的配置，可以先显式设置并验证：

export OMP_NUM_THREADS=128
export OMP_PLACES=cores
export OMP_PROC_BIND=spread

这是一组对照配置，是否合适需要实测。如果库使用自己的线程控制变量，应按该库的配置方式设置；例如实际链接 MKL 时，需要同时核对 MKL_NUM_THREADS。

还要检查 srun --cpu-bind=cores 给每个 task 分配的实际 CPU 集合。一个 task 使用 128 个核心时，如果节点包含多个 NUMA 域，线程与内存页的位置也可能影响性能。原代码并行初始化数组，有助于按首次写入位置分配内存页，但仍需确认初始化线程与计算线程的绑定是否一致。

可比较两种布局：

  • 每节点 1 个 MPI 进程，使用节点内全部核心。
  • 每节点多个 MPI 进程，每个进程使用更少线程，并结合实际 NUMA 拓扑绑定。

比较时，保持每节点使用的物理核心总数一致，并检查是否存在多个进程的线程争用同一组核心。增加 MPI 进程数可能改善局部性，也可能增加通信，哪一种更快需要实测。

建议的定位顺序

  1. 修复并验证程序：检查整数溢出、内存分配、描述符返回值和结果误差，记录准确的矩阵尺寸与性能单位。
  2. 固定运行环境：记录编译选项、MPI、ScaLAPACK、BLAS、线程数和 CPU 绑定，确保单节点与双节点使用相同配置。
  3. 测量本地计算：在相同线程和绑定条件下测试本地 DGEMM，尽量覆盖分布式计算中可能出现的矩阵形状，而不只测一个大型方阵。
  4. 测量跨节点通信：验证 MPI 通道和实际带宽，再用 MPI 性能分析工具观察 pdgemm_ 中的通信与等待时间。
  5. 比较进程布局与网格：测试每节点不同 MPI 进程数、对应线程数，以及可用的 BLACS 网格。
  6. 最后调整 block size：已尝试多个 block size 且现象一致，继续只改这个参数，暂时难以缩小原因范围。

各配置应进行预热和多次测量，报告中位数与波动范围。若本地 BLAS 性能稳定、跨节点后 MPI 等待明显增加，应继续检查通信路径与网格；若本地计算本身就慢，则优先排查线程配置、绑定和 NUMA。现有代码尚不足以把性能下降归因于某一个因素。

备注：内容仅供参考。