Overhauled MPIMCI and added method to load seeds from file

Author: Ithanil

examples/ex2/main.cpp

@@ -66,6 +66,9 @@ int main() {
     const int ndim = 1;
     MCI * mci = new MCI(ndim);
 
+    // if you want to seed the different MPI threads from a file, for reproducible results
+    MPIMCI::setSeed(mci, "rseed.txt");
+
     if (myrank == 0) cout << "ndim = " << mci->getNDim() << endl;
 
     // set the integration range to [-1:3]
@@ -84,7 +87,6 @@ int main() {
 
     if (myrank == 0) cout << "initial walker position = " << mci->getX(0) << endl;
 
-
     // initial MRT2 step
     double * step = new double[ndim];
     step[0] = 0.5;

examples/ex2/rseed.txt

@@ -0,0 +1 @@
+../../res/rseed.txt

src/MCIntegrator.cpp

@@ -526,7 +526,7 @@ void MCI::setIRange(const double * const * irange)
 }
 
 
-void MCI::setSeed(const long seed)
+void MCI::setSeed(const uint_fast64_t seed) // fastest unsigned integer which is at least 64 bit (as expected by rgen)
 {
     _rgen.seed(seed);
 }

src/MCIntegrator.hpp

@@ -88,7 +88,7 @@ class MCI
     ~MCI();  //Destructor
 
     // --- Setters
-    void setSeed(const long seed);
+    void setSeed(const uint_fast64_t seed);
 
     void setIRange(const double * const * irange);
 

src/MPIMCI.hpp

@@ -2,22 +2,63 @@
 #define MPIMCI
 
 #include <stdexcept>
+#include <string>
+#include <fstream>
 #include <mpi.h>
 #include "MCIntegrator.hpp"
 
 namespace MPIMCI
 {
-    int init() // return mpi rank of process
+    // return my rank
+    int rank()
+    {
+        int myrank;
+        MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+        return myrank;
+    }
+
+    // return size
+    int size()
+    {
+        int size;
+        MPI_Comm_size(MPI_COMM_WORLD, &size);
+        return size;
+    }
+
+    // init MPI and return rank of process
+    int init()
     {
         int isinit;
         MPI_Initialized(&isinit);
         if (isinit==1) throw std::runtime_error("MPI already initialized!");
         MPI_Init(NULL, NULL);
-        int myrank;
-        MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
-        return myrank;
+        return rank();
+    }
+
+    // set different random seeds per thread from a file
+    void setSeed(MCI * const mci, const std::string filename)
+    {
+        int myrank = rank();
+        int nranks = size();
+
+        uint_fast64_t seeds[nranks];
+        uint_fast64_t myseed;
+        if (myrank == 0) {
+            std::ifstream seedfile;
+            seedfile.open(filename);
+
+            for (int i=0; i<nranks; ++i) {
+                seedfile >> seeds[i];
+            }
+            seedfile.close();
+        }
+
+        MPI_Scatter(seeds, 1, MPI_UNSIGNED_LONG, &myseed, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
+        mci->setSeed(myseed);
     }
 
+
+    // integrate in parallel and accumulate results
     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) {
@@ -28,8 +69,7 @@ namespace MPIMCI
             MPI_Finalized(&isfinal);
             if (isfinal==1) throw std::runtime_error("MPI already finalized!");
 
-            int nranks;
-            MPI_Comm_size(MPI_COMM_WORLD, &nranks);
+            int nranks = size();
 
             // the results are stored in myAverage/Error and then reduced into the same average/error for all processes
             double myAverage[mci->getNObsDim()];
@@ -52,6 +92,7 @@ namespace MPIMCI
         }
     }
 
+    // finalize MPI
     void finalize()
     {
         // make sure the user has MPI in the correct state