# 并行程序设计实验报告
陈翊辉
SA19011116
[toc]
## 公共部分
### 并行求和
#### 蝶式求和
##### 并行算法描述
蝶形运算是并行计算中一种基本的方法，快速傅里叶变换就通常使用蝶形运算，而蝶形求和也是最简单的蝶形计算，可以使$2^n$个处理器经过$n$步得到全和。其并行算法如下：
```pascal
输入：n个数，每个处理器上各有一个，n是2的幂，n=2^M
输出：n个数的和，每个处理器上都有
Begin
for k=1 to m do
计算对应处理器号t
按蝶形计算向对应处理器t发送自己的部分和（第一次是初始数）
按蝶形计算从对应处理器t接收部分和并加到自己部分和（第一次是初始数）
end for
End
```
其中对应处理器号`t`计算
```c++
base = 2 ** k;
group = rank / base;
offset = rank % base;
target = (group % 2 ? group - 1 : group + 1) * base + offset;
```
##### MPI实现
在蝶形运算中，共有$log_2n$步，每一步需要每个处理器给对应处理器发送，并从对应处理器接收，如果都先发或先收则必定造成死锁，这里用简单的rank和target比较判断，也可以考虑异步收发或者用`MPI_Sendrecv`。
##### 核心MPI代码
```c++
for (int i = 0; i < steps; ++i) {
int group = rank / base;
int offset = rank % base;
int target = (group % 2 ? group - 1 : group + 1) * base + offset;
if (rank < target) {
MPI_Send(reinterpret_cast(&send), 1, MPI_INT, target, target, MPI_COMM_WORLD);
MPI_Recv(reinterpret_cast(&recv), 1, MPI_INT, target, rank, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
} else {
MPI_Recv(reinterpret_cast(&recv), 1, MPI_INT, target, rank, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
MPI_Send(reinterpret_cast(&send), 1, MPI_INT, target, target, MPI_COMM_WORLD);
}
send += recv;
base *= 2;
}
```
*完整代码见附件*
##### 性能结果
| 问题规模 | 任务数 | 总运行时间(sec) |
| -------- | ------------- | --------------- |
| 4 | 4 | 0.021 |
| 8 | 8 | 0.025 |
| 16 | 16 | 0.020 |
| 32 | 32 | 0.168 |
| 64 | 64（跨机器） | 0.397 |
| 128 | 128（跨机器） | 0.447 |
#### 二叉树求和
##### 并行算法描述
二叉树求和使$2^n$个处理器经过$2n$步得到全和，前$n$步是树上的归约求和过程，后$n$步是发送结果过程。其并行算法如下：
```pascal
输入：n个数，每个处理器上各有一个，n是2的幂，n=2^M
输出：n个数的和，每个处理器上都有
Begin
for k=1 to m do
叶子结点向根节点发送部分和（第一次是初始数）
根结点求和，根结点继续下一轮，非根结点结束
end for
for k=1 to m do
根结点向两个孩子结点发送全和
end for
End
```
二叉树根的选取可以任意，一种比较简单的二叉树构建方法是按处理器编号0,1;2,3;...;构建，每次编号小的作为根结点，然后进行下一轮运算。
##### MPI实现
一种比较简单的确定二叉树方法是按处理器编号0,1;2,3;...;构建，每次编号小的作为根结点，然后进行下一轮运算。比如开始时，所有结点都参与，即mod1=0的参与，mod2=0的作为根结点；第2步，mod2=0的参与，mod4=0的作为根结点......发送结果时反过来。
根结点和非根节点一个收，一个发，只需要用最简单的`MPI_Send`和`MPI_Recv`即可。
##### 核心MPI代码
```c++
int base = 1;
for (int i = 0; i < steps; ++i) {
if (rank % base == 0) {
int b = base * 2;
if (rank % b == 0) { // recv
MPI_Recv(&recv, 1, MPI_INT, rank + base, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
} else { // send
MPI_Send(&send, 1, MPI_INT, rank - base, 0, MPI_COMM_WORLD);
}
}
send += recv;
base *= 2;
}
base /= 2;
for (int i = 0; i < steps; ++i) {
if (rank % base == 0) {
int b = base * 2;
if (rank % b == 0) { // recv
MPI_Send(&send, 1, MPI_INT, rank + base, 0, MPI_COMM_WORLD);
} else { // send
MPI_Recv(&send, 1, MPI_INT, rank - base, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
}
}
base /= 2;
}
```
*完整代码见附件*
##### 性能结果
| 问题规模 | 任务数 | 总运行时间(sec) |
| -------- | ------------- | --------------- |
| 4 | 4 | 0.022 |
| 8 | 8 | 0.013 |
| 16 | 16 | 0.048 |
| 32 | 32 | 0.156 |
| 64 | 64（跨机器） | 0.394 |
| 128 | 128（跨机器） | 0.466 |
### FOX矩阵乘法
##### 并行算法描述
和其他并行矩阵乘法类似，FOX矩阵乘法基本原理也是矩阵分块乘法，不同的是使用了不同的矩阵块传送策略。
```pascal
输入：n*n矩阵A，n*n矩阵B，初始时A和B分成p个子矩阵，处理器P[i,j]存有块A[i,j]和B[i,j](q * q = p，处理器编号0,0;0,1;...;q-1,q-1)
输出：n*n矩阵C，每个处理器P[i,j]存有块C[i,j]
Begin
for k = 1 to q do
处理器P[i,(i+k) mod q]向所在行广播其A子块
各处理器将接收到的A子块和自己的B子块进行矩阵乘，并加到结果部分和
B子块向上循环移动
end for
End
```
##### MPI实现
MPI实现主要有几个部分：
* 矩阵A和B读取和分发，由进程0完成
* 运算过程
* 广播A子块
* 部分和矩阵乘
* B子块循环移动
* 矩阵C的接收，输出
其中读取，输出，部分矩阵乘虽然是比较基础的操作，但考虑到有广播分发，和接收的操作，需要考虑其数据摆放（layout）。
输入数据矩阵A和B按正常顺序放在文本文件中。
```
4
0 1 2 3
4 5 6 7
8 9 a b
c d e f
```
如果直接按顺序读入内存，则在散播子块时比较不方便，因为每一子块的内存不是连续的，这样要么需要再申请一块内存搬动数据，要么就需要分部分传输每个子块。因而考虑读入时就按子块连续的方式存放。
比如4*4矩阵，分到4个处理器，读入数据摆放如下，这样直接使用一个`MPI_Scatter`就可完成子块分发。
```
0 1 4 5
2 3 6 7
8 9 c d
a b e f
```
在运算前，需要将处理器按行和列分到不同的通信域中，因为处理器初始的编号是一维的，而算法需要二维的编号，并且一些在行内，列内的操作也需要各自的通信域。
运算过程中，广播A子块使用`MPI_Bcast`，在行通信域内广播
B子块的循环移动，可以使用`MPI_Sendrecv`，每个处理器的源是下一行，目标是上一行（环形考虑）
矩阵C的接收和矩阵AB的分发正好相反，使用`MPI_Gather`汇集到0处理器上，这时接收到的C矩阵不是正常顺序，需要根据一定的顺序输出。
##### 核心MPI代码
数据摆放变换，矩阵乘法
```c++
class MatWrap {
private:
float* data_;
int n_;
public:
MatWrap() {}
MatWrap(float* data, int n) : data_(data), n_(n) {}
float* operator[](size_t n) const { return data_ + n * n_; }
float* split_map(int sqrt_q, int i, int j) const {
int sub_n = n_ / sqrt_q;
return data_ + (i / sub_n * sqrt_q + j / sub_n) * sub_n * sub_n +
i % sub_n * sub_n + j % sub_n;
}
friend void MatMultAdd(const MatWrap& a, const MatWrap& b, MatWrap& c);
};
void MatMultAdd(const MatWrap& a, const MatWrap& b, MatWrap& c) {
assert(a.n_ == b.n_);
assert(a.n_ == c.n_);
int n = a.n_;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
for (int k = 0; k < n; ++k) {
c[i][j] += a[i][k] * b[k][j];
}
}
}
}
```
*完整代码见附件*
并行计算部分
```c++
// broadcast sub matrix
MPI_Scatter(mat_a, sub_n * sub_n, MPI_FLOAT, sub_mat_a, sub_n * sub_n,
MPI_FLOAT, 0, MPI_COMM_WORLD);
MPI_Scatter(mat_b, sub_n * sub_n, MPI_FLOAT, sub_mat_b, sub_n * sub_n,
MPI_FLOAT, 0, MPI_COMM_WORLD);
// split comm in col and row
MPI_Comm col_world, row_world;
int col_rank = rank % sqrt_q;
int row_rank = rank / sqrt_q;
MPI_Comm_split(MPI_COMM_WORLD, col_rank, row_rank, &col_world);
MPI_Comm_split(MPI_COMM_WORLD, row_rank, col_rank, &row_world);
// compute
for (int i = 0; i < sqrt_q; ++i) {
// broadcast sub_a
int send_root = (row_rank + i) % sqrt_q;
if (col_rank == (row_rank + i) % sqrt_q) {
memcpy(sub_mat_comm, sub_mat_a, sub_n * sub_n * sizeof(float));
}
MPI_Bcast(sub_mat_comm, sub_n * sub_n, MPI_FLOAT, send_root, row_world);
// calculate sub mat gemm
MatMultAdd(sub_comm, sub_b, sub_c);
// swap sub_b
MPI_Sendrecv_replace(
sub_mat_b, sub_n * sub_n, MPI_FLOAT, (row_rank + sqrt_q - 1) % sqrt_q,
1, (row_rank + 1) % sqrt_q, 1, col_world, MPI_STATUS_IGNORE);
}
// gather result
MPI_Gather(sub_mat_c, sub_n * sub_n, MPI_FLOAT, mat_c, sub_n * sub_n,
MPI_FLOAT, 0, MPI_COMM_WORLD);
```
##### 性能结果
| 问题规模 | 并行数 | 总运行时间(sec) |
| -------- | ------ | --------------- |
| 128 | 1 | 0.123 |
| 256 | 4 | 0.288 |
| 512 | 16 | 0.784 |
| 1024 | 64 | 4.856 |
| 问题规模 | 并行数 | 总运行时间(sec) | 加速比 |
| -------- | ------ | --------------- | ------ |
| 512 | 1 | 2.551 | 1.00 |
| 512 | 4 | 1.099 | 2.32 |
| 512 | 16 | 0.784 | 3.25 |
| 512 | 64 | 1.248 | 2.04 |
### 参数服务器系统
##### 并行算法描述
设系统中总计有N个进程，其中P个进程作为参数服务器进程，而Q个进程作为工作进程（N = P + Q， 且 0 < P << Q）。工作进程和服务器进程的互动过程如下：
1. 第i个工作进程首先产生一个随机数，发送给第i%P个参数服务器进程。然后等待并接收它对应的参数服务器进程发送更新后的数值，之后，再产生随机数，再发送……。
2. 每个参数服务器进程等待并接收来自它对应的所有工作进程的数据，在此之后，经通信，使所有的参数服务器获得所有工作进程发送数据的平均值。
3. 每个参数服务器发送该平均值给它对应的所有工作进程，然后再等待……
```pascal
do while true
工作进程i产生随机数
i mod p参数服务器接收进程i的随机数
所有服务器通信计算所有数据平均值
参数服务器发送平均值给工作进程
end do
```
##### MPI实现
为了方便各种通信，需要划分通信域：
* 所有服务器之间的域：服务器间通信，计算平均值
* 服务器和各自负责的工作进程的域：服务器接收随机数，发送平均值
划分好通信域之后就比较简单了，
1. 服务器从每个工作进程接收随机数
因为最终只需要计算一个平均值，并且要等待每个工作进程，直接使用MPI的归约操作
`MPI_Reduce`，其中的操作`MPI_SUM`，这样每个服务器都有着负责的工作进程的数字的和。
2. 所有服务器通信计算平均值
这里需要每个服务器都获得所有工作进程的平均，而现在每个服务器只有自己负责工作进程的部分和，考虑使用`MPI_Allreduce`，其中的操作仍为`MPI_SUM`，这时每个服务器都有所有工作进程的和，再除以工作进程数即可获得进程的平均数。
3. 服务器发送平均值给工作进程
每个服务器把平均数广播给工作进程即可。
##### 核心MPI代码
```c++
// server or client
int SorC = rank >= SERVER_NUM;
MPI_Comm serverClient, service;
MPI_Comm_split(MPI_COMM_WORLD, SorC, -1, &serverClient);
int service_num = rank % SERVER_NUM;
MPI_Comm_split(MPI_COMM_WORLD, service_num, SorC, &service);
int service_rank;
MPI_Comm_rank(service, &service_rank);
// start
float value;
float sum;
for (int i = 0; i < LOOP; ++i) {
if (SorC) { // client
value = dis(gen);
// printf("c:%d, %d, %f\n", rank, service_rank, value);
} else { // server
value = 0;
}
MPI_Reduce(&value, &sum, 1, MPI_FLOAT, MPI_SUM, 0, service);
if (!SorC) {
// printf("s:%d, %d, %f\n", rank, service_rank, sum);
MPI_Allreduce(&sum, &value, 1, MPI_FLOAT, MPI_SUM, serverClient);
value /= (size - SERVER_NUM);
// printf("s:%d, %d, %f\n", rank, service_rank, value);
}
MPI_Bcast(&value, 1, MPI_FLOAT, 0, service);
// printf("c:%d, %f\n", rank, value);
}
```
*完整代码见附件*
##### 性能结果
| 交换次数 | 服务器数 | 工作进程数 | 总运行时间(sec) |
| -------- | -------- | ---------- | --------------- |
| 1024 | 1 | 8 | 0.036 |
| 1024 | 2 | 16 | 0.087 |
| 1024 | 4 | 24 | 0.268 |
| 1024 | 4 | 32 | 0.292 |
| 交换次数 | 服务器数 | 工作进程数 | 总运行时间(sec) | 加速比 |
| -------- | -------- | ---------- | --------------- | ------ |
| 1024 | 1 | 64 | 0.441 | - |
| 1024 | 2 | 64 | 1.496 | - |
| 1024 | 4 | 64 | 0.714 | - |
| 1024 | 8 | 64 | 0.684 | - |
从结果上看，当工作进程数增加时，服务器成比例增加也难以保证总运行时间，因为每次运行需要所有服务器间通信一次。
而当工作进程数固定下，增加服务器数并不能使服务时间线性下降，服务器为1时是一个特例，因为服务器间不再需要通信。因而在服务器负载均衡时，尽量减少服务器间的通信，有时不均衡的微服务设计会好于非常均匀的负载均衡。
## 个人实验
##### 问题描述
###### LSTM长短时记忆
长短期记忆（英语：Long Short-Term Memory，LSTM）是一种时间递归神经网络（RNN），论文首次发表于1997年。由于独特的设计结构，LSTM适合于处理和预测时间序列中间隔和延迟非常长的重要事件。LSTM可以应用在语音识别，机器翻译，自然语言处理，手写识别等机器学习应用。
##### 串行算法描述
###### 计算公式
$$
i_t=\text{sigmoid}(W_{ii}x_t+b_{ii}+W_{hi}h_{(t-1)}+b_{hi})\\
f_t=\text{sigmoid}(W_{if}x_t+b_{if}+W_{hf}h_{(t-1)}+b_{hf})\\
g_t=\text{tanh}(W_{ig}x_t+b_{ig}+W_{hc}h_{(t-1)}+b_{hg})\\
o_t=\text{sigmoid}(W_{io}x_t+b_{io}+W_{ho}h_{(t-1)}+b_{ho})\\
c_t=f_t\cdot c_{(t-1)}+i_t\cdot g_t\\
h_t=o_t\cdot \text{tanh}(c_t)
$$
对于上面公式的矩阵乘加计算，使用增广矩阵形式，化为矩阵乘。
$$
i_t=\text{sigmoid}(W_{ii}x_t+W_{hi}h_{(t-1)})\\
f_t=\text{sigmoid}(W_{if}x_t+W_{hf}h_{(t-1)})\\
g_t=\text{tanh}(W_{ig}x_t+W_{hc}h_{(t-1)})\\
o_t=\text{sigmoid}(W_{io}x_t+W_{ho}h_{(t-1)})\\
c_t=f_t\cdot c_{(t-1)}+i_t\cdot g_t\\
h_t=o_t\cdot \text{tanh}(c_t)
$$
###### 串行伪代码
```pascal
输入：增广输入向量x[1..T]，增广矩阵形式权值W_ii,W_if,W_ig,W_io
输出：预测向量h[1..T]
for t = 1 to T
i[t] = sigmoid(W_ii * x[t] + W_hi * h[t-1])
f[t] = sigmoid(W_if * x[t] + W_hf * h[t-1])
g[t] = tanh(W_ig * x[t] + W_hc * h[t-1])
o[t] = sigmoid(W_io * x[t] + W_ho * h[t-1])
c[t] = f[t] .* c[t-1] + i[t] .* g[t]
h[t] = o[t] .* tanh(c[t])
end for
```
##### 并行算法描述
###### PCAM设计
1）划分
从功能上看，上面算法中主要可分为`it,ft,gt,ot,ct,ht`
从任务负载来看，上面算法中运行时间比重较大的是矩阵乘加运算——`it,ft,gt`和`ot`的工作量相同，而`ct,ht`的工作量稍少。
在每个矩阵乘加运算中，可以按照分块矩阵乘法对矩阵乘加进一步划分。
2）通信
计算`it,ft,gt,ot`需要`xt,ht-1`，而`ct`需要`ft,it,gt`，`ht`需要`ot,ct`
3）组合
考虑`ct,ht`工作量较少，且需要传输较多数据，将`ct,ht`和`it`合并
4）映射
`it,ct,ht; ft; gt; ot`分别放在4个处理机上，这时这4个处理机任务量非常接近，且此时传输数据量较少。
而`it,ct,ht; ft; gt; ot`各自计算矩阵乘时，再将矩阵分块，分配到各自通信域中的处理机上。
###### 依赖关系分析
从上面公式和串行伪代码来看，依赖关系非常明显，输出$c_t,h_t$流依赖于$c_{(t-1)},f_t,g_t,i_t,o_t$，而$i_t,f_t,g_t,o_t$又六依赖于$h_{(t-1)}$，对于这种随时间的迭代计算，不同时间$t$之间不能并行计算，因而考虑$i_t,f_t,g_t,o_t$可以并行计算，而在$i_t,f_t,g_t,o_t$内有矩阵乘加计算，也可以使用分块矩阵的并行计算。
###### MPI+OpenMP设计
对于LSTM的各个时间点的计算，无法并行
对于LSTM的一个时间点内的计算，考虑在同一个节点上使用MPI分为$i_t,f_t,g_t,o_t$4个部分在4个节点上计算，需要广播$x_t,h_{(t-1)}$，最后由0号节点收集结果，再在0号节点上计算$c_t,h_t$
$i_t,f_t,g_t,o_t$每个部分内的矩阵乘法使用OpenMP计算
##### MPI+OpenMP实现
###### 矩阵乘，矩阵加，激活，OpenMP实现
```c++
void MatrixMult(float *a, float *b, float *c, int n, bool init = false) {
#pragma omp parallel for num_threads(OMP_THREADS)
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
if (init) {
c[i * n + j] = 0.0f;
}
for (int k = 0; k < n; ++k) {
c[i * n + j] += a[i * n + k] * b[k * n + j];
}
}
}
}
void MatrixDot(float *a, float *b, float *c, int n) {
#pragma omp parallel for num_threads(OMP_THREADS)
for (int i = 0; i < n * n; ++i) {
c[i] = a[i] * b[i];
}
}
void MatrixAdd(float *a, float *b, float *c, int n) {
#pragma omp parallel for num_threads(OMP_THREADS)
for (int i = 0; i < n * n; ++i) {
c[i] = a[i] + b[i];
}
}
void tanh(float *x, float *y, int n) {
#pragma omp parallel for num_threads(OMP_THREADS)
for (int i = 0; i < n; ++i) {
y[i] = std::tanh(x[i]);
}
}
void sigmoid(float *x, float *y, int n) {
#pragma omp parallel for num_threads(OMP_THREADS)
for (int i = 0; i < n; ++i) {
y[i] = 1.0f / (1.0f + std::exp(-x[i]));
}
}
```
###### 随时间迭代计算，MPI实现
```c++
// computing
for (int i = 0; i < seq_len; ++i) {
// broadcast x
if (rank == 0) {
memcpy(sub_x, x + i * n * n, n * n * sizeof(float));
}
MPI_Bcast(sub_x, n * n, MPI_FLOAT, 0, MPI_COMM_WORLD);
// broadcast h
if (i == 0) {
for (int i = 0; i < n * n; ++i) {
sub_h[i] = 0.0f;
}
} else {
MPI_Bcast(sub_h, n * n, MPI_FLOAT, 0, MPI_COMM_WORLD);
}
// matrix multply
MatrixMult(sub_w, sub_x, ifgo, n, true);
MatrixMult(sub_w + n * n, sub_h, ifgo, n);
// active
if (rank == 2) {
tanh(ifgo, ifgo, n * n);
} else {
sigmoid(ifgo, ifgo, n * n);
}
// gather result
if (rank == 0) {
MPI_Recv(ft, n * n, MPI_FLOAT, 1, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
MPI_Recv(gt, n * n, MPI_FLOAT, 2, 2, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
MPI_Recv(ot, n * n, MPI_FLOAT, 3, 3, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
} else if (rank == 1) {
MPI_Send(ifgo, n * n, MPI_FLOAT, 0, 1, MPI_COMM_WORLD);
} else if (rank == 2) {
MPI_Send(ifgo, n * n, MPI_FLOAT, 0, 2, MPI_COMM_WORLD);
} else if (rank == 3) {
MPI_Send(ifgo, n * n, MPI_FLOAT, 0, 3, MPI_COMM_WORLD);
}
// compute ct ht
if (rank == 0) {
MatrixDot(ft, c, c, n);
MatrixDot(ifgo, gt, gt, n);
MatrixAdd(c, gt, c, n);
tanh(c, sub_h, n * n);
MatrixDot(ot, sub_h, sub_h, n);
memcpy(h + i * n * n, sub_h, n * n * sizeof(float));
}
}
```
*完整代码见附件*
##### 性能结果
| 问题规模 | OpenMP线程数 | 总运行时间(sec) |
| -------- | ------------ | --------------- |
| 64 | 1 | 0.492 |
| 128 | 2 | 1.725 |
| 256 | 4 | 6.681 |
| 512 | 8 | 46.498 |
| 问题规模 | OpenMP线程数 | 总运行时间(sec) | 加速比 |
| -------- | ------------ | --------------- | ------ |
| 256 | 1 | 19.848 | 1.00 |
| 256 | 2 | 11.068 | 1.79 |
| 256 | 4 | 7.162 | 2.77 |
| 256 | 8 | 6.252 | 3.17 |
## 分组实验
| 章 | 15 | 19 | 22 |
| ------------ | ------- | ----- | ---- |
| 算法程序数 | 4 | 2 | 2 |
| 分到的程序号 | 1 | 1 | 1 |
| 分到的程序 | closure | gauss | fft |
### closure-MPI
##### 性能结果
1
| 问题规模 | MPI进程数 | 总运行时间(sec) |
| -------- | --------- | --------------- |
| 64 | 1 | 0.122 |
| 128 | 4 | 0.267 |
| 256 | 16 | 0.276 |
| 512 | 64 | 4.923 |
2
| 问题规模 | MPI进程数 | 总运行时间(sec) | 加速比 |
| -------- | --------- | --------------- | --------------- |
| 256 | 1 | 1.014 | 1.00 |
| 256 | 4 | 0.811 | 1.25 |
| 256 | 16 | 1.396 | 0.72 |
| 256 | 64 | 4.923 | 0.21 |
##### 改进
* 优化了数据传输，使用`MPI_Bcast, MPI_Scatter, MPI_Gather`替换了`for-MPI_Send-MPI_Recv`
* 使用管道重定向输入输出
* 自动化性能测试脚本
*完整代码见附件*
### closure-Hybrid-omp-mpi
##### MPI+OpenMP实现
原有程序为MPI行划分实现，很难再在行内新增OpenMP并行。因而考虑重新设计并行方案。closure运算的本质是布尔矩阵乘，因而考虑使用分块并行矩阵乘法，使用MPI分块并行矩阵乘法，在每个子矩阵块乘法再使用OpenMP行划分矩阵并行乘法。
###### 核心代码
```c++
void MatMultAdd(const MatWrap& a, const MatWrap& b, MatWrap& c) {
assert(a.n_ == b.n_);
assert(a.n_ == c.n_);
int n = a.n_;
#pragma omp parallel for num_threads(OMP_THREADS)
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
for (int k = 0; k < n; ++k) {
c[i][j] = static_cast(c[i][j] + a[i][k] * b[k][j]);
}
}
}
}
```
```c++
for (int k = 0; k <= std::log2f(n); ++k) {
for (int i = 0; i < sub_n; ++i) {
for (int j = 0; j < sub_n; ++j) {
sub_c[i][j] = 0;
}
}
// broadcast sub matrix
MPI_Scatter(mat_a, sub_n * sub_n, MPI_INT, sub_mat_a, sub_n * sub_n,
MPI_INT, 0, MPI_COMM_WORLD);
MPI_Scatter(mat_a, sub_n * sub_n, MPI_INT, sub_mat_b, sub_n * sub_n,
MPI_INT, 0, MPI_COMM_WORLD);
// split comm in col and row
MPI_Comm col_world, row_world;
int col_rank = rank % sqrt_q;
int row_rank = rank / sqrt_q;
MPI_Comm_split(MPI_COMM_WORLD, col_rank, row_rank, &col_world);
MPI_Comm_split(MPI_COMM_WORLD, row_rank, col_rank, &row_world);
// compute
for (int i = 0; i < sqrt_q; ++i) {
// broadcast sub_a
int send_root = (row_rank + i) % sqrt_q;
if (col_rank == (row_rank + i) % sqrt_q) {
memcpy(sub_mat_comm, sub_mat_a, sub_n * sub_n * sizeof(int));
}
MPI_Bcast(sub_mat_comm, sub_n * sub_n, MPI_INT, send_root, row_world);
// calculate sub mat gemm
MatMultAdd(sub_comm, sub_b, sub_c);
// swap sub_b
MPI_Sendrecv_replace(
sub_mat_b, sub_n * sub_n, MPI_INT, (row_rank + sqrt_q - 1) % sqrt_q,
1, (row_rank + 1) % sqrt_q, 1, col_world, MPI_STATUS_IGNORE);
}
// gather result
MPI_Gather(sub_mat_c, sub_n * sub_n, MPI_INT, mat_a, sub_n * sub_n, MPI_INT,
0, MPI_COMM_WORLD);
// print result
if (rank == 0) {
MatWrap mc(mat_a, n);
printf("loop:%d\n", k);
mc.print(true, sqrt_q);
}
}
```
*完整代码见附件*
##### 性能结果
| 问题规模 | OpenMP线程数 | MPI进程数 | 总运行时间(sec) | 加速比 |
| -------- | ------------ | --------- | --------------- | ------ |
| 64 | 1 | 1 | 0.122 | 1.00 |
| 64 | 2 | 1 | 0.090 | 1.36 |
| 64 | 4 | 1 | 0.067 | 1.82 |
| 256 | 1 | 4 | 1.255 | 1.00 |
| 256 | 2 | 4 | 0.834 | 1.50 |
| 256 | 4 | 4 | 0.647 | 1.94 |
### gauss-MPI
##### 性能结果
1
| 问题规模 | 任务数 | 总运行时间 | 分发数据时间 | 并行计算时间 |
| -------- | ------ | ---------- | ------------ | ------------ |
| 512 | 1 | 0.387903 | 0.112848 | 0.275056 |
| 1024 | 2 | 1.265614 | 0.335754 | 0.929859 |
| 2048 | 4 | 4.077509 | 1.161536 | 2.915974 |
| 4096 | 8 | 19.073884 | 4.172863 | 14.901021 |
| 8192 | 16 | 100.614923 | 16.937205 | 83.677718 |
2
| 问题规模 | 任务数 | 总运行时间 | 分发数据时间 | 并行计算时间 | 加速比 |
| -------- | ------ | ---------- | ------------ | ------------ | ------ |
| 2048 | 1 | 5.461998 | 1.217323 | 4.244675 | 1.00 |
| 2048 | 2 | 4.983760 | 1.240886 | 3.742874 | 1.13 |
| 2048 | 4 | 4.045694 | 1.148780 | 2.896914 | 1.47 |
| 2048 | 8 | 4.017666 | 1.161650 | 2.856016 | 1.49 |
| 2048 | 16 | 4.135822 | 1.286142 | 2.849679 | 1.48 |
##### 改进
* 使用Python，numpy生成输入数据
* 重定向输出，使用管道重定向输出
* 自动化性能测试脚本
*完整代码见附件*
### fft-MPI
##### 性能结果
1
| 问题规模 | MPI进程数 | 总运行时间(sec) |
| -------- | --------- | --------------- |
| 512 | 1 | 0.036836 |
| 1024 | 2 | 0.073674 |
| 2048 | 4 | 0.148747 |
| 4096 | 8 | 0.347353 |
2
| 问题规模 | MPI进程数 | 总运行时间(sec) | 加速比 |
| -------- | --------- | --------------- | ------ |
| 4096 | 1 | 2.323001 | 1.00 |
| 4096 | 2 | 1.162240 | 2.00 |
| 4096 | 4 | 0.605259 | 3.84 |
| 4096 | 8 | 0.349958 | 6.64 |
##### 改进
* 使用Python，numpy生成输入数据
* 重定向输出，使用管道重定向输出
* 自动化性能测试脚本
* 修改了输入数量规模上限
*完整代码见附件*