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

简述

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

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

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

这里的方法二，

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

数学表达式

解法一的数学表达

Thread(i):∑j=0n−1Aij∗XjThread(i):\sum_{j=0}^{n-1}{A_{ij}*X_j}Thread(i):j=0∑n−1Aij∗Xj
result={j∈{0,1,..,n−1}∣Thread(j)}result = {\{j\in\{0,1,..,n-1\}| Thread(j)\}}result={j∈{0,1,..,n−1}∣Thread(j)}

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

解法二的数学表达式

Thread(i):{j∈{0,1,..,n−1}∣Aij∗Xi}Thread(i):{\{j\in\{0,1,..,n-1\}| A_{ij}*X_i\}}Thread(i):{j∈{0,1,..,n−1}∣Aij∗Xi}

result=∑i=0n−1Thread(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;// initializeif (my_rank == 0) {a = new double*[N];for (int i = 0; i < N; ++i) {a[i] = new double[N];}// read afor (int i = 0; i < N; ++i)for (int j = 0; j < N; ++j)cin >> a[j][i]; // transpose ks = new double[N];// read ksfor (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 NMPI_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;// finalizeMPI_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

计算过程

HINTS

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

关于MPI_Reduce : 我会之后补一份介绍。