public

report.md

@@ -0,0 +1,358 @@
+# 并行程序设计实验报告
+
+<p align=right>陈翊辉</p>
+<p align=right>SA19011116</p>
+[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<void*>(&send), 1, MPI_INT, target, target, MPI_COMM_WORLD);
+      MPI_Recv(reinterpret_cast<void*>(&recv), 1, MPI_INT, target, rank, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    } else {
+      MPI_Recv(reinterpret_cast<void*>(&recv), 1, MPI_INT, target, rank, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Send(reinterpret_cast<void*>(&send), 1, MPI_INT, target, target, MPI_COMM_WORLD);
+    }
+    send += recv;
+    base *= 2;
+  }
+```
+
+*完整代码见附件*
+
+##### 性能结果
+
+
+
+#### 二叉树求和
+
+##### 并行算法描述
+
+二叉树求和使$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;
+  }
+```
+
+*完整代码见附件*
+
+##### 性能结果
+
+
+
+### 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);
+```
+
+##### 性能结果
+
+### 参数服务器系统
+
+##### 并行算法描述
+
+设系统中总计有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);
+  }
+```
+
+*完整代码见附件*
+
+##### 性能结果
+
+
+## 个人实验
+
+##### 并行算法描述
+
+
+##### MPI实现
+
+##### 核心MPI代码
+
+##### 性能结果
+
+## 分组实验
+
+| 章           | 15      | 19    | 22   |
+| ------------ | ------- | ----- | ---- |
+| 算法程序数   | 4       | 2     | 2    |
+| 分到的程序号 | 1       | 1     | 1    |
+| 分到的程序   | closure | gauss | fft  |
+