记MPI_Scatterv 和 MPI_Gatherv学习过程

最近在自学《并行程序设计导论》，这两天才看完第三章，正在做题，然后在一道题上花费了前所未有的四个小时，故记下此文…

题目：通过使用MPI_Scatterv以及MPI_Gatherv来实现向量之间的和以及点乘，且向量长度n可以不为线程数comm_sz的倍数。

这道题是对于MPI_Scater以及MPI_Gather的扩展。

首先让我们看一下最初的MPI_Scatter和MPI_Gather这两个函数的具体用法

int MPI_Scatter(const void    *sendbuf, 
				int          sendcount, 
				MPI_Datatype  sendtype,
                void          *recvbuf, 
                int          recvcount, 
                MPI_ Datatype recvtype,
                int               root,
               	MPI_Comm         comm);

在分析每一个参数之前，我们需要知道的是MPI_Scatter作为一个常用的数据分发函数，它的作用是将一个数组均分给所有的线程，所以很适合块划分模式。
举个例子：
我们现在有一个数组int A[6]={1,2,3,4,5,6}，现在要计算3·A

数组A：

假设我们给程序分配三个线程：

如果我们使用块划分的话，那么分配情况是这样的：

也就是以长度为2的块划分整个数组，然后顺序分配给所有线程。
了解了MPI_Scatter的功能，接下来再来看看它的参数。

int MPI_Scatter(const void    *sendbuf, 
				int          sendcount, 
				MPI_Datatype  sendtype,
                void          *recvbuf, 
                int          recvcount, 
                MPI_ Datatype recvtype,
                int               root,
               	MPI_Comm         comm);

• *sendbuf
  很明显，这个参数是我们要发送的数据的地址，也就是要均分的那个地址。

• sendcount
  这个参数代表的是我们要发送的数据的数据量

• sendtype
  这个参数很有意思，虽然它仅仅是用来说明我们要传送数据的数据类型，但是它的类型MPI_Datatype是mpi库自己开发的一套派生类型。

• *recvbuf,recvcount,recvtype
  它们与上面对应，分别代表的是接收的数据要存放的地址，它们的数据量和数据类型。

• root
  这个参数说明的是你的根线程，因为你必须要有一个线程来掌控全局，来进行数据的分发，一般都使用线程0。

MPI_Gather方法与之相对，即将MPI_Scatter分配的数据再发送回根进程中，在此不作阐述，仅提供参数详情，有兴趣可以动手实践或者自行查阅官方文档。

• MPI_Gather

int MPI_Gather(const void    *sendbuf, 
               int          sendcount,
               MPI_Datatype  sendtype,
               void          *recvbuf, 
               int          recvcount,             
               MPI_Datatype  recvtype,
               int               root,
               MPI_Comm         comm);

接下来展示一个简单的实例，最基础的向量点乘，只不过我们必须控制线程数可以被向量长度整除

//
// Created by 林庚 on 2021/4/21.
// File_name: dot.c
//
#include<stdlib.h>
#include<stdio.h>
#include<mpi.h>
void parallel_vector_sum(
        double  local_x[],
        double  local_y[],
        double  local_z[],
        int     local_n);
void read_vector(
        double  local_a[],
        int     local_n,
        int     n,
        char    vec_name[],
        int     my_rank,
        MPI_Comm comm       );
void print_vector(
        double  local_b[],
        int     local_n,
        int     n,
        char    title[],
        int     my_rank,
        MPI_Comm comm
);
int main(){
    int comm_sz=0;
    int my_rank=0;

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

    int n=12;
    double x[n];
    double y[n];
    double z[n];
    read_vector(x,n/comm_sz,n,"X",my_rank,MPI_COMM_WORLD);
    read_vector(y,n/comm_sz,n,"Y",my_rank,MPI_COMM_WORLD);

    parallel_vector_sum(x,y,z,n/comm_sz);

    print_vector(z,n/comm_sz,n,"result",my_rank,MPI_COMM_WORLD);

    MPI_Finalize();
    return 0;
}
/*----------------------------------------------------*/
void read_vector(
        double  local_a[],
        int     local_n,
        int     n,
        char    vec_name[],
        int     my_rank,
        MPI_Comm comm       ){
    double * a=NULL;
    int i;

    if (my_rank == 0){
        a = malloc(n*sizeof(double));
        printf("enter the vector %s\n",vec_name);
        fflush(stdout);
        for (i = 0; i < n; i++) {
            scanf("%lf",&a[i]);
            fflush(stdin);
        }
        MPI_Scatter(a,local_n,MPI_DOUBLE,local_a,local_n,MPI_DOUBLE,0,comm);
        free(a);
    } else{
        MPI_Scatter(a,local_n,MPI_DOUBLE,local_a,local_n,MPI_DOUBLE,0,comm);
        fflush(stdin);
    }
}
/*----------------------------------------------------*/
void parallel_vector_sum(
        double  local_x[],  //input
        double  local_y[],  //input
        double  local_z[],  //output
        int     local_n){
    int local_i;
    for (local_i = 0; local_i < local_n; ++local_i) {
        local_z[local_i]=local_x[local_i]+local_y[local_i];
    }
}
/*----------------------------------------------------*/
void print_vector(
        double  local_b[],
        int     local_n,
        int     n,
        char    title[],
        int     my_rank,
        MPI_Comm comm
){

    double *b=NULL;
    int i;
    if (my_rank==0){
        b=malloc(n*sizeof(double ));
        MPI_Gather(local_b,local_n,MPI_DOUBLE,b,local_n,MPI_DOUBLE,0,comm);
        printf("%s\n",title);
        fflush(stdout);
        for (i = 0; i <n ;i ++) {
            printf("%f ",b[i]);
            fflush(stdout);
        }
        printf("\n");
        fflush(stdout);
        free(b);
    } else{
        MPI_Gather(local_b,local_n,MPI_DOUBLE,b,local_n,MPI_DOUBLE,0,comm);
    }
}

在这个C程序中，我们将数组x,y使用块划分分享到各个线程当中，最后计算结果将输出到数组z中。
每一个进程只需要各调用一次parallel_vector_sum方法，输入参数对应的是MPI_Scatter分配的等量却不同的数据，结果保存在数组z中，最后在print_vector方法中调用MPI_Gather方法将结果收回根线程，并且输出。
最后运行结果如下：

(base)  % mpicc -o ./a.out dot.c
(base)  % mpiexec -n 4 ./a.out     
enter the vector X1 
2 3 4 5 6 7 8 9 10 11 12
enter the vector Y1 
2 3 4 5 6 7 8 9 10 11 12
result
2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 18.0 20.0 22.0 24.0 

接下来回到本文标题MPI_Scatterv

这个方法的作用是为了解决数据分配过程中只能均分的缺陷，话不多说我们直接来看参数详情

int MPI_Scatterv(const void      *sendbuf, 
                 const int   sendcounts[], 
                 const int       displs[],
                 MPI_Datatype    sendtype, 
                 void            *recvbuf, 
                 int            recvcount,
                 MPI_Datatype    recvtype,
                 int                 root, 
                 MPI_Comm           comm);

可以很明显看到，MPI_Scatterv方法中相较于之前的MPI_Scatter多了两个参数 sendcounts[] 和 displs[] ，并且少了 recvcount。
首先阐述一个简单的推理思路：
为了实现块划分数据分配且让线程数不必强制满足要被向量长度整除的条件，我们要做的就是可以让块划分变得不平均。那么我们就有必要将每一个块的长度记录下来，也就是每一次发送的数据量记录下来，这也就是sendcounts[]的作用。而displs[]则相对抽象一些，它里面储存的是偏移量，有了它可以使数据分配更加灵活。接下来通过举例来更好的说明它们所代表的的意义。
比如：
我们有数组A={1,2,3,4,5,6,7,8} 有三个线程0，1，2。
为了满足某个神秘操作，我们需要给0号进程分配{1,2,3},给1号进程分配{3,4,5,6}，给2号进程分配{4,5,6,7,8}。

那么我们可以设置传入到MPI_Scatterv的参数为： recvcounts={3,4,5} displs={0,2,3}
所以很容易知道，要保证数据分配过程中不出现问题，recvcounts和displs的长度必须等于线程的数量。

理解了recvcounts[]以及displs[]的作用后，对于MPI_Gatherv的作用理解也就是水到渠成的事情了。

在此对于自己耗费的四个小时感到痛心，只能说看官方文档的时候不仔细思考，就想着一遍遍测试，最后debug一直de了四个小时…

---

最后附上这道花了我4个h的题目的代码
来自于《并行程序设计导论》第三章习题3.13

改用MPI_Scatterv和MPI_Gatherv来实现向量的点乘，且向量长度n可以不为线程数comm_sz的倍数。

//
// Created by 林庚 on 2021/4/23.
//File_name 3.13.c
//
#include <stdio.h>
#include <mpi.h>
#include <stdlib.h>
void read_vector(int ,int[],int [],int);
void get_slides(int[],int[],int);
void dot(int[],int[],int[],int);
void print_result(int[],int,int);
void get_n(int*,int*);
int comm_sz;
int my_rank;
MPI_Comm comm;
int main(){
    MPI_Init(NULL,NULL);
    comm=MPI_COMM_WORLD;


    MPI_Comm_size(comm,&comm_sz);
    MPI_Comm_rank(comm,&my_rank);


    int n;
    int local_num;
    get_n(&n,&local_num);
    int local_a[local_num];
    int local_b[local_num];
    int local_result[local_num];


    read_vector(n,local_a,local_b,local_num);
    dot(local_a,local_b,local_result,local_num);
    print_result(local_result,n,local_num);


    MPI_Finalize();

}
/*----------------------------------------------------*/
void read_vector(int n,int local_a[],int local_b[],int local_num){
    int local_i;
    int *a=NULL;
    int *b=NULL;
    int *slides=NULL;
    int *counts=NULL;
    if (my_rank==0) {
        a = malloc(n * sizeof(int));
        b = malloc(n * sizeof(int));
        slides = malloc(comm_sz * sizeof(int));
        counts = malloc(comm_sz * sizeof(int));
        get_slides(slides, counts, n);
        printf("please input vec a\n");
        for (local_i = 0; local_i < n; ++local_i) {
            scanf("%d", a + local_i);
            fflush(stdin);
        }
        printf("please input vec b\n");
        for (local_i = 0; local_i < n; ++local_i) {
            scanf("%d", b + local_i);
            fflush(stdin);
        }
        MPI_Scatterv(a, counts, slides, MPI_INT, local_a, local_num, MPI_INT, 0, MPI_COMM_WORLD);
        MPI_Scatterv(b, counts, slides, MPI_INT, local_b, local_num, MPI_INT, 0, MPI_COMM_WORLD);
        printf("sent successfully\n");
        free(a);
        free(b);
        free(slides);
        free(counts);
    } else{
        MPI_Scatterv(a, counts, slides, MPI_INT, local_a, local_num, MPI_INT, 0, MPI_COMM_WORLD);
        MPI_Scatterv(b, counts, slides, MPI_INT, local_b, local_num, MPI_INT, 0, MPI_COMM_WORLD);
    }
    MPI_Barrier(MPI_COMM_WORLD);
}
/*----------------------------------------------------*/
void get_slides(int slides[],int counts[],int n){
    int local_i;
    int offset=0;
    for (local_i = 0; local_i < comm_sz; ++local_i) {
        if (local_i==comm_sz-1){
            counts[local_i]=n/comm_sz+n%comm_sz;
            slides[local_i]=offset;
            break;
        }
        counts[local_i]=n/comm_sz;
        slides[local_i]=offset;
        offset+=counts[local_i];
    }
}
/*----------------------------------------------------*/
void get_n(int *n,int *local_num){
    if (my_rank==0){
    printf("input n\n");
    scanf("%d",n);}
    MPI_Bcast(n,1,MPI_INT,0,MPI_COMM_WORLD);

    if (my_rank!=comm_sz-1){
        *local_num=*n/comm_sz;
    } else{
        *local_num=*n/comm_sz+*n%comm_sz;
    }
    printf("local_num of thread%d is %d\n",my_rank,*local_num);
}
/*----------------------------------------------------*/
void dot(int local_a[],int local_b[],int local_result[],int local_num){
    int local_i;
    for ( local_i = 0; local_i < local_num; ++local_i) {
        local_result[local_i]=local_a[local_i]*local_b[local_i];
    }
    printf("calculating of %d has finished... first of result is %d\n",my_rank,*local_result);
}
/*----------------------------------------------------*/
void print_result(int local_result[],int n,int local_num){
    int* a=NULL;
    int* counts=NULL;
    int* slides=NULL;
    int local_i;
    if (my_rank==0){
        a=malloc(n*sizeof(int));
        slides=malloc(comm_sz*sizeof(int));
        counts=malloc(comm_sz * sizeof(int));
        get_slides(slides, counts, n);
        MPI_Gatherv(local_result, local_num, MPI_INT, a, counts, slides, MPI_INT, 0, MPI_COMM_WORLD);
        printf("result:\n");
        for (local_i = 0; local_i < n; ++local_i) {
            printf("%d ",a[local_i]);
        }
        free(a);
        free(slides);
        free(counts);
    } else{
        MPI_Gatherv(local_result, local_num, MPI_INT, a, counts, slides, MPI_INT, 0, MPI_COMM_WORLD);
    }
}

最后的运行结果：

(base) % mpicc -o ./a.out 3.13.c
(base) % mpiexec -n 4 ./a.out   
input n
12
local_num of thread0 is 3
local_num of thread2 is 3
local_num of thread1 is 3
local_num of thread3 is 3
please input vec a
1 2 3 4 5 6 7 8 9 10 11 12
please input vec b
2 2 2 2 2 2 2 2 2 2 2 2 
sent successfully
calculating of 0 has finished... first of result is 2
calculating of 1 has finished... first of result is 8
calculating of 2 has finished... first of result is 14
calculating of 3 has finished... first of result is 20
result:
2 4 6 8 10 12 14 16 18 20 22 24 %     

---

本人水平不高，如若内容有误或者有疑问想要探讨，可以在评论区内留言