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

朱雀 2022-04-22 14:44 193阅读 0赞

## 简述 ##

有留心过的朋友可能会发现，其实我没写过解法一。

因为解法一就是大家最直观的感受的一种解法。

*  将矩阵按照行划分，之后，再每个线程都用整个向量跟这个块做内积。之后，传回来一个数值。

这里的方法二，

*  将矩阵按照列进行划分。然后将列向量传递给所有的线程，之后，再传对应的向量的对应的元素，给这个线程。之后，这个元素（数值），乘上传给它的向量的每一个元素。最后，再做一个Reduce的操作，将所有的向量加起来。就得到了一个新的向量。这个向量就是解。

## 数学表达式 ##

### 解法一的数学表达 ###

T h r e a d ( i ) : ∑ j = 0 n − 1 A i j ∗ X j Thread(i):\\sum\_\{j=0\}^\{n-1\}\{A\_\{ij\}\*X\_j\} Thread(i):j=0∑n−1Aij∗Xj  
 r e s u l t = \{ j ∈ \{ 0 , 1 , . . , n − 1 \} ∣ T h r e a d ( j ) \} result = \{\\\{j\\in\\\{0,1,..,n-1\\\}| Thread(j)\\\}\} result=\{ j∈\{ 0,1,..,n−1\}∣Thread(j)\}

*  第一步得到的是一个数值
 *  result的过程是将数值组合成向量的过程。

### 解法二的数学表达式 ###

T h r e a d ( i ) : \{ j ∈ \{ 0 , 1 , . . , n − 1 \} ∣ A i j ∗ X i \} Thread(i):\{\\\{j\\in\\\{0,1,..,n-1\\\}| A\_\{ij\}\*X\_i\\\}\} Thread(i):\{ j∈\{ 0,1,..,n−1\}∣Aij∗Xi\}

r e s u l t = ∑ i = 0 n − 1 T h r e a d ( i ) result = \\sum\_\{i=0\}^\{n-1\}\{Thread(i)\} result=i=0∑n−1Thread(i)

*  第一步得到的是一个向量
 *  result求和的过程是向量的求和过程。

从数学表达式上，可以很清楚地知道这个结果是一致的。（这里假设了这个矩阵是n\*n）

## 在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>
    
    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;
    	}
    	ifstream cin("D:\\C++\\VS\\repo\\MPI-DEMO\\MPI-DEMO\\input.txt");
    	int N = 0;
    	double *b = NULL, *local, **a, *ks, k, *local_arr;
    	cin >> N;
    	int rank_len = 0;
    
    	// initialize
    	if (my_rank == 0) { 
    		a = new double*[N];
    		for (int i = 0; i < N; ++i) { 
    			a[i] = new double[N];
    		}
    		// read a
    		for (int i = 0; i < N; ++i)
    			for (int j = 0; j < N; ++j)
    				cin >> a[j][i]; // transpose 
    		ks = new double[N];
    		// read ks
    		for (int i = 0; i < N; ++i) cin >> ks[i];
    		// send data to other.
    		rank_len = N / comm_sz;
    		for (int i = 1; i < comm_sz; ++i) { 
    			// send N
    			MPI_Send(&N, 1, MPI_INT, i, 0, MPI_COMM_WORLD);
    			for (int j = 0; j < rank_len; ++j) { 
    				MPI_Send(a[i * rank_len + j], N, MPI_DOUBLE, i, 2 * j + 1, MPI_COMM_WORLD);
    				MPI_Send(&ks[i * rank_len + j], 1, MPI_DOUBLE, i, 2 * j + 2, MPI_COMM_WORLD);
    			}
    		}
    		local_arr = new double[N];
    		local = new double[N];
    		for (int i = 0; i < N; ++i) 
    			local_arr[i] = 0;
    		for (int i = 0; i < rank_len; ++i) { 
    			for (int j = 0; j < N; ++j) { 
    				local_arr[j] += (a[i][j] * ks[i]);
    			}
    		}
    		MPI_Reduce(local_arr, local, N, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
    		for (int i = 0; i < N; ++i) cout << local[i] << " ";
    		cout << endl;
    
    		// release 
    		for (int i = 0; i < N; ++i) delete[]a[i];
    		delete[]a;
    		delete[]ks;
    	}
    	else { 
    		MPI_Recv(&N, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, MPI_STATUSES_IGNORE);
    		local = new double[N];
    		local_arr = new double[N];
    		for (int i = 0; i < N; ++i) local_arr[i] = 0;
    		rank_len = N / comm_sz;
    		for (int j = 0; j < rank_len; ++j) { 
    			MPI_Recv(local, N, MPI_DOUBLE, 0, 2 * j + 1, MPI_COMM_WORLD, MPI_STATUSES_IGNORE);
    			MPI_Recv(&k, 1, MPI_DOUBLE, 0, 2 * j + 2, MPI_COMM_WORLD, MPI_STATUSES_IGNORE);
    			for (int i = 0; i < N; ++i) local_arr[i] += (k * local[i]);
    		}
    		MPI_Reduce(local_arr, NULL, N, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
    	}
    
    	delete[]local_arr;
    	delete[] local;
    
    	// 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

### 计算过程 ###

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

### HINTS ###

**注意，这里调用的线程数量，必须可以整除n**

*  关于MPI\_Reduce : 我会之后补一份介绍。

[20181126174851574.png]: /images/20220414/2ae81daa96944e2abe35dcec4abb63c2.png