【发布时间】：2017-10-10 12:23:16
【问题描述】：
根据评论编辑
我正在学习 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 ;
}
这是完整的“工作”代码：
#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 ;
}
【问题讨论】：
-
有条件的
MPI_Allgather?如果没有在所有进程中调用它，那不会阻止您的程序吗？ -
你把它从 if 移除块中取出来是对的，不幸的是程序仍然无法工作
-
那么新问题是什么？
-
当p改变total_E的值改变（如果ones改变随机数生成器的种子为固定数就清楚了），当p = 2时它仍然类似于p = 1但是当p =4 完全不同
-
我并没有真正得到您的化学或物理资料，但如果您想对所有流程的不可预测部分进行收集，您可以简单地让其余流程发送一些虚拟数据。有点浪费但直截了当。
标签： c mpi physics montecarlo markov-chains