Reduce->Allreduce, necessary for VMC

Author: Ithanil

examples/ex2/main.cpp

@@ -63,7 +63,7 @@ int main() {
         cout << "We start by initializing MPI and setting the MCI:" << endl;
     }
 
-    const int ndim = 1;    
+    const int ndim = 1;
     MCI * mci = new MCI(ndim);
 
     if (myrank == 0) cout << "ndim = " << mci->getNDim() << endl;
@@ -119,12 +119,9 @@ int main() {
 
     // integrate
     const long Nmc = 1000000;
-    double * average;
-    double * error;
-    if (myrank == 0) { // allocate only for root
-        average = new double[mci->getNObsDim()];
-        error = new double[mci->getNObsDim()];
-    }
+    double * average = new double[mci->getNObsDim()];
+    double * error = new double[mci->getNObsDim()];
+
     MPIMCI::integrate(mci, Nmc, average, error);
 
     if (myrank == 0) {
@@ -156,6 +153,16 @@ int main() {
     // integrate
     MPIMCI::integrate(mci, Nmc, average, error);
 
+    if (myrank == 0) {
+        cout << "The integral gives as result = " << average[0] << "   +-   " << error[0] << endl;
+        cout << "--------------------------------------------------------" << endl << endl;
+
+
+        // final comments
+        cout << "Using a sampling function in this case gives worse performance. In fact, the error bar is larger." << endl;
+        cout << "This implies that the variance of the re-factored f(x) written for introducing a sampling function, is larger than the original f(x)." << endl;
+    }
+
     // deallocate per-thread allocations
     delete sf;
 
@@ -170,22 +177,13 @@ int main() {
     delete[] irange[0];
     delete[] irange;
 
-    MPIMCI::finalize(mci); // also deletes all mci
-
-    // from now on just root thread
-    if (myrank == 0) {
-        cout << "The integral gives as result = " << average[0] << "   +-   " << error[0] << endl;
-        cout << "--------------------------------------------------------" << endl << endl;
-
+    delete[] average;
+    delete[] error;
 
-        // final comments
-        cout << "Using a sampling function in this case gives worse performance. In fact, the error bar is larger." << endl;
-        cout << "This implies that the variance of the re-factored f(x) written for introducing a sampling function, is larger than the original f(x)." << endl;
+    delete mci;
 
-        // deallocate
-        delete[] average;
-        delete[] error;
-    }
+    // finalize MPI
+    MPIMCI::finalize();
 
     // end
     return 0;

src/MPIMCI.hpp

@@ -15,46 +15,42 @@ namespace MPIMCI
         return MPI::COMM_WORLD.Get_rank();
     }
 
-    void integrate(MCI * const mci, const long &Nmc, double * average, double * error, bool findMRT2step=true, bool initialdecorrelation=true, bool use_mpi=true) // by setting use_mpi to false you can use this without requiring MPI 
+    void integrate(MCI * const mci, const long &Nmc, double * average, double * error, bool findMRT2step=true, bool initialdecorrelation=true, bool use_mpi=true) // by setting use_mpi to false you can use this without requiring MPI
     {
         if (use_mpi) {
             // make sure the user has MPI in the correct state
             if (!MPI::Is_initialized()) throw std::runtime_error("MPI not initialized!");
             if (MPI::Is_finalized()) throw std::runtime_error("MPI already finalized!");
 
-            const int myrank = MPI::COMM_WORLD.Get_rank();
             const int nranks = MPI::COMM_WORLD.Get_size();
 
-            // the results are stored in myAverage/Error and then reduced into average/error for root process
+            // the results are stored in myAverage/Error and then reduced into the same average/error for all processes
             double myAverage[mci->getNObsDim()];
             double myError[mci->getNObsDim()];
 
             mci->integrate(Nmc, myAverage, myError, findMRT2step, initialdecorrelation);
 
             for (int i=0; i<mci->getNObsDim(); ++i) {
-                MPI::COMM_WORLD.Reduce(&myAverage[i], &average[i], 1, MPI::DOUBLE, MPI::SUM, 0);
+                MPI::COMM_WORLD.Allreduce(&myAverage[i], &average[i], 1, MPI::DOUBLE, MPI::SUM);
 
                 myError[i] *= myError[i];
-                MPI::COMM_WORLD.Reduce(&myError[i], &error[i], 1, MPI::DOUBLE, MPI::SUM, 0);
+                MPI::COMM_WORLD.Allreduce(&myError[i], &error[i], 1, MPI::DOUBLE, MPI::SUM);
 
-                if (myrank == 0) {
-                    average[i] /= nranks;
-                    error[i] = sqrt(error[i]) / nranks;
-                }
+                average[i] /= nranks;
+                error[i] = sqrt(error[i]) / nranks;
             }
         }
         else {
             mci->integrate(Nmc, average, error, findMRT2step, initialdecorrelation); // regular single thread call
         }
     }
 
-    void finalize(MCI * mci)
+    void finalize()
     {
         // make sure the user has MPI in the correct state
         if (!MPI::Is_initialized()) throw std::runtime_error("MPI not initialized!");
         if (MPI::Is_finalized()) throw std::runtime_error("MPI already finalized!");
 
-        delete mci;
         MPI::Finalize();
     }
 };