并行程序设计导论 第三章习题

3.1
在问候程序中，如果strlen(greeting)代替strlen(greeting)+1来计算进程1、2、…、comm_sz-1发送消息的长度，会发生什么情况？如果用MAX_STRING代替strlen(greeting)+1又会是什么结果？你可以解释这些结果吗？

解答：
①这样的替换并不影响显示，不过传递的是两个不大一样的字符串。
【引用 58页】在我们的程序，参数msg_size是消息字符串加上C语言中字符串结束符\0所占的字符数量。参数msg_type的值是MPI_CHAR。这两个参数一起告知系统整个消息有strlen(geeting)+1个字符。
②这个问题与上面的问题完全相反。这里MAX_STRING是绝对大于strlen(greeting)+1的，所以传递过去的字符数量会是MAX_STRING，会将接收区的缓存占满。不过，对于打印输出来说，并没有什么影响。

3.2
改变梯形积分法，使其能够在comm_sz无法被n整除的情况下，正确估算积分值(假设n≥comm_sz)

解答：

#include <stdio.h>
#include <mpi.h>

double f(double x){
  return x*x;
}

double Trap(
  double left_endpt,   /*  in  */
  double right_endpt,  /*  in  */
  int trap_count,      /*  in  */
  double base_len      /*  in  */
  ){
  double estimate, x;
  int i;

  estimate = (f(left_endpt) + f(right_endpt))/2.0;
  for (i = 0; i < trap_count - 1; i++){
    x = left_endpt + (i + 1) * base_len;
    estimate += f(x);
  }
  return estimate * base_len;
}

int main(){
  int my_rank, comm_sz, n = 64, local_n;
  double a = 0.0, b = 3.0, h, local_a, local_b, remain = 0.0;
  double local_int, total_int;
  int source;

  MPI_Init(NULL, NULL);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

  //Get_data(my_rank, comm_sz, &a, &b, &n);

  h = (b - a)/n; //h is the same for all processes

  local_n = n/comm_sz; // So is the number of trapezoids

  local_a = a + my_rank * local_n * h;
  local_b = local_a + local_n * h;
  local_int = Trap(local_a, local_b, local_n, h);



  if (my_rank != 0){
    MPI_Send(&local_int, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
  } else {

    total_int =local_int;
    for (source = 1; source < (comm_sz); source++){
      MPI_Recv(&local_int, 1, MPI_DOUBLE, source, 0, MPI_COMM_WORLD,
               MPI_STATUS_IGNORE);
      total_int += local_int;
    }

    // 这里是因为有一些矩阵没有计算到，所以，非整除的情况下结果会有比较大的缺失。
    // 这里将这些矩形在0号线程上补上。
    if (n%comm_sz){
      remain = Trap(h * local_n * comm_sz, b, n%comm_sz, h);
    }
    
    total_int += remain;
  }

  if (my_rank == 0){
    printf("With n = %d trapezoids, our estimate\n", n);
    printf("of the integral from %f to %f = %.15e\n", a, b, total_int);
  }

  MPI_Finalize();
  return 0;
}

3.3
梯形积分法程序中那些变量是局部的，那些是全局的？

解答：
【引用63页】注意，我们对标识符的选择，是为了区分局部变量与全局变量。局部变量只在使用它们的进程内有效。梯形积分法程序中的例子有：local_a,local_b和local_n。如果变量在所有进程中都有效，那么该变量就称为全局变量。该程序中的例子有：变量a,b和n。这与你在编程导论课上学到的用法不同。在编程导论课上，局部变量是指单个函数的私有变量，而全局变量是指所有函数都可以访问的变量。

3.4
mpi_output.c程序中，每个进程只打印一行输出。修改程序，使程序能够按进程号的顺序打印，即，进程0先输出，然后进程1，一次类推。

解答：
其实这里的实现和3.1中的问候程序差不多。

#include <stdio.h>
#include <string.h>
#include <mpi.h>

int main(){
  int my_rank, comm_sz, q;
  char message[100];

  MPI_Init(NULL, NULL);
  MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

  //printf("Proc %d of %d > Does anyone have a toothpick?\n", my_rank, comm_sz);

  if (my_rank != 0){
    sprintf(message, "Proc %d of %d > Does anyone have a toothpick?", my_rank, comm_sz);
    MPI_Send(message, strlen(message)+1, MPI_CHAR, 0, 0, MPI_COMM_WORLD);
  } else {
    printf("Proc %d of %d > Does anyone have a toothpick??\n", my_rank, comm_sz);
    for(q = 1; q < comm_sz; q++){
      MPI_Recv(message, 100, MPI_CHAR, q, 0,MPI_COMM_WORLD, MPI_STATUS_IGNORE);
      printf("%s\n", message);
    }
  }

  MPI_Finalize();
  return 0;
}

3.5
二叉树中有着从根到每个节点的最短路径，这条路径的长度称为该节点的深度。每个非叶子节点结点都有2个子节点深度都相同的满二叉树称为完全二叉树。用数学推导证明：如果T是一颗有着n个叶子节点的完全二叉树，那么叶子节点的深度为log(n)【底为2】

解答：
当n=1 => 0=log(1)时显然。
假设当T具有n个叶子结点时的完全二叉树的深度为log(2)，则
当n=2^k（以及2(k+1),…,2^{(k+p)）时，由归纳假设知前2}k个结点构成深度为log(n)的树，再由完全二叉树的定义知剩余的1（或2,…,2^k)个结点均填在第log(n)-1层上，故深度刚好增加了1。
故n<=2^(k+1)时命题成立。证毕。
（首先最好能先从直观上理解：完全二叉树中：
1层 有1个叶子结点；
2层 有2个叶子结点；
3层 有4个叶子结点；
……
k层 有2^(k-1)个叶子结点；)

3.6
假设comm_sz=4，x是一个拥有n=14个元素的向量
a 如何用块划分法在一个程序的进程间分发x的元素
b 如何用循环划分法在一个程序的进程间分发x的元素
c 如何用b=2大小的块用块-循环划分法在一个程序的进程间分发x的值
你的分配方法应该通用性足够强，无论comm_sz和n取何值都能分配。你应该使你的分配“公正”，如果q和r是任意的2个进程，分给q和r的x分量个数的差应尽可能小。

解答：
a 块划分

进程0  | 0  1  2  3
进程1  | 4  5  6  7
进程2  | 8  9  10 11
进程3  | 12 13

b 循环划分

进程0  | 0  4  8  12
进程1  | 1  5  9  13
进程2  | 2  6  10
进程3  | 3  7  11

c 块-循环划分

进程0  | 0  1  8  9  
进程1  | 2  3  10 11  
进程2  | 4  5  12 13
进程3  | 6  7

3.7
如果通信子只包含一个进程，不同的MPI集合通信函数分别会做什么。

解答：
这个验证起来很简单，就是将书中的程序，使用mpicc编译后，直接运行二进制文件就行了。
可以看到的是，所有的操作就相当于这个进程发送信息给自己，然后完成对应的操作。
个人感觉，当只有一个进程存在在通信子中时，就没有必要使用MPI来进行通讯了。

3.8
假定comm_sz=8, n=16
a 画一张图来说明进程0要分发n个元素的数组，怎样使用拥有comm_sz个进程的树形结构的通信来实现MPI_Scatter。
b 画一张图来说明原先分布在comm_sz个进程间的n个数组元素由进程0保存，怎样使用树形结构的通信来实现MPI_Gather.

解答：
这里个人感觉在细看一下73和74页的内容，就能知道这两个函数，是如何进行分配的了。
至于这里说的树形结构，个人没弄明白是怎么样的一个图。

3.9
编写一个MPI程序实现向量与标量相乘以及向量点积的功能。用户需要输入2个向量和1个标量，都有进程0读入并分配给其他进程，计算结果由进程0计算和保存，并最终由进程0打印出来。假定向量的秩n可以被comm_sz整除。

解答：

#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>

#define N 4

void Read_vector(
  double vector_a[],
  double vector_b[],
  double *scalar,
  int local_n,
  int n,
  int my_rank,
  MPI_Comm comm){

  double va[N], vb[N];
  int i;

  if (my_rank == 0){
    printf("Enter the vector_a:\n");
    for (i = 0; i < n; i++){
      scanf("%lf", &va[i]);
    }

    printf("Enter the vector_b:\n");
    for (i = 0; i < n; i++){
      scanf("%lf", &vb[i]);
    }

    printf("Enter the scalar:\n");
    scanf("%lf", scalar);

    MPI_Scatter(va, local_n, MPI_DOUBLE, vector_a, local_n, MPI_DOUBLE, 0, comm);
    MPI_Scatter(vb, local_n, MPI_DOUBLE, vector_b, local_n, MPI_DOUBLE, 0, comm);
  } else {
    MPI_Scatter(va, local_n, MPI_DOUBLE, vector_a, local_n, MPI_DOUBLE, 0, comm);
    MPI_Scatter(vb, local_n, MPI_DOUBLE, vector_b, local_n, MPI_DOUBLE, 0, comm);
  }
}


void vect_mult(
  double vector_a[],
  double vector_b[],
  double scalar,
  double v1_v2_buf[],
  double v1_s_buf[],
  double v2_s_buf[],
  int local_n,
  int n,
  MPI_Comm comm){

  double va[N], vb[N];
  int i;

  MPI_Allgather(vector_a, local_n, MPI_DOUBLE, va, local_n, MPI_DOUBLE, comm);
  MPI_Allgather(vector_b, local_n, MPI_DOUBLE, vb, local_n, MPI_DOUBLE, comm);

  for (i = 0; i < n; i++){
    v1_v2_buf[i] += va[i] * vb[i];
    v1_s_buf[i] += va[i] * scalar;
    v2_s_buf[i] += vb[i] * scalar;
  }
}

int main(){
  int my_rank, comm_sz, n = N, local_n;
  double vector1[N] = {0}, vector2[N] = {0}, scalar = 0.0, v1_v2 = 0.0;
  double v1_v2_buffer[N], v1_s_buffer[N], v2_s_buffer[N];

  MPI_Init(NULL, NULL);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

  local_n = n / comm_sz;

  Read_vector(vector1, vector2, &scalar, local_n, n, my_rank, MPI_COMM_WORLD);

  vect_mult(vector1, vector2, scalar, v1_v2_buffer, v1_s_buffer, v2_s_buffer, local_n, n, my_rank, MPI_COMM_WORLD);

  int i = 0;
  if (my_rank == 0){
#if 0
// just for debug
    printf("\nvector1 is :\n");
    for (i = 0; i < n; i++){
      printf("%g, ", vector1[i]);
    }

    printf("\nvector2 is :\n");
    for (i = 0; i < n; i++){
      printf("%g, ", vector2[i]);
    }

    printf("\nscalar is %g\n", scalar);

    for(i = 0; i < n; i++){
      printf("v1_v2_buffer[%d] is %g\n", i, v1_v2_buffer[i]);
      v1_v2 += v1_v2_buffer[i];
    }
#endif
    printf("\nvector1 * vector2 is %lf.\n", v1_v2);

    printf("\nvector1 * scalar:\n");
    for (i = 0; i < n; i++){
      printf("%lf ", v1_s_buffer[i]);
    }

    printf("\n\nvector2 * scalar:\n");
    for (i = 0; i < n; i++){
      printf("%lf ", v2_s_buffer[i]);
    }

    printf("\n");
  }

  MPI_Finalize();

  return 0;
}

3.11
求n个数和的表达式：x0+x1+…+x[n-1]
这n个数值的前缀和是n个部分和：x0, x0+x1, x0+x1+x2,…, x0+x1+…+x[n-1]
求前缀和事实上是求数值总和的一般化，它不只求出n个值的总和。
a. 设计一个串行算法来计算n个元素的前缀和。
b. 在有n个进程的系统和说那个并行化(a)小题中设计的串行程序，每个程序存储一个x_i值。
c. 设n=2^k, k为正整数。设计一个串行算法，然后将该算法并行化，使得这个并行算法仅需要k个通讯阶段。
d. MPI提供一个集合通信函数MPI_Scan，用来计算前缀和。该函数对有count个元素的数组进行操作，sendbuf_p和recvbuf_p都应该指向有count个datatype类型元素的数据块。op参数和MPI_Reduce中的op一样。编写一个MPI程序使每个MPI进程生成有count个元素的随机数数组，计算其前缀和并打印结果。

解答：
a

for (int i = 1; i < n; i++){
  a[i] += a[i - 1];
}

b

//如何将x值分发到每个线程上就不写了
//这里my_rank是进程号,假设my_rank号线程就是x[my_rank]
for (int i = 0; i < my_rank;i++){
  x_my_rank += x[i];
}

c
有k次的串行稍稍绕一点，开始以为是k次循环，后来是k个循环。
下面下一个n为8的情况

for(int i = 1; i < n; i += 2){
  x[i] += x[i - 1];
}

for(int i = 3; i < n; i += 4){
  x[i] += x[i - 2];
  x[i - 1] += x[i - 2];
}

for(int i = 7;i < n; i += 8){
  x[i] += x[i - 4];
  x[i - 1] += x[i - 4];
  x[i - 2] += x[i - 4];
  x[i - 3] += x[i - 4];
}

当然这里也可以用递归来写。

使用MPI的话，一个for循环算作1次通讯，这样就满足了题目的条件。

d.
使用MPI_Scan并不容易，困扰了好几天。
最后参考这里，并对原始代码进行修改，得到的如下代码。
这里和题目中不大一样的就是每一个进程都只有1个元素的随机数组。
这里暂时就不去纠结多个元素了。

#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>

/* these variables are used by prefix sum also, hence here */

int main(int argc, char **argv)
{
  int i, num, psum, ksum, *arr, *input;
  int P, PID;

  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &PID);
  MPI_Comm_size(MPI_COMM_WORLD, &P);

  arr = (int*)malloc(P*sizeof(int));

  /* 0 creates input, sends */
  if (PID == 0) {
    srand(getpid());
    input = (int*)malloc(P*sizeof(int));
    printf("Input is:");
    for (i = 0; i < P; i++) {
      input[i] = rand()%10;
      printf(" %d", input[i]);
    }
    printf("\n");
  }

  MPI_Scatter(input, 1, MPI_INT, &num, 1, MPI_INT, 0, MPI_COMM_WORLD);

  /* prefix sum over all processes */
  MPI_Scan(&num, &psum, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
  printf("%d: Scan: %d\n", PID, psum);

  free(arr);

  MPI_Finalize();

  exit(0);

}
/* prefix sum */

3.12
可以用环形传递来替代蝶形结构的全归约，在环形传递结构体中，如果有p个线程，每个进程q向进程q+1发送数据（进程p-1向进程0发送数据）。这一过程持续循环直至所有进程都获得理想的结果。我们可以用以下代码实现全归约：
sum = temp_val = my_val;
for (i=1;i < p; i++) {
MPI_Sendrecv_replace(&temp_val, 1, MPI_INT, dest, sendtag, source, recvtag, comm, &status);
sum += temp_val;
}
a. 编写一个MPIc恒旭实现这一算法。与蝶形结构的全归约相比，它的性能如何？
b. 修改刚才的程序，实现前缀求和的运算。

解答:
A (1) 先来实现这个算法：(用c++实现)

#include "mpi.h"

#include <vector>
#include <array>
#include <iostream>

int loop_pass(int argc, char* argv[])
{
  int myid, numprocs, left, right;

  //MPI_Request request;
  MPI_Status status;

  std::vector<int> buffer{ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };

  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
  MPI_Comm_rank(MPI_COMM_WORLD, &myid);

  right = (myid + 1) % numprocs;
  left = myid - 1;
  if (left < 0) {
    left = numprocs - 1;
  }

  double start = MPI_Wtime();
  int temp = buffer.at(myid);

  for (int i = 1; i < numprocs; ++i) {
    MPI_Sendrecv_replace(
      &temp,
      1,
      MPI_INT,
      right, 1,
      left, 1,
      MPI_COMM_WORLD,
      &status);
    buffer[myid] += temp;
  }
  double end = MPI_Wtime();

  std::cout << "Thread ID[ " << myid << "] : " << buffer.at(myid) << std::endl;

  MPI_Finalize();

  if (myid == 0) {
    std::cout << "Runtime = " << end - start << std::endl;
  }

  return 0;
}

int main(int argc, char* argv[]) {
  loop_pass(argc, argv);
}

B

#include "mpi.h"

#include <vector>
#include <array>
#include <iostream>

int loop_pass(int argc, char* argv[])
{
  int myid, numprocs, left, right;

  //MPI_Request request;
  MPI_Status status;

  std::vector<int> buffer{ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };

  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
  MPI_Comm_rank(MPI_COMM_WORLD, &myid);

  right = (myid + 1) % numprocs;
  left = myid - 1;
  if (left < 0) {
    left = numprocs - 1;
  }

  double start = MPI_Wtime();
  int temp = buffer.at(myid);

  for (int i = 1; i < numprocs; ++i) {
    MPI_Sendrecv_replace(
      &temp,
      1,
      MPI_INT,
      right, 1,
      left, 1,
      MPI_COMM_WORLD,
      &status);
    if (myid >= i) {
      buffer[myid] += temp;
    }
  }


  double end = MPI_Wtime();

  std::cout << "Thread ID[ " << myid << "] : " << buffer.at(myid) << std::endl;

  MPI_Finalize();

  if (myid == 0) {
    std::cout << "Runtime = " << end - start << std::endl;
  }

  return 0;
}

int main(int argc, char* argv[]) {
  loop_pass(argc, argv);
}

(未完)