☰
矩阵乘法 MPI 并行实战:从按行分块到性能优化
2026/10/15 5:09:02 网站建设 项目流程

简介:这份资源面向并行计算初学者与科学计算方向的学习者,聚焦基于MPI的矩阵乘法并行实现,帮助读者理解分布式内存架构下如何将矩阵运算拆分到多处理器协同完成。压缩包内共1个docx文件,约15KB,以文档形式整理算法原理、通信函数用法与代码实现思路,便于边看边对照练习。内容围绕MPI_Init、MPI_Comm_rank、MPI_Comm_size、MPI_Bcast、MPI_Scatterv、MPI_Gatherv、MPI_Finalize等核心接口展开,讲解矩阵A分发、矩阵B广播、各进程局部计算与结果聚合的完整流程,并涉及BLOCK_SIZE等优化思路。目前已有443人学习,适合需要掌握MPI矩阵乘法并行计算、理解分布式内存架构与通信函数配合方式的读者参考,也可作为课程实验与项目实践的辅助材料。

1. 矩阵乘法 MPI:从单机瓶颈到多进程并行的真实分水岭

很多人第一次接触 MPI,都是被矩阵乘法带进来的。原因很直接:一个 2048×2048 的方阵乘法,单线程朴素实现要跑几十秒甚至几分钟,而一旦切到 MPI 多进程,同样的机器上时间能掉到原来的几分之一。这个落差足够让人产生“我是不是摸到了高性能计算的门”的错觉。但矩阵乘法 MPI 真正要解决的问题,不是“能不能并行”,而是“怎么切数据、怎么传数据、怎么让通信不把计算吃干净”。它适合两类人:一类是刚学完 MPI 点对点通信、想找一个完整例子把 Send/Recv、广播、归约串起来的学生;另一类是要把已有串行数值代码搬到集群上、但不确定分块策略和通信开销怎么权衡的工程师。这篇文章就按“先能跑、再跑对、最后跑快”的顺序,把矩阵乘法 MPI 的落地路径拆开讲。

2. 矩阵乘法 MPI 的分块逻辑:为什么按行切是最稳的起点

2.1 并行矩阵乘法的三种切法,先选对再谈优化

矩阵乘法 C = A × B,最朴素的串行写法是三层循环。并行化的本质是把 C 的计算任务分给多个进程,每个进程只算一部分结果。常见切法有三种:

第一种是按行切 A。每个进程拿 A 的若干行,和完整的 B 相乘,算出 C 的对应行。这种切法通信量最小,因为 B 只需要广播一次,A 的各行本来就分散在各进程本地。缺点是每个进程都要存完整的 B,内存占用不随进程数下降。

第二种是按列切 B。每个进程拿 B 的若干列,需要完整的 A。通信模式和按行切对称,但 A 通常比 B 更需要复用,所以实际中用得少一些。

第三种是二维分块,把 A、B、C 都切成棋盘格,每个进程负责一个子块。这种切法内存扩展性最好,但通信逻辑复杂,需要处理行广播和列广播两轮通信,新手容易在进程网格映射上翻车。

我一般建议从按行切 A 开始。原因不是它最快,而是它的通信模式最清晰:一次广播 B,一次收集 C,中间全是本地计算。你先把这条链路跑通,再去看二维分块,才知道复杂度花在哪里。

2.2 按行切分的数学表达与进程映射

假设 A 是 M×K,B 是 K×N,C 是 M×N。用 p 个进程,把 A 按行分成 p 块,每块大约 M/p 行。进程 rank 负责 A 的第 rank 块,记为 A_local,形状是 (M/p)×K。B 是 K×N,所有进程都需要完整的一份。

计算过程:

  • 进程 0 读入 A 和 B。
  • 把 B 广播给所有进程。
  • 把 A 按行散射给各进程,或者各进程自己从文件读自己那几行。
  • 每个进程算 C_local = A_local × B。
  • 把所有 C_local 收集到进程 0,拼成完整的 C。

这里有一个容易忽略的点:如果 M 不能被 p 整除,分块要处理余数。常见做法是前 r 个进程多拿一行,r = M mod p。这个细节不处理,小矩阵测试时看不出来,一上大矩阵就会越界或者丢行。

2.3 最小可运行版本:广播 B + 本地计算 + 收集 C

下面是一个能直接跑的 MPI 矩阵乘法最小版本,用 C 写,按行切 A,B 用 MPI_Bcast 广播,C 用 MPI_Gatherv 收集。为了聚焦通信逻辑,矩阵用一维数组按行优先存储。

#include <mpi.h> #include <stdio.h> #include <stdlib.h> // 生成矩阵,按行优先存成一维 void init_matrix(double *mat, int rows, int cols, int seed) { srand(seed); for (int i = 0; i < rows * cols; i++) { mat[i] = (double)(rand() % 10) / 10.0; } } int main(int argc, char **argv) { MPI_Init(&argc, &argv); int rank, size; MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &size); const int M = 1024; // A 的行数 const int K = 1024; // A 的列数 / B 的行数 const int N = 1024; // B 的列数 double *A = NULL; double *B = (double *)malloc(K * N * sizeof(double)); double *C = NULL; // 进程 0 准备完整数据 if (rank == 0) { A = (double *)malloc(M * K * sizeof(double)); C = (double *)malloc(M * N * sizeof(double)); init_matrix(A, M, K, 1); init_matrix(B, K, N, 2); } // 广播 B 给所有进程 MPI_Bcast(B, K * N, MPI_DOUBLE, 0, MPI_COMM_WORLD); // 计算每个进程分到多少行 int base = M / size; int rem = M % size; int local_rows = base + (rank < rem ? 1 : 0); // 计算散射时的偏移和计数 int *sendcounts = NULL; int *displs = NULL; if (rank == 0) { sendcounts = (int *)malloc(size * sizeof(int)); displs = (int *)malloc(size * sizeof(int)); int offset = 0; for (int i = 0; i < size; i++) { int rows = base + (i < rem ? 1 : 0); sendcounts[i] = rows * K; displs[i] = offset; offset += rows * K; } } double *A_local = (double *)malloc(local_rows * K * sizeof(double)); // 把 A 按行散射到各进程 MPI_Scatterv(A, sendcounts, displs, MPI_DOUBLE, A_local, local_rows * K, MPI_DOUBLE, 0, MPI_COMM_WORLD); // 本地矩阵乘法 double *C_local = (double *)calloc(local_rows * N, sizeof(double)); for (int i = 0; i < local_rows; i++) { for (int k = 0; k < K; k++) { double a = A_local[i * K + k]; for (int j = 0; j < N; j++) { C_local[i * N + j] += a * B[k * N + j]; } } } // 收集 C 到进程 0 int *recvcounts = NULL; int *recvdispls = NULL; if (rank == 0) { recvcounts = (int *)malloc(size * sizeof(int)); recvdispls = (int *)malloc(size * sizeof(int)); int offset = 0; for (int i = 0; i < size; i++) { int rows = base + (i < rem ? 1 : 0); recvcounts[i] = rows * N; recvdispls[i] = offset; offset += rows * N; } } MPI_Gatherv(C_local, local_rows * N, MPI_DOUBLE, C, recvcounts, recvdispls, MPI_DOUBLE, 0, MPI_COMM_WORLD); if (rank == 0) { printf("C[0][0] = %f\n", C[0]); free(A); free(C); free(sendcounts); free(displs); free(recvcounts); free(recvdispls); } free(B); free(A_local); free(C_local); MPI_Finalize(); return 0; }

这段代码的逻辑说明:MPI_Bcast 把 B 从进程 0 广播到所有进程,这是整个程序里唯一一次全局通信。MPI_Scatterv 把 A 按行分给各进程,用 sendcounts 和 displs 处理不整除的情况。本地三重循环算 C_local,注意这里把 k 放在中间层,是为了让 B 的访问按行连续,缓存命中率更好。MPI_Gatherv 把各进程的 C_local 拼回进程 0。

参数说明:M、K、N 是矩阵维度,改这三个值就能测不同规模。base 和 rem 控制分块,rank < rem 的进程多拿一行。sendcounts 的单位是元素个数,不是行数,所以是 rows * K。displs 是元素偏移,不是行偏移。这两个单位搞混是新手最常见的翻车点。

编译和运行:

mpicc -O2 -o matmul_mpi matmul_mpi.c mpirun -np 4 ./matmul_mpi

如果机器上没有 mpicc,先确认 MPI 实现装好了。常见的是 MPICH 或 OpenMPI,用 mpicc --version 能看到版本信息。运行时的 -np 4 表示开 4 个进程,建议先设成 1 验证串行逻辑,再逐步加到 2、4、8。

3. 通信开销与性能拐点:什么时候 MPI 反而更慢

3.1 加速比不涨的三个原因

很多人跑完上面的代码,发现 4 进程比 1 进程快不了多少,甚至更慢。这不是代码写错了,而是矩阵乘法 MPI 的性能拐点到了。加速比上不去,通常有三个原因。

第一,通信量随进程数增长。按行切 A 时,B 的广播量是 K×N,这个量不随进程数变化,但每个进程都要接收完整 B。进程数越多,广播的同步开销越大。当 K×N 很大时,广播本身就要花掉可观时间。

第二,计算量被切得太碎。如果 M 只有 256,开 16 个进程,每个进程只算 16 行,本地计算时间极短,通信占比就上去了。经验值是每个进程至少分到 64 行以上,MPI 的收益才明显。

第三,内存带宽成为瓶颈。本地三重循环里,B 的每一行被反复读取。如果 B 放不进缓存,每次读 B 都要走内存,进程数再多也救不了。这时候要做的是分块乘法,把 B 切成小块,让每块能塞进 L2 或 L3 缓存。

3.2 用计时函数定位通信与计算占比

不要靠感觉判断瓶颈,直接加计时。MPI_Wtime 返回墙钟时间,精度足够。

double t_start = MPI_Wtime(); MPI_Bcast(B, K * N, MPI_DOUBLE, 0, MPI_COMM_WORLD); double t_after_bcast = MPI_Wtime(); // ... 本地计算 ... double t_after_compute = MPI_Wtime(); MPI_Gatherv(C_local, local_rows * N, MPI_DOUBLE, C, recvcounts, recvdispls, MPI_DOUBLE, 0, MPI_COMM_WORLD); double t_end = MPI_Wtime(); if (rank == 0) { printf("bcast: %.4f s\n", t_after_bcast - t_start); printf("compute: %.4f s\n", t_after_compute - t_after_bcast); printf("gather: %.4f s\n", t_end - t_after_compute); }

逻辑说明:把广播、计算、收集三段分别计时。如果 bcast 占比超过 20%,说明 B 太大或者进程数太多,考虑换二维分块减少广播量。如果 compute 占比高但加速比低,说明本地循环没优化好,先做缓存分块。如果 gather 占比高,说明 C 的收集方式有问题,可以改成各进程自己写文件,避开集中收集。

参数说明:MPI_Wtime 的单位是秒,返回 double。注意它测的是墙钟时间,不是 CPU 时间,所以能反映真实的等待开销。计时点要放在通信调用前后,不要放在循环内部,否则计时本身会干扰性能。

3.3 二维分块什么时候值得上

按行切 A 的瓶颈在 B 的广播。如果 K×N 是 4096×4096,双精度就是 128MB,广播一次在千兆网络上要一秒以上。这时候二维分块的优势就出来了:每个进程只需要 B 的一个子块,广播量降到原来的 1/sqrt(p)。

但二维分块的代价是逻辑复杂度。进程要映射成 sqrt(p) × sqrt(p) 的网格,A 按行广播、B 按列广播,各做一轮。进程网格不是正方形时,还要处理非方阵映射。我的建议是:先用按行切跑通,测出 bcast 占比。如果 bcast 超过总时间 30%,再考虑二维分块。否则按行切足够用,没必要为了理论上的通信优化把代码复杂度翻三倍。

4. 避坑与排查:矩阵乘法 MPI 最容易翻车的五个地方

4.1 现象:小矩阵结果对,大矩阵结果错

原因:分块时没有处理 M 不能被进程数整除的情况。小矩阵刚好整除,大矩阵有余数,某些进程少算或多算了几行,收集时错位。

解决:分块逻辑统一用 base + (rank < rem ? 1 : 0),sendcounts、displs、recvcounts、recvdispls 全部按这个规则算。写完先拿 M=10、p=3 这种不整除的组合测一遍,确认 C 的每一行都对。

4.2 现象:程序卡死,不报错也不退出

原因:MPI_Scatterv 或 MPI_Gatherv 的 counts 和 displs 在非 0 进程上传了 NULL,但 MPI 在某些实现里仍然会读这些参数。或者 Bcast 的 root 设错,所有进程都在等一个不存在的发送者。

解决:非 0 进程的 sendcounts、displs 可以传 NULL,但 recvcounts、recvdispls 在 Gatherv 里只有 root 需要,其他进程传 NULL 是合法的。如果卡死,先用 -np 1 跑,再 -np 2 跑,定位是哪个通信调用挂住。加 MPI_Barrier 分段,看卡在哪一段。

4.3 现象:加速比是负的,进程越多越慢

原因:矩阵太小,通信开销盖过计算收益。或者 B 的广播量太大,每个进程都在等网络。

解决:先算一下每个进程的计算量。如果 local_rows × K × N 小于 10^7 次浮点运算,MPI 基本没收益。把矩阵放大到每个进程至少 64 行,再看加速比。如果 B 太大,换二维分块或者把 B 也切分。

4.4 现象:结果有微小误差,和串行版本对不上

原因:浮点累加顺序变了。串行是 i→k→j 顺序累加,MPI 分块后每个进程算自己那部分,累加顺序不同,浮点舍入误差就不同。

解决:这是正常现象,不是 bug。用相对误差判断,|C_mpi - C_serial| / |C_serial| 小于 1e-10 就算对。如果误差大到 1e-3,那才是逻辑错误,回去查分块和收集的偏移。

4.5 现象:mpirun 报错 “not enough slots available”

原因:本机核心数不够,或者 MPI 配置的 slot 数小于 -np 指定的进程数。

解决:用 mpirun -np 4 --oversubscribe ./matmul_mpi 允许超订。或者检查机器实际核心数,nproc 看一下。超订会带来上下文切换开销,测性能时不要用,调试时可以临时开。

5. 从能跑到跑快:缓存分块与验证脚本的落地技巧

把 MPI 矩阵乘法跑通只是第一步。真正让它跑快,关键在本地计算那一层。上面代码里的三重循环,B 的访问是 k 行 j 列,每次内层循环读 B[k*N+j],如果 N 很大,B 的一行读不完就换行,缓存命中率很低。一个立竿见影的改法是做缓存分块:把 K 和 N 各切一块,比如 64×64,让 B 的子块能塞进 L1 或 L2。

#define BS 64 // 分块大小,按 CPU 缓存调 for (int ii = 0; ii < local_rows; ii += BS) { for (int kk = 0; kk < K; kk += BS) { for (int jj = 0; jj < N; jj += BS) { int i_max = (ii + BS < local_rows) ? ii + BS : local_rows; int k_max = (kk + BS < K) ? kk + BS : K; int j_max = (jj + BS < N) ? jj + BS : N; for (int i = ii; i < i_max; i++) { for (int k = kk; k < k_max; k++) { double a = A_local[i * K + k]; for (int j = jj; j < j_max; j++) { C_local[i * N + j] += a * B[k * N + j]; } } } } } }

逻辑说明:三层分块循环把大矩阵切成 BS×BS 的小块,每个小块在缓存里完成乘加。BS 取 64 是常见起点,具体值要看 CPU 的 L1/L2 大小。改完这个,单进程的本地计算时间通常能降 30% 到 50%,MPI 的加速比也会跟着好看。

参数说明:BS 不是越大越好。BS 太大,子块塞不进缓存,分块就白做了。BS 太小,循环开销占比高。建议从 32 开始,试 32、64、128 三档,看哪个最快。这个调参没有理论最优,跟机器相关,实测为准。

验证方面,不要只打印 C[0][0]。写一个串行版本,跑同样的 M、K、N,把两个结果逐元素比。下面这个校验片段可以直接嵌到主程序里:

// 只在 rank 0 执行,和串行结果对比 double max_err = 0.0; for (int i = 0; i < M * N; i++) { double diff = fabs(C[i] - C_serial[i]); double ref = fabs(C_serial[i]) > 1e-12 ? fabs(C_serial[i]) : 1e-12; double rel = diff / ref; if (rel > max_err) max_err = rel; } printf("max relative error: %.3e\n", max_err);

逻辑说明:逐元素算相对误差,取最大值。如果 max_err 在 1e-10 量级,说明 MPI 版本和串行版本一致。如果超过 1e-6,回去查分块偏移和收集逻辑。

参数说明:1e-12 是防止除零的兜底值。C_serial 是串行版本算出来的完整 C,只在 rank 0 上生成和比较。这个校验脚本建议每次改完分块或通信逻辑都跑一遍,比肉眼盯结果可靠得多。

最后一个习惯:每次改完代码,先用 -np 1 跑,确认串行逻辑没坏;再用 -np 2 跑,确认通信逻辑没坏;最后才上 -np 4、-np 8 测性能。跳过前两步直接上多进程,出了问题你连是计算错还是通信错都分不清。这个顺序帮我省过很多次后悔药。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询