2017-05-11 2 views
0

im Licht der KOMMENTAR EDITEDFehler bei Parallelisierung MPI_Allgather

I MPI lerne und ich einige Übungen mache einige Aspekte der es zu verstehen. Ich habe einen Code geschrieben, der ein einfaches Monte-Carlo ausführen soll.

Es gibt zwei Hauptschleifen, die durchgeführt werden müssen: eine über die Zeitschritte T und eine kleinere über die Anzahl der Moleküle N. Nachdem ich versucht habe, jedes Molekül zu bewegen, geht das Programm zum nächsten Zeitschritt.

Ich habe versucht, es zu parallelisieren, indem Sie die Operationen auf den Molekülen auf den verschiedenen Prozessoren teilen. Leider gibt der Code, der für 1 Prozessor funktioniert, die falschen Ergebnisse für aus, wenn p> 1 ist. Das Problem liegt wahrscheinlich in der folgenden Funktion und genauer ist durch einen Anruf an MPI_Allgather(local_r,n,MPI_DOUBLE,r,n,MPI_DOUBLE,MPI_COMM_WORLD);

Ich verstehe völlig nicht warum. Was mache ich falsch? (Neben einer primitiven Parallelisierungsstrategie)

Meine Logik war, dass ich für jeden Zeitschritt die Bewegungen auf den Molekülen auf den verschiedenen Prozessoren berechnen konnte. Leider, während ich mit den lokalen Vektoren local_r auf den verschiedenen Prozessoren arbeite, um die Energiedifferenz local_DE zu berechnen, brauche ich den globalen Vektor r, da die Energie des i-ten Moleküls von allen anderen abhängt. Daher dachte ich, MPI_Allgather zu nennen, da ich sowohl den globalen Vektor als auch die lokalen Vektoren aktualisieren muss.

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 ; 

} 

Hier ist es die komplette „arbeiten“ Code:

#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 ; 

} 
+1

Ein bedingter 'MPI_Allgather'? Würde das nicht Ihr Programm blockieren, falls es nicht in allen Prozessen aufgerufen wird? –

+0

Sie haben Recht, es aus dem zu entfernen, wenn der Block entfernt wird, leider funktioniert das Programm immer noch nicht – Fra

+0

Also, was ist das neue Problem? –

Antwort

2

Sie können nicht erwarten, die gleiche Energie unterschiedliche Anzahl von MPI Prozesse aus einem einfachen Grund gegeben zu bekommen - die Konfigurationen erzeugt sind sehr unterschiedlich, je nachdem wie viele Prozesse es gibt. Der Grund ist nicht MPI_Allgather, sondern die Art, wie die Monte-Carlo-Sweeps ausgeführt werden.

Bei einem Prozess versuchen Sie, Atom 1, dann Atom 2, dann Atom 3 usw. zu bewegen, bis Sie Atom N erreichen. Jeder Versuch sieht die Konfiguration aus der vorherigen, was in Ordnung ist.

Bei zwei Prozessen versuchen Sie, Atom 1 zu verschieben, während Sie gleichzeitig versuchen, Atom N/2 zu verschieben. Weder Atom 1 sieht die eventuelle Verschiebung von Atom N/2 noch umgekehrt, aber dann sehen die Atome 2 und N/2 + 1 die Verschiebung von Atom 1 und Atom N/2. Sie erhalten zwei Teilkonfigurationen, die Sie einfach mit all-gather zusammenführen. Dies ist nicht entspricht dem vorherigen Fall, wenn ein einzelner Prozess alle MC versucht. Gleiches gilt für den Fall von mehr als zwei Prozessen.

Es gibt eine weitere Quelle von Differenz - die Pseudozufallszahl (PRN) Sequenz. Die Sequenz, die durch die wiederholten Aufrufe von lrand48() in einem Prozess erzeugt wird, ist nicht dieselbe wie die kombinierte Sequenz, die durch mehrere unabhängige Aufrufe von lrand48() in verschiedenen Prozessen erzeugt wird. Daher unterscheidet sich die Akzeptanz auch bei sequenziellen Tests aufgrund der lokal unterschiedlichen PRN Sequenzen.

Vergessen Sie die spezifischen Werte der nach jedem Schritt erzeugten Energie. In einer richtigen MC-Simulation sind diese nicht signifikant. Was zählt, ist der Durchschnittswert über eine große Anzahl von Schritten. Diese sollten gleich sein (innerhalb eines bestimmten Bereichs proportional zu 1/sqrt(N)), unabhängig vom verwendeten Aktualisierungsalgorithmus.

+0

Bedeutet dies, dass dieser Code inhärent seriell ist? Schließlich sollten Atome nacheinander behandelt werden. – Shibli

+1

MCMC ist inhärent seriell, aber es gibt Methoden, die es ermöglichen, mehrere Markov-Ketten auszuführen und periodisch Informationen zwischen ihnen auszutauschen, wobei das Endergebnis sehr ähnlich zu dem ist, was eine sequentielle MCMC ergibt. Siehe [hier] (https://arxiv.org/pdf/1312.7479v1.pdf) für eine detaillierte Beschreibung einer Technik oder [hier] (https://www2.stat.duke.edu/~scs/Parallel.pdf)) für die vereinfachte Version. –

0

Es ist seit dem letzten Mal ziemlich lange ich MPI verwendet, aber es scheint, dass Ihr Programm anhält, wenn Sie versuchen, zu „sammeln“ und Aktualisieren Sie die Daten in mehreren aller Prozesse, und es ist unvorhersehbar, welche Prozesse für das Sammeln erforderlich sind.

In diesem Fall ist eine einfache Lösung, den Rest der Prozesse einige Dummy-Daten zu senden, damit sie einfach von anderen ignoriert werden können. Zum Beispiel,

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); 
} 

Dummy-Daten könnten einige außergewöhnliche falsche Zahlen sein, die nicht in den Ergebnissen sein sollten, so dass andere Prozesse sie herausfiltern könnten.

Aber wie ich bereits erwähnt habe, ist dies ziemlich verschwenderisch, da Sie nicht so viele Daten von allen Prozessen erhalten müssen, und wir möchten sie vor allem dann vermeiden, wenn sehr viele Daten gesendet werden müssen.

In diesem Fall können Sie einige "Flags" von anderen Prozessen sammeln, damit wir wissen, welche Prozesse eigene Daten senden.

// 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(...) 

ich mit MPI_Allgatherv erinnern, könnten Sie die Anzahl der Elemente angeben, von einem bestimmten Verfahren zu erhalten. Hier ist ein Beispiel: http://mpi.deino.net/mpi_functions/MPI_Allgatherv.html

Aber bedenken Sie, dies könnte ein Overkill sein, wenn das Programm nicht gut parallelisiert ist. In Ihrem Fall wird dies beispielsweise innerhalb einer Schleife platziert, sodass diese Prozesse ohne Daten immer noch auf das nächste Sammeln der Flags warten müssen.

+0

Ich habe versucht, Ihre Lösung, aber es scheint nicht zu funktionieren – Fra

+0

@Fra Dann schlage ich vor, Sie testen einige einfachere Beispiele, um zu sehen, ob es irgendwelche Probleme mit Ihrer Umgebung an erster Stelle. –

0

Sie sollten außerhalb für Schleife. Ich habe mit dem folgenden Code getestet, aber beachten Sie, dass ich die Zeilen mit RAND_MAX geändert habe, um konsistente Ergebnisse zu erhalten. Als Ergebnis gibt der Code die gleiche Antwort für die Anzahl der Prozessoren 1, 2 und 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 ; 
} 
+0

Dies ändert die Physik und erzeugt ein völlig anderes Ensemble, da das Atom "i + 1" die neue Position des Atoms "i" nicht sieht. –

+0

Ja, Sie haben Recht. – Shibli

Verwandte Themen