对一并行程序实验的简单理解

写在前面

使用一维数组

在写并行程序时，会遇到很多针对矩阵、张量的计算，尽管它们都是以多维数组的方式组合在一起，但一般情况下，我更加倾向于使用一维数组来表示这些多维数组的数据。我认为，这样做的好处有：

1. 方便函数参数类型的统一。如果程序中既声明了二维数组（这意味着它是一个 pointer to pointer 的类型），又存在一维数组，那么在调用函数时（比如需要生成随机数据，或者检查运行结果是否正确），会出现类型不一致，就需要实现两个功能类似的函数，增加了代码冗余。当然了，这里我们只讨论 C 语言的情况 :sweat:。
2. 方便索引的计算。计算偏移量是并行编程中常做的事情，然而在多维数组中，指针偏移是很容易出错且吃力的事情（很多时候单看 p+1 你甚至不知道指针到底偏移了多少！）。而如果全部使用一维数组，偏移的量就很容易理解，你再也不会被多维指针搞得晕头转向:smile:。
3. 更接近机器底层。对计算机系统稍做了解，就会发现所谓的“多维数组”其实是程序员抽象出来的东西，计算机真正认识的只有一维数组。使用一维数组写并行程序，不仅对 MPI 更加友好（大部分 MPI 函数只接受一维指针），还会更有助于我们理解机器行为，挖掘并行潜力。

使用一维数组虽然有如上好处，但也有明显的弊端：索引表达不方便。当然也有办法可以解决这些问题：

1. 使用宏函数。注意一定要给每个运算部分加上括号。

   #define INDEX(i, j) (((i)*(N))+(j))

2. 如果你对宏函数嗤之以鼻，那么可以使用 inline 函数代替宏函数，性能不会降低且更加安全。

使用堆空间并注意内存泄露

并行编程往往涉及非常大的数据量，因此推荐将数据放入堆空间中。不过，这会带来一个非常棘手的问题：内存泄露。C 语言中没有自动的内存回收工具，内存能不能被回收利用全靠程序员自觉，但良好的编程习惯和适当的技巧可以尽可能地避免这类问题：

do {
  // 分配堆空间
  pa = malloc(N);
    // 若失败 break 退出
    if (pa==NULL)
        break;
}while(0);
if (pa!=NULL) {
    free(pa);
	pa = NULL;
}

用循环 do ... while(0) 将程序主体包裹起来，这样一旦分配内存出错就可以直接 break 跳出来，然后将其释放掉。还避免使用了 goto 语句。另外需要注意的是，释放内存前一定要检查指针是否已经指向 NULL ，这可避免二次释放造成的运行时错误。

使用串行程序进行验算

使用串行程序检验最终计算结果对并行编程的重要性不言而喻。但在调试时，串行程序也可以贡献出自己的一份力量。若遇到复杂的程序设计，可能需要分步完成程序，每一步都用串行程序验证正确性。稳扎稳打，步步为营才是王道。

题目描述

矩阵 $A$ 和 $B$ 均为 $N \times N$ 的双精度矩阵，现有 $P$ 个处理器，针对以下程序片段，分别采用按行块连续划分和棋盘划分的方式，给出相应的 MPI 程序实现。

for (i = 1; i < N-1; i++)
    for(j = 1; j < N-1; j++)
        B[i][j] = (A[i-1][j]+A[i][j-1]+A[i+1][j]+A[i][j-1]) / 4.0;

按行连续划分

方法概要

按行连续划分，即将矩阵以行为单位分割，并交给不同的处理器进行处理。下图就是我们对矩阵 $A$ 行连续划分的示意图。除主处理器外，每个处理器轮流地取出连续的三行。在程序一开始，假设只有0号处理器拥有所有的数据。

我们的程序需要完成以下步骤：

1. 在某个处理器中生成随机数据。为了计算的方便，我假设编号最大的处理器中握有数据，并且还负责最后结果的汇总。
2. 将生成的随机数据分发给不同的处理器。从上面的程序片段中看出，计算一个矩阵 $B$ 元素需要三行的矩阵 $A$ 元素，因此至少需要给一个处理器三行数据。
3. 各个处理器完成自己的任务。
4. 将计算好的结果从各个处理器发出，由主处理器收集整理。

细节处理

细节决定成败，这对并行编程来说尤为如此。特别地，需要注意

• 从处理器会接到多个来自主处理器的数据，必须要用明确的标签区分这些数据。
• 处理好余数情况，某几个从处理器收到的数据可能会比其他从处理器多一个，需要一个私有变量记录处理器需要计算多少行矩阵 $B$。当然，也可以采用对矩阵 $A$ 补 0 对齐的方法
• 收集计算结果时，要将起点定在矩阵 $B$ 的（1，1）中，还要弄清楚到底哪一行计算结果存放在哪一个处理器中。
• 时刻注意索引的计算

部分代码解析

先是程序的初始化和随机生成数据，并串行地执行计算。

double *A, *B, *B2;
// malloc A, B, B2 with N*N*sizeof(double)
// init MPI
int id_procs, num_procs, num_1;
MPI_Status status;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
MPI_Comm_rank(MPI_COMM_WORLD, &id_procs);

num_1 = num_procs -1;
// Proc#N-1 randomize the data, and compute B2 for verification
if (id_procs == num_1) {
  random_array(A, N*N);
  comp(A, B2, N*N);
}
// wait Proc#N-1 for computing
MPI_Barrier(MPI_COMM_WORLD);

随后，开始真正的并行计算：先让拥有数据的线程将数据广播到相应的线程中。注意，根据图上的示例，需要发送连续的三层数据。需要注意到，N 不一定能整除 (num_procs-1)，需要考虑留下的余数。ctn 变量记录收到数据的份数，表明了后面的计算需要循环的次数。

// ********************************************* //
// Parallel Computing Part
// Proc#N-1 broadcast 3 lines of A to each Proc
int ctn = 0;
for(int i = 0; i < N-2; i++) {
  if (id_procs == num_1) {
    int dest = i % num_1;
    int tag = i / num_1;
    MPI_Send(&A[INDEX(i, 0)], N*3, MPI_DOUBLE, dest, tag, MPI_COMM_WORLD);
  }
}

for(int i = 0; i < (N-2)/num_1; i++) {
  if (id_procs != num_1) {
    MPI_Recv(&A[INDEX(3*ctn, 0)], 3*N, MPI_DOUBLE, num_1, ctn, MPI_COMM_WORLD, &status);
    ctn++;
  }
}
if (id_procs < (N-2) % num_1) {
  MPI_Recv(&A[INDEX(ctn*3, 0)], 3*N, MPI_DOUBLE, num_1, ctn, MPI_COMM_WORLD, &status);
  ctn++;
}

分发数据完成！接下来就是每个线程的计算部分。这里的实现和串行程序十分类似。注意，ctn 变量记录了该线程受到了几组数据，因此外层循环需要循环 ctn 次：

// compute
if (id_procs != num_1) {
  for(int i = 1; i <= ctn; i++) {
    for(int j = 1; j < N-1; j++) {
      B[INDEX(i, j)] = (A[INDEX(i-1, j)]+A[INDEX(i, j+1)]+A[INDEX(i+1, j)]+A[INDEX(i, j-1)]) / 4.0;
    }
  }
}

把计算得到的矩阵$B$全都发给 Proc#num-1。

// Gather
for(int i = 0; i < N-2; i++) {
    if (id_procs == num_1) {
        int src = i % num_1;
        MPI_Recv(&B[INDEX(i+1, 1)], N-2, MPI_DOUBLE, src, i/num_1+N, MPI_COMM_WORLD, &status);
    }
    else {
        for(int j = 0; j < ctn; j++)
            MPI_Send(&B[INDEX(j+1, 1)], N-2, MPI_DOUBLE, num_1, j+N, MPI_COMM_WORLD);
    }
}

棋盘划分

方法概要

先以块棋盘划分来实现并行。在程序片段中，矩阵 $B$ 的有效结果范围为 [1~N-2, 1~N-2] 。若将矩阵 B 按如下示意图划分为棋盘状，每个小格子由一个处理器单独完成（以下就用“格子”数据类型称呼它们）。

此时，程序需要完成以下步骤：

1. 在某个处理器中生成随机数据。为方便，这次选择编号0的处理器生成随机数据。
2. 将生成的随机数据分发给不同的处理器。按照总共有的处理器数量，合理分配每行处理器的数量以及总行数，分配完成后，将随机数据发给不同的处理器。遇到余数的问题，用补 0 对齐的方法解决（虚线表示）。
3. 各个处理器计算自己“格子”内的数据
4. 将计算好的结果从各个处理器发出，由编号0的处理器收集整理。

细节处理

除了之前在行连续划分中提到的外，还要特别注意到：棋盘划分出的“格子”在一维数组中是不连续的。把一个格子发送给某个处理器，如果不采用一些小技巧，是很难做到的（你需要一个循环，并计算格子中每一行的索引），但如果使用 MPI 提供的创建数据类型功能，就很容易办到。通过下面的三条语句，就可以轻松地创建出类型为“格子”的数据类型。

MPI_Datatype SubMat;
MPI_Type_vector(a+2, b+2, N, MPI_DOUBLE, &SubMat);
MPI_Type_commit(&SubMat);

对于矩阵 $A$ 来说，一个处理器所需要的数据，即所谓的棋盘格子，事实上就是一个有 a+2 个数据段，每个段长度为 b+2 ，段与段之间相隔 N 个单位的“结构体”，使用 MPI_Type_vector 函数就可以轻松描述它，如果你想了解更多，请参考 MPICH 官网。

部分代码解析

double *A, *B, *B2;
// init MPI as usual
int id_procs, num_procs;		MPI_Status status;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
MPI_Comm_rank(MPI_COMM_WORLD, &id_procs);

创建新的数据类型，用于描述前面提到的“格子”。

int rows = sqrt(num_procs);
int cols = num_procs / rows;
int a = (N-2 + rows-1) / rows;
int b = (N-2 + cols-1) / cols;
int alloc_num = (a+1)*(b+1)*num_procs;

/* Proc#0 randomize the data */
/* ..... */
// Proc#0 broadcast (a+2)x(b+2) mat
MPI_Datatype SubMat;
MPI_Type_vector(a+2, b+2, N, MPI_DOUBLE, &SubMat);
MPI_Type_commit(&SubMat);

有了这个数据类型，就不需要纠结于繁琐的边界问题，按照索引计算规则，直接将矩阵$A$ 切割好发送出去。当然索引计算什么的仍然很不友好 :sweat:，不过接收方的代码确实简单不少。

if (id_procs == 0) {
  for(int i = 0; i < rows; i++) {
    for(int j = 0; j < cols; j++) {
      if (i == 0 && j == 0)
        continue;
      MPI_Send(A+i*a*N+b*j, 1, SubMat, j+cols*i, 0, MPI_COMM_WORLD);
    }
  }
} else {
  MPI_Recv(A, 1, SubMat, 0, 0, MPI_COMM_WORLD, &status);
}

开始计算，这边的计算完全和串行程序一致，因为我们只是将矩阵$A$ 和 $B$ 变小了而已。

// compute same as the single version
comp(A, B, a, b);

收集计算好的数据。注意到，此时数据类型变化了，因为矩阵$A$ 的格子比$B$的要大一圈，用一个数据类型就无法描述了。

// Gather result
MPI_Datatype SubMat_B;
MPI_Type_vector(a, b, N, MPI_DOUBLE, &SubMat_B);
MPI_Type_commit(&SubMat_B);
if (id_procs == 0) {
  for(int i = 0; i < rows; i++) {
    for(int j = 0; j < cols; j++) {
      if (i == 0 && j == 0)
        continue;
      MPI_Recv(&B[INDEX(a*i+1, b*j+1)], 1, SubMat_B, i*cols+j, 1, MPI_COMM_WORLD, &status);
    }
  }
} else {
  int x = id_procs / cols;
  int y = id_procs % cols;
  MPI_Send(&B[INDEX(1, 1)], 1, SubMat_B, 0, 1, MPI_COMM_WORLD);
}

当然，还有很多不同的做法：可以使用 MPI_Pack 函数对“格子”类型的数据进行打包；也可以用面向对象的思想，直接封装一个 class ；当然也可以用循环这种最笨拙的办法传递“格子”

我的代码实现见 github 仓库。