并行矩阵乘法（C++ mpi 并行实现）

Myth丶恋晨 2024-03-17 10:53 152阅读 0赞

矩阵乘法有2种思路，我最先想到的是第一种思路，但是时间、空间复杂度都比较高。后面参考了一些资料，实现了第二种思路。

一、思路1：按行、列分块

矩阵乘法有一个很好的性质，就是结果矩阵的每个元素是不互相依赖的，因此我们可以很好地实现并行。

假如想要计算的矩阵如下：

在这里插入图片描述

如果有8个进程，0号进程用于数据分发，其它进程用于计算，具体的分块如下图所示：

在这里插入图片描述

但是这种思路实现后会报错 Out of memory，具体原因后面分析。

在这里插入图片描述

思路1因为只传递了单行的数组，而二维vector单行数组是连续的，因此这里直接用二维vector表示矩阵就可以，代码如下：

  1. #include <stdio.h>
  2. #include <mpi.h>
  3. #include<iostream>
  4. #include<vector>
  5. #include <cstdlib>
  6. using namespace std;
  8. int main(int argc, char* argv[])
  10. {
  12.     int myrank, processNum; // myrank: 进程号; processNum: 进程总数
  13.     char processor_name[MPI_MAX_PROCESSOR_NAME];
  14.     int namelen;
  16.     MPI_Init(&argc, &argv);
  17.     MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
  18.     MPI_Comm_size(MPI_COMM_WORLD, &processNum);
  19.     MPI_Get_processor_name(processor_name, &namelen); // 获取处理器名称
  21.     double startTime, endTime; // 开始时间和结束时间
  23.     // 如果是 0 号进程，随机初始化矩阵数据
  24.     int row = 3000, middleColumn = 200, column = 300; // 矩阵 A 的行数、矩阵 A 的列数和 B 的行数、矩阵 B 的行数
  25.     if (myrank == 0) {
  26.      // 地主进程
  27.         vector<vector<int>> A(row, vector<int>(middleColumn));
  28.         vector<vector<int>> B(middleColumn, vector<int>(column));
  29.         // 随机初始化矩阵 A 和 B
  30.         for (int i = 0; i < row; i++) {
  32.             for (int j = 0; j < middleColumn; j++) {
  34.                 A[i][j] = rand() % 10;
  35.             }
  36.         }
  37.         for (int i = 0; i < middleColumn; i++) {
  39.             for (int j = 0; j < column; j++) {
  41.                 B[i][j] = rand() % 10;
  42.             }
  43.         }
  45.         // 打印矩阵 A 和 B
  46.   //      cout << "A:" << endl;
  47.   //      for (int i = 0; i < row; i++) {
  49.         //    cout << "    ";
  50.   //          for (int j = 0; j < middleColumn; j++) {
  52.         //        cout << A[i][j] << " ";
  53.         //    }
  54.         //    cout << endl;
  55.         //}
  57.   //      cout << "B:" << endl;
  58.   //      for (int i = 0; i < middleColumn; i++) {
  60.   //          cout << "    ";
  61.   //          for (int j = 0; j < column; j++) {
  63.         //        cout << B[i][j] << " ";
  64.         //    }
  65.   //          cout << endl;
  66.   //      }
  69.         startTime = MPI_Wtime(); // 开始计时
  70.         int farmersProcessNum = processNum - 1; // 农民进程数
  71.         // 分发矩阵 A 和 B
  72.         for (int i = 0; i < row; i++) {
  74.             for (int j = 0; j < column; j++) {
  76.                 int destProcess = (i * row + j) % farmersProcessNum + 1; // 目标进程号
  77.                 MPI_Send(A[i].data(), middleColumn, MPI_INT, destProcess, 0, MPI_COMM_WORLD);
  78.                 vector<int> BColumn(middleColumn);
  79.                 for (int k = 0; k < middleColumn; k++) {
  81.                     BColumn[k] = B[k][j];
  82.                 }
  83.                 MPI_Send(BColumn.data(), middleColumn, MPI_INT, destProcess, 0, MPI_COMM_WORLD);
  84.             }
  85.         }
  86.     }
  87.     else {
  88.      // 农民进程 (进程号 > 0 的进程)
  89.         int count = (row * column) / (processNum - 1);      // 每个进程计算的数量
  90.         int remainder = (row * column) % (processNum - 1);  // 均分后多出来的待计算数
  91.         if (myrank <= remainder) {
  93.             count++;
  94.         }
  96.         for (int i = 0; i < count; i++) {
  98.             // 接收矩阵 ARow 和 BColumn
  99.             vector<int> ARow(middleColumn);
  100.             vector<int> BColumn(middleColumn);
  101.             MPI_Recv(ARow.data(), middleColumn, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
  102.             MPI_Recv(BColumn.data(), middleColumn, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
  103.             // 计算矩阵乘法
  104.             int result = 0;
  105.             for (int j = 0; j < middleColumn; j++) {
  107.                 result += ARow[j] * BColumn[j];
  108.             }
  109.             // 发送计算结果
  110.             MPI_Send(&result, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);
  111.         }
  112.     }
  114.     // 收集计算结果
  115.     if (myrank == 0) {
  117.         vector<vector<int>> C(row, vector<int>(column));
  118.         for (int i = 0; i < row; i++) {
  120.             for (int j = 0; j < column; j++) {
  122.                 // 计算农民进程的数量
  123.                 int farmersProcessNum = processNum - 1;
  124.                 // 计算接收的进程号
  125.                 int sourceProcess = (i * row + j) % farmersProcessNum + 1;
  126.                 MPI_Recv(&C[i][j], 1, MPI_INT, sourceProcess, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
  127.             }
  128.         }
  129.         endTime = MPI_Wtime(); // 结束计时
  131.         // 打印矩阵 C
  132.         //cout << "C:" << endl;
  133.         //for (int i = 0; i < row; i++) {
  135.         //    cout << "    ";
  136.         //    for (int j = 0; j < column; j++) {
  138.         //        cout << C[i][j] << " ";
  139.         //    }
  140.         //    cout << endl;
  141.         //}
  143.         // 打印计算时间
  144.         cout << "Time: " << endTime - startTime << "s" << endl;
  145.     }
  148.     MPI_Finalize();
  150.     return 0;
  151. }

二、思路2：矩阵分块乘法

如果我们想要计算的矩阵如下：

在这里插入图片描述

同时，如果我们采用4个进程来计算这个矩阵乘法，0号进程用于数据分发，其它进程用于计算，如图：

在这里插入图片描述

1维数组表示2维数组

思路2的算法实现，由于我们进行MPI_Send发送数组的时候，必须提供一个连续地址的数组的首地址，而思路2不同于思路1只发送1行数据，它是要发送多行的。

这就要求我们不可以使用二维数组来实现，因为二维数组不同行的首尾地址并不是连续的。具体见这篇文章：https://blog.csdn.net/weixin_44560698/article/details/120566680.

所以，这里我们使用1维数组，表示2维数组，代码如下：

  1. #include<iostream>
  2. #include<vector>
  3. #include<mpi.h>
  5. using namespace std;
  7. // 输入：两个一维vector，表示矩阵，输出矩阵的row，middleColumn，column
  8. // 输出：一维vector，表示两个输入vector的乘积
  9. vector<int> matrixMultiplication(vector<int> A, vector<int> B, int row, int middleColumn, int column) {
  11.     vector<int> res(row * column);
  12.     for (int i = 0; i < row; i++) {
  14.         for (int j = 0; j < column; j++) {
  16.             int tmp = 0;
  17.             for (int k = 0; k < middleColumn; k++) {
  19.                 // A[i][k] * B[k][j]
  20.                 tmp += A[i * middleColumn + k] * B[k * column + j];
  21.             }
  22.             // res[i][j] = tmp
  23.             res[i * column + j] = tmp;
  24.         }
  25.     }
  26.     return res;
  27. }
  30. int main(int argc, char* argv[])
  32. {
  35.     int myrank, processNum;
  36.     char processor_name[MPI_MAX_PROCESSOR_NAME];
  37.     int namelen;
  39.     MPI_Init(&argc, &argv);
  40.     MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
  41.     MPI_Comm_size(MPI_COMM_WORLD, &processNum);
  42.     MPI_Get_processor_name(processor_name, &namelen);
  44.     // 如果是 0 号进程，随机初始化矩阵数据
  45.     int row = 3000, middleColumn = 2000, column = 3000;
  46.     int farmerProcessNum = processNum - 1; // farmer 进程的数量
  47.     int chunkSize = row / farmerProcessNum; // 每个 farmer 进程处理的行数
  48.     int chunkSizeRemainder = row % farmerProcessNum; // 均分后的余数
  49.     double startTime, endTime; // 记录开始和结束时间
  50.     if (myrank == 0) {
  51.      // 地主进程
  52.         vector<int> A(row * middleColumn);
  53.         vector<int> B(middleColumn * column);
  54.         // 随机初始化矩阵 A 和 B
  55.         for (int i = 0; i < row * middleColumn; i++) {
  57.             A[i] = rand() % 10;
  58.         }
  59.         for (int i = 0; i < middleColumn * column; i++) {
  61.             B[i] = rand() % 10;
  62.         }
  64.         // 打印矩阵 A 和 B
  65.         //cout << "A:" << endl;
  66.         //for (int i = 0; i < row; i++) {
  68.         //    cout << "    ";
  69.         //    for (int j = 0; j < middleColumn; j++) {
  71.         //        cout << A[i * middleColumn + j] << " ";
  72.         //    }
  73.         //    cout << endl;
  74.         //}
  76.         //cout << "B:" << endl;
  77.         //for (int i = 0; i < middleColumn; i++) {
  79.         //    cout << "    ";
  80.         //    for (int j = 0; j < column; j++) {
  82.         //        cout << B[i * column + j] << " ";
  83.         //    }
  84.         //    cout << endl;
  85.         //}
  87.         // 记录程序开始时间
  88.         startTime = MPI_Wtime();
  90.         int startRow = 0, endRow = 0; // 每个 farmer 进程处理的起始行和结束行
  92.         // 将矩阵 A 和 B 分发给其他进程
  93.         for (int i = 1; i < processNum; i++) {
  95.             // 计算下一个进程处理的结束行
  96.             endRow = startRow + chunkSize - 1;
  97.             if (i <= chunkSizeRemainder) {
  99.                 endRow++;
  100.             }
  102.             int tmpRow = endRow - startRow + 1;
  103.             // 将矩阵 A 部分的数据分发给其他进程
  104.             MPI_Send(A.data() + startRow * middleColumn, tmpRow * middleColumn, MPI_INT, i, 0, MPI_COMM_WORLD);
  105.             // 将矩阵 B 所有的数据分发给其他进程
  106.             MPI_Send(B.data(), middleColumn * column, MPI_INT, i, 0, MPI_COMM_WORLD);
  107.             // 更新下一个进程处理的起始行
  108.             startRow = endRow + 1;
  109.         }
  110.     }
  111.     else {
  112.      // 农民进程 (进程号 >= 1)
  113.         // 计算行数
  114.         int tmpRow = chunkSize;
  116.         if (myrank <= chunkSizeRemainder) {
  118.             tmpRow++;
  119.         }
  121.         // 计算列数
  122.         int tmpColumn = column;
  123.         // 计算中间列数
  124.         int tmpMiddleColumn = middleColumn;
  126.         // 接收矩阵 A 的数据
  127.         vector<int> A(tmpRow * tmpMiddleColumn);
  128.         MPI_Recv(A.data(), tmpRow * tmpMiddleColumn, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
  130.         // 接收矩阵 B 的数据
  131.         vector<int> B(tmpMiddleColumn * tmpColumn);
  132.         MPI_Recv(B.data(), tmpMiddleColumn * tmpColumn, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
  134.         // 计算矩阵乘法
  135.         vector<int> res = matrixMultiplication(A, B, tmpRow, tmpMiddleColumn, tmpColumn);
  137.         // 将计算结果返回给 0 号进程
  138.         MPI_Send(res.data(), tmpRow * tmpColumn, MPI_INT, 0, 0, MPI_COMM_WORLD);
  139.     }
  141.     // 0 号进程接收其他进程的计算结果
  142.     if (myrank == 0) {
  144.         vector<int> res(row * column);
  145.         int startRow = 0, endRow = 0; // 每个 farmer 进程处理的起始行和结束行
  146.         for (int i = 1; i < processNum; i++) {
  148.             // 计算下一个进程处理的结束行
  149.             endRow = startRow + chunkSize - 1;
  150.             if (i <= chunkSizeRemainder) {
  152.                 endRow++;
  153.             }
  154.             int tmpRow = endRow - startRow + 1;
  155.             // 将矩阵 A 部分的数据分发给其他进程
  156.             MPI_Recv(res.data() + startRow * column, tmpRow * column, MPI_INT, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
  157.             // 更新下一个进程处理的起始行
  158.             startRow = endRow + 1;
  159.         }
  161.         // 记录程序结束时间
  162.         endTime = MPI_Wtime();
  164.         // 打印计算结果
  165.         //cout << "res:" << endl;
  166.         //for (int i = 0; i < row; i++) {
  168.         //    cout << "    ";
  169.         //    for (int j = 0; j < column; j++) {
  171.         //        cout << res[i * column + j] << " ";
  172.         //    }
  173.         //    cout << endl;
  174.         //}
  176.         // 打印计算时间
  177.         cout << "Time: " << endTime - startTime << "s" << endl;
  178.     }
  180.     MPI_Finalize();
  182.     return 0;
  184. }

三、算法性能对比

1. 时间复杂度

思路一 300 × 200 300\times 200 300×200 和 200 × 300 200 \times 300 200×300 相乘所耗费的时间为8.83s。

在这里插入图片描述

相同规模的矩阵相乘，思路2的时间是0.07s。

在这里插入图片描述

可以看到差距还是很大的。

我们再来看一下串行程序的运行时间：

在这里插入图片描述
串行程序的运行时间是0.321s，看来如果我们采用思路一的并行方法就真有点得不偿失了。

这里给出我写的串行代码：

  1. #include<iostream>
  2. #include<vector>
  4. using namespace std;
  6. // 输入：两个一维vector，表示矩阵，输出矩阵的row，middleColumn，column
  7. // 输出：一维vector，表示两个输入vector的乘积
  8. vector<int> matrixMultiplication(vector<int> A, vector<int> B, int row, int middleColumn, int column) {
  10.     vector<int> res(row * column);
  11.     for (int i = 0; i < row; i++) {
  13.         for (int j = 0; j < column; j++) {
  15.             int tmp = 0;
  16.             for (int k = 0; k < middleColumn; k++) {
  18.                 // A[i][k] * B[k][j]
  19.                 tmp += A[i * middleColumn + k] * B[k * column + j];
  20.             }
  21.             // res[i][j] = tmp
  22.             res[i * column + j] = tmp;
  23.         }
  24.     }
  25.     return res;
  26. }
  28. int main() {
  30.     // 矩阵 A 的行数、矩阵 A 的列数和 B 的行数、矩阵 B 的行数
  31.     int row = 300, middleColumn = 200, column = 300;
  32.     vector<int> A(row * middleColumn);
  33.     vector<int> B(middleColumn * column);
  34.     // 随机初始化矩阵 A 和 B
  35.     for (int i = 0; i < row * middleColumn; i++) {
  37.         A[i] = rand() % 10;
  38.     }
  39.     for (int i = 0; i < middleColumn * column; i++) {
  41.         B[i] = rand() % 10;
  42.     }
  43.     // 记录开始和结束时间
  44.     double startTime, endTime;
  45.     startTime = clock();
  46.     vector<int> res = matrixMultiplication(A, B, row, middleColumn, column);
  47.     endTime = clock();
  48.     cout << "Time: " << (endTime - startTime) / CLOCKS_PER_SEC << "s" << endl;
  49.     // 输出 res
  50.     //for (int i = 0; i < row; i++) {
  52.     //    for (int j = 0; j < column; j++) {
  54.     //        cout << res[i * column + j] << " ";
  55.     //    }
  56.     //    cout << endl;
  57.     //}
  58.     //cout << endl;
  59. }

我们将矩阵的大小扩大10倍，再来看：

思路1运行时间：

在这里插入图片描述

可以看到思路1直接out of memory了，其实是在意料之中的，具体原因可以看文章最后一部分 四、复杂度分析 。

思路2运行时间：

在这里插入图片描述

串行运行时间：

在这里插入图片描述

2. 空间复杂度

四、复杂度分析

首先，我们设：

r o w : 矩阵 A 行数 row:矩阵A行数 row:矩阵A行数

m i d d l e C o l u m n : 矩阵 A 列数和矩阵 B 行数 middleColumn:矩阵A列数和矩阵B行数 middleColumn:矩阵A列数和矩阵B行数

c o l u m n : 矩阵 B 列数 column:矩阵B列数 column:矩阵B列数

f a r m e r P r o c e s s N u m : 农民进程数 farmerProcessNum:农民进程数 farmerProcessNum:农民进程数

由于mpi主要是在并行程序之间进行通信。

对于思路1，是每次计算结果矩阵的一个元素，因此，它的总发送的数据量如下：

r o w × c o l u m n × m i d d l e C o l u m n 2 ( i n t ) row\times column \times middleColumn^2\quad (int) row×column×middleColumn2(int)

对于思路2，总发送的是A矩阵的所有+B矩阵的全部乘以进程数(因为每次都要发送B矩阵)，如下：

r o w × m i d d l e C o l u m n + f a r m e r P r o c e s s N u m × m i d d l e C o l u m n × c o l u m n row\times middleColumn+farmerProcessNum\times middleColumn\times column row×middleColumn+farmerProcessNum×middleColumn×column

我们带入 r o w = 3000 , m i d d l e C o l u m n = 2000 , c o l u m n = 3000 row=3000,middleColumn=2000,column=3000 row=3000,middleColumn=2000,column=3000来计算一下上述思路1需要传输的内存：
r o w × c o l u m n × m i d d l e C o l u m n 2 ( i n t ) = 3000 × 3000 × 2000 × 2000 ( i n t ) = 36 0000 0000 0000 ( i n t ) = 4 × 36 0000 0000 0000 ( B y t e ) ≈ 131 ( T B ) \begin{aligned} &row\times column \times middleColumn^2\quad (int) \\ =&3000\times 3000 \times 2000\times 2000\quad (int) \\ =&36\quad 0000\quad 0000\quad 0000\quad (int) \\ =&4\times 36 \quad 0000\quad 0000\quad 0000\quad (Byte)\\ \approx&131 \quad (TB) \end{aligned} ===≈row×column×middleColumn2(int)3000×3000×2000×2000(int)36000000000000(int)4×36000000000000(Byte)131(TB)

具体转化结果如下：
在这里插入图片描述

传输的数据不管是传给哪个进程都是要放在内存里的，因此这种方式的空间数量级根本不可取= =。

再来看下思路2：

r o w × m i d d l e C o l u m n + f a r m e r P r o c e s s N u m × m i d d l e C o l u m n × c o l u m n ( i n t ) = 3000 × 2000 + 7 × 2000 × 3000 ( i n t ) = 4800 0000 ( i n t ) = 19200 0000 ( B y t e ) \begin{aligned} &row\times middleColumn+farmerProcessNum\times middleColumn\times column \quad (int) \\ =&3000\times 2000+7 \times 2000 \times 3000 \quad (int) \\ =&4800\quad 0000 \quad (int) \\ =&19200 \quad 0000 \quad (Byte) \end{aligned} ===row×middleColumn+farmerProcessNum×middleColumn×column(int)3000×2000+7×2000×3000(int)48000000(int)192000000(Byte)

在这里插入图片描述

可以看到 思路2 消耗的内存仅有 183MB ，与思路1 的 131TB 的空间复杂度完全不在一个数量级。

结论：所以对于矩阵并行乘法，思路1是不可行的，至少也是应该是采用思路2来进行实现的。

发表评论

表情：
评论列表 （有 0 条评论，152人围观）

还没有评论，来说两句吧...

相关阅读