c - 并行化错误 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);

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

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

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 ;

}

这是完整的“工作”代码:

#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 进程的情况下获得相同的能量 - 生成的配置非常不同，具体取决于进程的数量。原因不是 MPI_Allgather，而是蒙特卡洛扫描的执行方式。

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

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

还有另一个差异来源 - 伪随机数 (PRN) 序列。在一个进程中重复调用lrand48()产生的序列与在不同进程中多次独立调用lrand48()产生的组合序列是不一样的，因此即使您按顺序进行试验，由于局部不同的 PRN 序列，接受度仍然会有所不同。

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