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

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

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

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

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

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

![在这里插入图片描述][a91851cd469a4bdb90456f3a8eb24275.png_pic_center_x_16]

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

![在这里插入图片描述][2d45e894949b4bb6b84e15811c812ae0.png_pic_center_x_16]

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

![在这里插入图片描述][94f2838d5b624de5b28ff7f9d89947b6.png_pic_center_x_16]

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

#include <stdio.h>
    #include <mpi.h>
    #include<iostream>
    #include<vector>
    #include <cstdlib>
    using namespace std;
    
    int main(int argc, char* argv[])
    
    {
        
        int myrank, processNum; // myrank: 进程号; processNum: 进程总数
        char processor_name[MPI_MAX_PROCESSOR_NAME];
        int namelen;
    
        MPI_Init(&argc, &argv);
        MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
        MPI_Comm_size(MPI_COMM_WORLD, &processNum);
        MPI_Get_processor_name(processor_name, &namelen); // 获取处理器名称
    
        double startTime, endTime; // 开始时间和结束时间
    
        // 如果是 0 号进程，随机初始化矩阵数据
        int row = 3000, middleColumn = 200, column = 300; // 矩阵 A 的行数、矩阵 A 的列数和 B 的行数、矩阵 B 的行数
        if (myrank == 0) {
         // 地主进程
            vector<vector<int>> A(row, vector<int>(middleColumn));
            vector<vector<int>> B(middleColumn, vector<int>(column));
            // 随机初始化矩阵 A 和 B
            for (int i = 0; i < row; i++) {
        
                for (int j = 0; j < middleColumn; j++) {
        
    				A[i][j] = rand() % 10;
    			}
    		}
            for (int i = 0; i < middleColumn; i++) {
        
                for (int j = 0; j < column; j++) {
        
                    B[i][j] = rand() % 10;
                }
            }
    
            // 打印矩阵 A 和 B
      //      cout << "A:" << endl;
      //      for (int i = 0; i < row; i++) {
        
    		//	cout << "    ";
      //          for (int j = 0; j < middleColumn; j++) {
        
    		//		cout << A[i][j] << " ";
    		//	}
    		//	cout << endl;
    		//}
    
      //      cout << "B:" << endl;
      //      for (int i = 0; i < middleColumn; i++) {
        
      //          cout << "    ";
      //          for (int j = 0; j < column; j++) {
        
    		//		cout << B[i][j] << " ";
    		//	}
      //          cout << endl;
      //      }
    
    
            startTime = MPI_Wtime(); // 开始计时
            int farmersProcessNum = processNum - 1; // 农民进程数
            // 分发矩阵 A 和 B
            for (int i = 0; i < row; i++) {
        
                for (int j = 0; j < column; j++) {
        
                    int destProcess = (i * row + j) % farmersProcessNum + 1; // 目标进程号
                    MPI_Send(A[i].data(), middleColumn, MPI_INT, destProcess, 0, MPI_COMM_WORLD);
                    vector<int> BColumn(middleColumn);
                    for (int k = 0; k < middleColumn; k++) {
        
    					BColumn[k] = B[k][j];
    				}
                    MPI_Send(BColumn.data(), middleColumn, MPI_INT, destProcess, 0, MPI_COMM_WORLD);
                }
            }
        }
        else {
         // 农民进程 (进程号 > 0 的进程)
            int count = (row * column) / (processNum - 1);      // 每个进程计算的数量
            int remainder = (row * column) % (processNum - 1);  // 均分后多出来的待计算数
            if (myrank <= remainder) {
        
    			count++;
    		}
            
            for (int i = 0; i < count; i++) {
        
                // 接收矩阵 ARow 和 BColumn
                vector<int> ARow(middleColumn);
                vector<int> BColumn(middleColumn);
                MPI_Recv(ARow.data(), middleColumn, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                MPI_Recv(BColumn.data(), middleColumn, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                // 计算矩阵乘法
                int result = 0;
                for (int j = 0; j < middleColumn; j++) {
        
                    result += ARow[j] * BColumn[j];
                }
                // 发送计算结果
                MPI_Send(&result, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);
            }
        }
    
        // 收集计算结果
        if (myrank == 0) {
        
    		vector<vector<int>> C(row, vector<int>(column));
            for (int i = 0; i < row; i++) {
        
                for (int j = 0; j < column; j++) {
        
                    // 计算农民进程的数量
                    int farmersProcessNum = processNum - 1;
                    // 计算接收的进程号
                    int sourceProcess = (i * row + j) % farmersProcessNum + 1;
                    MPI_Recv(&C[i][j], 1, MPI_INT, sourceProcess, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                }
            }
            endTime = MPI_Wtime(); // 结束计时
    
            // 打印矩阵 C
            //cout << "C:" << endl;
            //for (int i = 0; i < row; i++) {
        
            //    cout << "    ";
            //    for (int j = 0; j < column; j++) {
        
            //        cout << C[i][j] << " ";
            //    }
            //    cout << endl;
            //}
    
            // 打印计算时间
            cout << "Time: " << endTime - startTime << "s" << endl;
        }
    
    
        MPI_Finalize();
    
        return 0;
    }

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

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

![在这里插入图片描述][93f798bcb0c244f09f0af3dbb5f661fc.png_pic_center_x_16]

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

![在这里插入图片描述][14057950f4914fbbb4c1b23c46ef4be1.png_pic_center_x_16]

### 1维数组表示2维数组 ###

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

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

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

#include<iostream>
    #include<vector>
    #include<mpi.h>
    
    using namespace std;
    
    // 输入：两个一维vector，表示矩阵，输出矩阵的row，middleColumn，column
    // 输出：一维vector，表示两个输入vector的乘积
    vector<int> matrixMultiplication(vector<int> A, vector<int> B, int row, int middleColumn, int column) {
        
    	vector<int> res(row * column);
    	for (int i = 0; i < row; i++) {
        
    		for (int j = 0; j < column; j++) {
        
    			int tmp = 0;
    			for (int k = 0; k < middleColumn; k++) {
        
    				// A[i][k] * B[k][j]
    				tmp += A[i * middleColumn + k] * B[k * column + j];
    			}
    			// res[i][j] = tmp
    			res[i * column + j] = tmp;
    		}
    	}
    	return res;
    }
    
    
    int main(int argc, char* argv[])
    
    {
        
    
    	int myrank, processNum;
    	char processor_name[MPI_MAX_PROCESSOR_NAME];
    	int namelen;
    
    	MPI_Init(&argc, &argv);
    	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    	MPI_Comm_size(MPI_COMM_WORLD, &processNum);
    	MPI_Get_processor_name(processor_name, &namelen);
    
    	// 如果是 0 号进程，随机初始化矩阵数据
    	int row = 3000, middleColumn = 2000, column = 3000;
    	int farmerProcessNum = processNum - 1; // farmer 进程的数量
    	int chunkSize = row / farmerProcessNum; // 每个 farmer 进程处理的行数
    	int chunkSizeRemainder = row % farmerProcessNum; // 均分后的余数
    	double startTime, endTime; // 记录开始和结束时间
    	if (myrank == 0) {
         // 地主进程
    		vector<int> A(row * middleColumn);
    		vector<int> B(middleColumn * column);
    		// 随机初始化矩阵 A 和 B
    		for (int i = 0; i < row * middleColumn; i++) {
        
    			A[i] = rand() % 10;
    		}
    		for (int i = 0; i < middleColumn * column; i++) {
        
    			B[i] = rand() % 10;
    		}
    
    		// 打印矩阵 A 和 B
    		//cout << "A:" << endl;
    		//for (int i = 0; i < row; i++) {
        
    		//	cout << "    ";
    		//	for (int j = 0; j < middleColumn; j++) {
        
    		//		cout << A[i * middleColumn + j] << " ";
    		//	}
    		//	cout << endl;
    		//}
    
    		//cout << "B:" << endl;
    		//for (int i = 0; i < middleColumn; i++) {
        
    		//	cout << "    ";
    		//	for (int j = 0; j < column; j++) {
        
    		//		cout << B[i * column + j] << " ";
    		//	}
    		//	cout << endl;
    		//}
    
    		// 记录程序开始时间
    		startTime = MPI_Wtime();
    
    		int startRow = 0, endRow = 0; // 每个 farmer 进程处理的起始行和结束行
    
    		// 将矩阵 A 和 B 分发给其他进程
    		for (int i = 1; i < processNum; i++) {
        
    			// 计算下一个进程处理的结束行
    			endRow = startRow + chunkSize - 1;
    			if (i <= chunkSizeRemainder) {
        
    				endRow++;
    			}
    
    			int tmpRow = endRow - startRow + 1;
    			// 将矩阵 A 部分的数据分发给其他进程
    			MPI_Send(A.data() + startRow * middleColumn, tmpRow * middleColumn, MPI_INT, i, 0, MPI_COMM_WORLD);
    			// 将矩阵 B 所有的数据分发给其他进程
    			MPI_Send(B.data(), middleColumn * column, MPI_INT, i, 0, MPI_COMM_WORLD);
    			// 更新下一个进程处理的起始行
    			startRow = endRow + 1;
    		}
    	}
    	else {
         // 农民进程 (进程号 >= 1)
    		// 计算行数
    		int tmpRow = chunkSize;
    
    		if (myrank <= chunkSizeRemainder) {
        
    			tmpRow++;
    		}
    
    		// 计算列数
    		int tmpColumn = column;
    		// 计算中间列数
    		int tmpMiddleColumn = middleColumn;
    
    		// 接收矩阵 A 的数据
    		vector<int> A(tmpRow * tmpMiddleColumn);
    		MPI_Recv(A.data(), tmpRow * tmpMiddleColumn, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    
    		// 接收矩阵 B 的数据
    		vector<int> B(tmpMiddleColumn * tmpColumn);
    		MPI_Recv(B.data(), tmpMiddleColumn * tmpColumn, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    
    		// 计算矩阵乘法
    		vector<int> res = matrixMultiplication(A, B, tmpRow, tmpMiddleColumn, tmpColumn);
    
    		// 将计算结果返回给 0 号进程
    		MPI_Send(res.data(), tmpRow * tmpColumn, MPI_INT, 0, 0, MPI_COMM_WORLD);
    	}
    
    	// 0 号进程接收其他进程的计算结果
    	if (myrank == 0) {
        
    		vector<int> res(row * column);
    		int startRow = 0, endRow = 0; // 每个 farmer 进程处理的起始行和结束行
    		for (int i = 1; i < processNum; i++) {
        
    			// 计算下一个进程处理的结束行
    			endRow = startRow + chunkSize - 1;
    			if (i <= chunkSizeRemainder) {
        
    				endRow++;
    			}
    			int tmpRow = endRow - startRow + 1;
    			// 将矩阵 A 部分的数据分发给其他进程
    			MPI_Recv(res.data() + startRow * column, tmpRow * column, MPI_INT, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    			// 更新下一个进程处理的起始行
    			startRow = endRow + 1;
    		}
    
    		// 记录程序结束时间
    		endTime = MPI_Wtime();
    		
    		// 打印计算结果
    		//cout << "res:" << endl;
    		//for (int i = 0; i < row; i++) {
        
    		//	cout << "    ";
    		//	for (int j = 0; j < column; j++) {
        
    		//		cout << res[i * column + j] << " ";
    		//	}
    		//	cout << endl;
    		//}
    
    		// 打印计算时间
    		cout << "Time: " << endTime - startTime << "s" << endl;
    	}
    
    	MPI_Finalize();
    
    	return 0;
    
    }

## 三、算法性能对比 ##

### 1. 时间复杂度 ###

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

![在这里插入图片描述][c2420d599eb54d8bbf2e8a7f213ae714.png_pic_center_x_16]

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

![在这里插入图片描述][113211f4598a421d86cbd470798ceb9f.png_pic_center_x_16]

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

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

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

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

#include<iostream>
    #include<vector>
    
    using namespace std;
    
    // 输入：两个一维vector，表示矩阵，输出矩阵的row，middleColumn，column
    // 输出：一维vector，表示两个输入vector的乘积
    vector<int> matrixMultiplication(vector<int> A, vector<int> B, int row, int middleColumn, int column) {
        
    	vector<int> res(row * column);
    	for (int i = 0; i < row; i++) {
        
    		for (int j = 0; j < column; j++) {
        
    			int tmp = 0;
    			for (int k = 0; k < middleColumn; k++) {
        
    				// A[i][k] * B[k][j]
    				tmp += A[i * middleColumn + k] * B[k * column + j];
    			}
    			// res[i][j] = tmp
    			res[i * column + j] = tmp;
    		}
    	}
    	return res;
    }
    
    int main() {
        
    	// 矩阵 A 的行数、矩阵 A 的列数和 B 的行数、矩阵 B 的行数
    	int row = 300, middleColumn = 200, column = 300;
    	vector<int> A(row * middleColumn);
    	vector<int> B(middleColumn * column);
    	// 随机初始化矩阵 A 和 B
    	for (int i = 0; i < row * middleColumn; i++) {
        
    		A[i] = rand() % 10;
    	}
    	for (int i = 0; i < middleColumn * column; i++) {
        
    		B[i] = rand() % 10;
    	}
    	// 记录开始和结束时间
    	double startTime, endTime;
    	startTime = clock();
    	vector<int> res = matrixMultiplication(A, B, row, middleColumn, column);
    	endTime = clock();
    	cout << "Time: " << (endTime - startTime) / CLOCKS_PER_SEC << "s" << endl;
    	// 输出 res
    	//for (int i = 0; i < row; i++) {
        
    	//	for (int j = 0; j < column; j++) {
        
    	//		cout << res[i * column + j] << " ";
    	//	}
    	//	cout << endl;
    	//}
    	//cout << endl;
    }

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

思路1运行时间：

![在这里插入图片描述][566ca0a34e9f4f9b9d5e7f028f7c9825.png_pic_center_x_16]

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

思路2运行时间：

![在这里插入图片描述][30cbb3c9d2684f94acb97de8b5ee82a3.png_pic_center_x_16]

串行运行时间：

![在这里插入图片描述][a78dfe546d354838a4451d686bb86a03.png_pic_center_x_16]

### 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)

具体转化结果如下：  
![在这里插入图片描述][732e897a7f254ee6bc03eecf11d4a70e.png_pic_center_x_16]

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

再来看下思路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)

![在这里插入图片描述][38ec3509257447d3bde64f49e32b70d6.png_pic_center_x_16]

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

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

[a91851cd469a4bdb90456f3a8eb24275.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/71784a1988b542199d3af3bcf2913d99.png
[2d45e894949b4bb6b84e15811c812ae0.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/e91d44dfd4d247e5a3aa0e43e59f43d8.png
[94f2838d5b624de5b28ff7f9d89947b6.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/0dc6a1989ed54b0096e393717199cf1d.png
[93f798bcb0c244f09f0af3dbb5f661fc.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/b911b85e86904f158d4e52c4764caa5c.png
[14057950f4914fbbb4c1b23c46ef4be1.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/d4510ed3d5874866a00b3b1b5f399128.png
[https_blog.csdn.net_weixin_44560698_article_details_120566680]: https://blog.csdn.net/weixin_44560698/article/details/120566680
[c2420d599eb54d8bbf2e8a7f213ae714.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/5b7c567fff9346bea34843d85c617b30.png
[113211f4598a421d86cbd470798ceb9f.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/1d25d4ee7fae4eb9b97cdda6b710d211.png
[65bf13939ada4f9baa43e56afbfe5d9c.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/2b17b5222dca474f9cee9dfa5974d9ad.png
[566ca0a34e9f4f9b9d5e7f028f7c9825.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/c612126831e94cc78dd2f09a3e36db93.png
[30cbb3c9d2684f94acb97de8b5ee82a3.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/a7867a3d7dce4a87ae02b055ede80780.png
[a78dfe546d354838a4451d686bb86a03.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/dcbe80dd9a5f40649d6433685ffeb95b.png
[732e897a7f254ee6bc03eecf11d4a70e.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/2cd1ccfea78a489ebbc3f3118e6995d7.png
[38ec3509257447d3bde64f49e32b70d6.png_pic_center_x_16]: https://image.dandelioncloud.cn/pgy_files/images/2024/03/17/39b7b450acad4a95a3f0b3fb8ec33676.png