【MPI编程】矩阵向量乘法--解法三（子矩阵块分解）【高性能计算】

桃扇骨 2022-04-13 03:21 250阅读 0赞

## 简述 ##

子矩阵分解，就是说，将原来的矩阵给分解为更小的矩阵块。

*  让所有的线程都共享有向量（但不共享矩阵）

为了试验简单，这里做了几个简单的假设

*  矩阵为方阵(n\*n)
 *  comm\_sz（线程数量）和n满足下面的关系
    
     *  n被 c o m m \_ s z \\sqrt\{comm\\\_sz\} comm\_sz整除
     *  comm\_sz为平方数

我们基于之前的版本改进

*  [【MPI编程】矩阵向量乘法–解法二（高性能计算）][MPI]
 *  数据也是用一样的。（所以我们知道结果也应该是一样的。）

## 算法分析 ##

*  用一个一维数据来记录数据
 *  线程0来读取数据，并分发给其他线程。最后大家一起来算，最后再汇总给线程0。
 *  每个线程算自己对应的矩阵，和部分的向量的内积。

**汇总的过程：**

*  汇总的过程中，处于同一行的其他线程将数据发给同行的第一个
 *  最后第0列的所有元素按顺序发给线程0。之后，按顺序输出数据。

**难点：**

*  定位上，通过一维数组来表示，便于传递数据，但是在设计程序上稍显麻烦。可以用宏定义的方式来减轻负担。不过这次的代码只有160行。**不太需要这些操作。（可以用来锻炼思维）**

## 代码 ##

### 在VS上可以编译代码 ###

#include<stdio.h>
    #include<string.h>
    #include<mpi.h>
    #pragma warning(disable : 4996)
    #define MAX_STRING 100
    using namespace std;
    #include <fstream>
    #include <iostream>
    #include <cmath>
    
    int main(void) { 
    	int comm_sz;
    	int my_rank;
    	MPI_Init(NULL, NULL);
    	MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
    	MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
    
    
    	// 只有一个线程的时候不操作
    	if (comm_sz <= 1) { 
    		MPI_Finalize();
    		return 0;
    	}
    
    	int N = 0, local_N, sqrtCommSZ = sqrt(comm_sz);
    	double *b = NULL, *local, *a, *ks, *local_a, *local_ans;
    	int rank_len = 0, local_len;
    	// initialize
    	if (my_rank == 0) { 
    		ifstream cin("D:\\C++\\VS\\repo\\MPI-DEMO\\MPI-DEMO\\input.txt");
    		cin >> N;
    		// Send N to all 
    		MPI_Bcast(&N, 1, MPI_INT, 0, MPI_COMM_WORLD);
    		// create a
    		a = new double[N*N];
    		// read a
    		for (int i = 0; i < N; ++i)
    			for (int j = 0; j < N; ++j)
    				cin >> a[i * N + j];
    		// create ks
    		ks = new double[N];
    		// read ks
    		for (int i = 0; i < N; ++i) cin >> ks[i];
    		// send data to other.
    		MPI_Bcast(ks, N, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    		local_N = N * N / comm_sz;
    		local = new double[local_N];
    		local_a = new double[local_N];
    		local_len = N / sqrtCommSZ;
    		for (int i = 0; i < sqrtCommSZ; ++i) { 
    			for (int j = 0; j < sqrtCommSZ; ++j) { 
    				if (i == 0 && j == 0) {  // record from Mat a
    					for (int k = 0; k < local_len; ++k) { 
    						for (int t = 0; t < local_len; ++t) { 
    							local_a[k * local_len + t] =
    								a[k * N + t];
    						}
    					}
    				}
    				else {  // copy and send!
    					// (i * sqrtCommSZ + j) is the target rank.
    					for (int k = 0; k < local_len; ++k) { 
    						for (int t = 0; t < local_len; ++t) { 
    							local[k * local_len + t] = // keep 
    								a[i * local_len * N + j * local_len // 定位到左上角
    								+ k * N + t];
    							// 左上角坐标：[i * local_len, j * local_len]
    						}
    					}
    					MPI_Send(local, local_len * local_len, MPI_DOUBLE,
    						i * sqrtCommSZ + j, 0, MPI_COMM_WORLD);
    				}
    			}
    		}
    		local_ans = new double[local_len];
    		for (int i = 0; i < local_len; ++i) { 
    			local_ans[i] = 0;
    			for (int j = 0; j < local_len; ++j) { 
    				local_ans[i] += local_a[i * local_len + j] * ks[j];
    			}
    		}
    
    		// recv from the same row
    		double * temp_local = new double[local_len];
    		for (int i = 1; i < sqrtCommSZ; ++i) { 
    			MPI_Recv(temp_local, local_len, MPI_DOUBLE,
    				my_rank + i, 1, MPI_COMM_WORLD, MPI_STATUSES_IGNORE);
    			for (int j = 0; j < local_len; ++j) { 
    				local_ans[j] += temp_local[j];
    			}
    		}
    		// recv from different rows
    
    		for (int i = 0; i < local_len; ++i) { 
    			cout << local_ans[i] << " ";
    		}
    
    		for (int i = 1; i < sqrtCommSZ; ++i) { 
    			MPI_Recv(temp_local, local_len, MPI_DOUBLE,
    				my_rank + i * sqrtCommSZ, 2, MPI_COMM_WORLD, MPI_STATUSES_IGNORE);
    			for (int j = 0; j < local_len; ++j) { 
    				cout << temp_local[j] << " ";
    			}
    		}
    		cout << endl;
    
    		delete[] temp_local;
    		delete[]a;
    		delete[]local;
    	}
    	else {  // other rank
    		// recv data
    		MPI_Bcast(&N, 1, MPI_INT, 0, MPI_COMM_WORLD);
    		ks = new double[N];
    		MPI_Bcast(ks, N, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    		local_len = N / sqrtCommSZ;
    		local_a = new double[local_len * local_len];
    		MPI_Recv(local_a, local_len * local_len, MPI_DOUBLE, 
    			0, 0, MPI_COMM_WORLD, MPI_STATUSES_IGNORE);
    		// 本地做乘法
    		local_ans = new double[local_len];
    		for (int i = 0; i < local_len; ++i) { 
    			local_ans[i] = 0;
    			for (int j = 0; j < local_len; ++j) { 
    				local_ans[i] += local_a[i * local_len + j] * 
    					ks[j + local_len * (my_rank % sqrtCommSZ)];
    			}
    		}
    
    		
    
    		if (my_rank % sqrtCommSZ == 0) {  // recv and send to 0,0
    			double * temp_local = new double[local_len];
    			for (int i = 1; i < sqrtCommSZ; ++i) { 
    				MPI_Recv(temp_local, local_len, MPI_DOUBLE,
    					my_rank + i, 1, MPI_COMM_WORLD, MPI_STATUSES_IGNORE);
    				for (int j = 0; j < local_len; ++j) { 
    					local_ans[j] += temp_local[j];
    				}
    			}
    			delete[] temp_local;
    			// send to 0, 0
    			MPI_Send(local_ans, local_len, MPI_DOUBLE,
    				0, 2, MPI_COMM_WORLD);
    		}
    		else {  // send
    			MPI_Send(local_ans, local_len, MPI_DOUBLE,
    				my_rank - my_rank % sqrtCommSZ, 1, MPI_COMM_WORLD);
    		}
    		
    	}
    	delete[] local_ans;
    	delete[]ks;
    	delete[] local_a;
    
    	// finalize
    	MPI_Finalize();
    
    	return 0;
    }

### 测试数据 ###

10
    3 10 6 0 8 5 6 4 1 6 
    0 1 10 7 2 10 8 10 3 2 
    6 1 10 10 2 5 5 5 5 2 
    3 7 1 2 7 2 3 8 2 7 
    3 7 8 4 8 7 4 4 4 1 
    5 4 8 3 4 10 7 5 5 0 
    2 4 7 10 5 7 9 6 5 10 
    9 9 1 4 10 0 1 8 1 5 
    8 0 3 2 10 9 10 8 10 7 
    6 9 0 4 4 9 1 0 0 6 
    2 1 0 2 5 8 7 4 6 4

### 测试结果 ###

![在这里插入图片描述][20181126225316591.png]

很明显跟之前的结果是一样的。

### Hints ###

*  这里用到了广播。就是将一份数据copy之后，发给所有的其他线程上。
 *  可以参照下面的链接
 *  [【MPI编程】MPI\_Bcast广播讲解和使用][MPI_MPI_Bcast]

[MPI]: https://blog.csdn.net/a19990412/article/details/84553635
[20181126225316591.png]: /images/20220413/9dd328277981480ea73492af22ae329b.png
[MPI_MPI_Bcast]: https://blog.csdn.net/a19990412/article/details/84283205