并行化错误 MPI_Allgather

Question

根据评论进行编辑

我正在学习 MPI，我正在做一些练习以了解它的某些方面。我已经编写了一个应该执行简单蒙特卡洛的代码。

其中有两个主要循环必须完成：一个关于时间步长 T 和一个较小的关于分子数量 N 的循环。因此，在我尝试移动每个分子后，程序会进入下一个时间步。

我试图通过在不同处理器上划分对分子的操作来并行化它。不幸的是，适用于 1 个处理器的代码会在 p>1 时打印 total_E 的错误结果。问题可能在于以下函数，更准确地说是通过调用 MPI_Allgather(local_r,n,MPI_DOUBLE,r,n,MPI_DOUBLE,MPI_COMM_WORLD);

我完全不明白为什么。我究竟做错了什么？（除了原始的并行化策略）

我的逻辑是，对于每个时间步长，我都可以计算不同处理器上分子的移动。不幸的是，当我在各种处理器上使用局部向量 local_r 来计算能量差 local_DE 时，我需要全局向量 r，因为第 i 个分子的能量取决于在所有其他人身上。因此我想调用 MPI_Allgather 因为我必须更新全局向量和本地向量。

void Step(double (*H)(double,double),double* local_r,double* r,double *E_,int n,int my_rank){

    int i;
    double* local_rt = calloc(n,sizeof(double));
    double local_DE;

    for(i=0;i<n;i++){

        local_rt[i] = local_r[i] + delta*((double)lrand48()/RAND_MAX-0.5);
            local_rt[i] = periodic(local_rt[i]);

            local_DE = E_single(H,local_rt,r,i,n,my_rank) - E_single(H,local_r,r,i,n,my_rank);

        if ( local_DE <= 0.0 || exp(-local_DE) > (double) lrand48()/RAND_MAX  ) {               
            (*E_) += local_DE;
            local_r[i] = local_rt[i];

        }
MPI_Allgather(local_r,n,MPI_DOUBLE,r,n,MPI_DOUBLE,MPI_COMM_WORLD);

    }


    return ;

}

这是完整的 "working" 代码：

#define _XOPEN_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <mpi.h>

#define N 100
#define L 5.0
#define T_ 5000
#define delta 2.0


void Step(double (*)(double,double),double*,double*,double*,int,int);
double H(double ,double );
double E(double (*)(double,double),double* ,double*,int ,int );
double E_single(double (*)(double,double),double* ,double*,int ,int ,int);
double * pos_ini(void);
double periodic(double );
double dist(double , double );
double sign(double );

int main(int argc,char** argv){

    if (argc < 2) {
        printf("./program <outfile>\n");
        exit(-1);
    }

    srand48(0);
        int my_rank;
    int p;  
    FILE* outfile = fopen(argv[1],"w");
    MPI_Init(&argc,&argv);
    MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
    MPI_Comm_size(MPI_COMM_WORLD,&p);
    double total_E,E_;
    int n;
    n = N/p;
    int t;
    double * r =  calloc(N,sizeof(double)),*local_r = calloc(n,sizeof(double));

    for(t = 0;t<=T_;t++){
        if(t ==0){
            r = pos_ini();
            MPI_Scatter(r,n,MPI_DOUBLE, local_r,n,MPI_DOUBLE, 0, MPI_COMM_WORLD);
                E_ = E(H,local_r,r,n,my_rank);
        }else{
            Step(H,local_r,r,&E_,n,my_rank);
        }

        total_E = 0;
        MPI_Allreduce(&E_,&total_E,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);

            if(my_rank == 0){
                fprintf(outfile,"%d\t%lf\n",t,total_E/N);
            }

    }

    MPI_Finalize();

    return 0;

}



double sign(double a){

    if(a < 0){
        return -1.0 ;
    }else{
        return  1.0 ;
    }

}

double periodic(double a){

    if(sqrt(a*a) > L/2.0){
        a = a - sign(a)*L;
    }

    return a;
}


double dist(double a, double b){

    double d = a-b;
    d = periodic(d);

    return sqrt(d*d);
}

double * pos_ini(void){

  double * r  = calloc(N,sizeof(double));
  int i;

  for(i = 0;i<N;i++){
    r[i] =  ((double) lrand48()/RAND_MAX)*L - L/2.0;
  }

  return r;

}

double H(double a,double b){

      if(dist(a,b)<2.0){
        return  exp(-dist(a,b)*dist(a,b))/dist(a,b);
    }else{

    return 0.0;

    }
}



double E(double (*H)(double,double),double* local_r,double*r,int n,int my_rank){

    double local_V = 0;
    int i;

    for(i = 0;i<n;i++){
             local_V += E_single(H,local_r,r,i,n,my_rank);
     }
    local_V *= 0.5;

    return local_V; 
}

double E_single(double (*H)(double,double),double* local_r,double*r,int i,int n,int my_rank){

    double local_V = 0;
    int j;

      for(j = 0;j<N;j++){

        if( (i + n*my_rank) != j ){
            local_V+=H(local_r[i],r[j]);
        }

    }

    return local_V; 
}


void Step(double (*H)(double,double),double* local_r,double* r,double *E_,int n,int my_rank){

    int i;
    double* local_rt = calloc(n,sizeof(double));
    double local_DE;

    for(i=0;i<n;i++){

        local_rt[i] = local_r[i] + delta*((double)lrand48()/RAND_MAX-0.5);
            local_rt[i] = periodic(local_rt[i]);

            local_DE = E_single(H,local_rt,r,i,n,my_rank) - E_single(H,local_r,r,i,n,my_rank);

        if ( local_DE <= 0.0 || exp(-local_DE) > (double) lrand48()/RAND_MAX  ) {               
            (*E_) += local_DE;
            local_r[i] = local_rt[i];

        }
MPI_Allgather(local_r,n,MPI_DOUBLE,r,n,MPI_DOUBLE,MPI_COMM_WORLD);

    }


    return ;

}

Answer 1

自从我上次使用 MPI 以来已经很长时间了，但是当您尝试 "gather" 并更新所有进程中的数据时，您的程序似乎停止了，并且无法预测哪些进程会需要聚会。

所以在这种情况下，一个简单的解决方案是让其余进程发送一些虚拟数据，这样它们就可以被其他进程忽略。例如，

if ( local_DE <= 0.0 || exp(-local_DE) > (double) lrand48()/RAND_MAX  ) {               
    (*E_) += local_DE;
    local_r[i] = local_rt[i];
    MPI_Allgather(local_r,n,MPI_DOUBLE,r,n,MPI_DOUBLE,MPI_COMM_WORLD);
    // filter out the dummy data out of "r" here
} else {
    MPI_Allgather(dummy_sendbuf, n, MPI_DOUBLE, dummy_recvbuf, n, MPI_DOUBLE, MPI_COMM_WORLD);
}

虚拟数据可能是一些不应该出现在结果中的异常错误数字，因此其他进程可以将它们过滤掉。

但正如我所提到的，这非常浪费，因为您实际上不需要从所有进程接收那么多数据，我们希望避免这种情况，尤其是当有大量数据要发送时。

在这种情况下，您可以从其他进程收集一些 "flags" 以便我们可以知道哪些进程拥有要发送的数据。

// pseudo codes
// for example, place 1 at local_flags[my_rank] if it's got data to send, otherwise 0
MPI_Allgather(local_flags, n, MPI_BYTE, recv_flags, n, MPI_BYTE, MPI_COMM_WORLD)
// so now all the processes know which processes will send
// receive data from those processes
MPI_Allgatherv(...)

我记得MPI_Allgatherv, you could specify the number of elements to receive from a specific process. Here's an example: http://mpi.deino.net/mpi_functions/MPI_Allgatherv.html

但请记住，如果程序没有很好地并行化，这可能会有点矫枉过正。例如，在您的情况下，这被放置在一个循环中，因此那些没有数据的进程仍然需要等待下一次收集标志。

Answer 2

你应该把 MPI_Allgather() 放在 for 循环之外。我使用以下代码进行了测试，但请注意，我修改了涉及 RAND_MAX 的行以获得一致的结果。结果，代码对处理器 1、2 和 4 的数量给出了相同的答案。

void Step(double (*H)(double,double),double* local_r,double* r,double *E_,int n,int my_rank){

    int i;
    double* local_rt = calloc(n,sizeof(double));
    double local_DE;

    for(i=0;i<n;i++){

        //local_rt[i] = local_r[i] + delta*((double)lrand48()/RAND_MAX-0.5);
        local_rt[i] = local_r[i] + delta*((double)lrand48()-0.5);
        local_rt[i] = periodic(local_rt[i]);

        local_DE = E_single(H,local_rt,r,i,n,my_rank) - E_single(H,local_r,r,i,n,my_rank);

        //if ( local_DE <= 0.0 || exp(-local_DE) > (double) lrand48()/RAND_MAX  )
        if ( local_DE <= 0.0 || exp(-local_DE) > (double) lrand48()  )
        {
            (*E_) += local_DE;
            local_r[i] = local_rt[i];
        }

    }

    MPI_Allgather(local_r,n,MPI_DOUBLE,r,n,MPI_DOUBLE,MPI_COMM_WORLD);

    return ;
}

Answer 3

您不能指望在给定不同数量的 MPI 进程的情况下获得相同的能量，原因很简单 - 生成的配置非常不同，具体取决于进程的数量。原因不是 MPI_Allgather，而是蒙特卡洛扫描的执行方式。

给定一个过程，您尝试移动原子 1，然后是原子 2，然后是原子 3，依此类推，直到到达原子 N。每次尝试都会看到前一次的配置结果，这很好。

给定两个进程，您尝试移动原子 1，同时尝试移动原子 N/2。原子 1 看不到原子 N/2 的最终位移，反之亦然，但是原子 2 和 N/2+1 看不到原子 1 和原子 N/2 的最终位移。您最终得到两个部分配置，您只需将它们与全收集合并即可。这 not 等同于前一种情况，即单个进程完成所有 MC 尝试。两个以上进程的情况也是如此。

还有另一个差异来源 - 伪随机数 (PRN) 序列。在一个进程中重复调用 lrand48() 产生的序列与在不同进程中多次独立调用 lrand48() 产生的组合序列不同，因此即使你对试验进行顺序化，验收仍然会因本地不同的PRN序列而不同。

忘记每一步后产生的能量的具体值。在适当的 MC 模拟中，这些是微不足道的。重要的是大量步骤的平均值。无论使用何种更新算法，它们都应该相同（在与 1/sqrt(N) 成比例的一定范围内）。

并行化错误 MPI_Allgather

Error in parallelization MPI_Allgather

c

mpi

physics

montecarlo

markov-chains