hmc_tm: add real reproducible random numbers mode which will keep the RNG state consistent

across executions of different numbers of trajectories (say you have a serial run doing
1000 trajectories at a time in 10 goes and a parallel run doing 10000 trajectories
in one run, with this mode, the two runs will still have exactly the
same random numbers because the RNG state is saved to file between executions of hmc_tm)

Author: Bartosz Kostrzewa

default_input_values.h

@@ -54,8 +54,9 @@
 #define _default_g_beta 6.0
 #define _default_g_N_s 20
 #define _default_g_dflgcr_flag 0
+#define _default_save_ranlux_state 0
 #define _default_random_seed 123456
-#define _default_rlxd_level 1
+#define _default_rlxd_level 2
 #define _default_solver_flag 1
 #define _default_startoption 0
 #define _default_Ntherm 0
@@ -64,7 +65,8 @@
 #define _default_write_cp_flag 1
 #define _default_cp_interval 5
 #define _default_nstore 0
-#define _default_rlxd_input_filename "last_state"
+#define _default_rlxd_input_filename "rlxd_state"
+#define _default_rlxs_input_filename "rlxs_state"
 #define _default_gauge_input_filename "conf"
 #define _default_read_source_flag 0
 #define _default_source_filename "source"

doc/input.tex

@@ -79,7 +79,41 @@ \subsection{Input parameter for main program}
   See {\ttfamily NrXProcs}.
 
 \item {\ttfamily seed}:\\
-  The seed for the random number generator. Default value is $123456$.
+  The seed for the random number generator. Possible values: Integer 
+  $i \in \left[ 1, 2^{31}-1 \right] $ or {\ttfamily statefile}. Default value is $123456$.
+  The seeds for different MPI processes are computed from this number XOR'ed
+  with the store counter and a regular pattern according to the process number.
+  When {\ttfamily statefile} is specified here, ranlux will be initialised
+  from the state stored in the files {\ttfamily rlxd\_state} and {\ttfamily rlxs\_state}.
+  Also, {\ttfamily hmc\_tm} will save the ranlux state at the end of each
+  trajectory. Note that when using MPI, {\ttfamily seed = statefile} should only ever be used 
+  when in {\ttfamily ReproduceRandomNumbers} mode because all processes 
+  will be initialised with the same random number generator. Serial or OpenMP
+  executions of {\ttfamily hmc\_tm} have only one random number generator and
+  so this option allows running a simulation with a single chain of
+  random numbers fully specified by the initial seed set for the first
+  execution.
+
+\item {\ttfamily SaveRanluxState}:\\
+  Store the state of the random number generator of MPI process 0 
+  to the files {\ttfamily rlxd\_state} and {\ttfamily rlxs\_state} 
+  at the end of each trajectory. This is useful to make a fully reproducible
+  chain of random numbers even across multiple executions of {\ttfamily hmc\_tm}.
+  When using MPI, this only makes sense if {\ttfamily ReproduceRandomNumbers} is
+  used. (see {\ttfamily seed}) 
+
+\item {\ttfamily ReproduceRandomNumbers}:\\
+  Possible values are {\ttfamily yes} or {\ttfamily no}, default is {\ttfamily yes}.
+  When set to {\ttfamily yes}, when random numbers are requested the random
+  number generators of all processes are set to the same state. Then each
+  process generates random numbers for the whole volume but uses only those
+  which belong to its local volume. This ensures that a parallelised run
+  is equivalent to a serial one as far as the random numbers are concerned and
+  is useful for testing purposes.
+  This can also be extended by setting {\ttfamily seed = statefile }
+  to produce a single chain of random numbers across multiple executions
+  of the program (such as using {\ttfamily StartCondition = continue} and
+  {\ttfamily InitialStoreCounter = readin}, even when using MPI.
 
 \item {\ttfamily kappa}:\\
   The $\kappa$ value. Default is $0.12$. For the {\ttfamily hmc\_tm}

hmc_tm.c

@@ -290,7 +290,15 @@ int main(int argc,char *argv[]) {
 #endif
 
   /* Initialise random number generator */
-  start_ranlux(rlxd_level, random_seed^trajectory_counter);
+  /* if running in reproducible mode we can initialize ranlux using 
+   * a saved state */
+  if(reproduce_randomnumber_flag && random_seed == -1) {
+    start_ranlux_from_file(rlxd_input_filename,rlxs_input_filename);
+  } else if (random_seed == -1) {
+    fatal_error("Initializing RANLUX from file only works in reproducible random numbers mode. Aborting!","hmc_tm");
+  } else {  
+    start_ranlux(rlxd_level, random_seed^trajectory_counter);
+  }
 
   /* Set up the gauge field */
   /* continue and restart */

read_input.h

@@ -53,13 +53,15 @@ extern "C"
   extern int nstore;
   extern int crylov_space_dim;
   extern char rlxd_input_filename[500];
+  extern char rlxs_input_filename[500];
   extern char gauge_input_filename[500];
   extern int subforwilson_flag;
   extern int eigenvalue_method_flag;
   extern int eigenvalue_max_iterations;
   extern double eigenvalue_precision;
   extern int index_start;
   extern int index_end;
+  extern int save_ranlux_state;
   extern int random_seed;
   extern int rlxd_level;
   extern double X0, X1, X2, X3;

read_input.l

@@ -116,8 +116,10 @@ inline void rmQuotes(char *str){
   int nstore;
   int index_start, index_end;
   int random_seed;
+  int save_ranlux_state;
   int rlxd_level;
   char rlxd_input_filename[500];
+  char rlxs_input_filename[500];
   char gauge_input_filename[500];
   int read_source_flag;
   int return_check_flag, return_check_interval;
@@ -191,6 +193,7 @@ inline void rmQuotes(char *str){
 %x MU
 %x CSW
 %x SEED
+%x SAVERANLUXSTATE
 %x RLXDLEVEL
 %x NSAVE
 %x RLXDINPUTFILE
@@ -336,6 +339,7 @@ inline void rmQuotes(char *str){
 ^NoEigenvalues{EQL}                BEGIN(NOEV);
 ^EigenvaluePrecision{EQL}          BEGIN(PRECEV);
 ^seed{EQL}                         BEGIN(SEED);
+^SaveRanluxState{EQL}              BEGIN(SAVERANLUXSTATE);
 ^StartCondition{EQL}               BEGIN(STARTCOND);
 ^ThermalisationSweeps{EQL}         BEGIN(THERMSWEEPS);
 ^Measurements{EQL}                 BEGIN(NMEAS);
@@ -1622,10 +1626,23 @@ inline void rmQuotes(char *str){
   dfl_poly_iter=atoi(yytext);
   if(myverbose!=0) printf("dfl_poly_iter = %s \n", yytext);
 }
+<SAVERANLUXSTATE>yes {
+  save_ranlux_state = 1;
+  if(myverbose!=0) printf("Save RANLUX state at end of trajectory.\n");
+}
+<SAVERANLUXSTATE>no {
+  save_ranlux_state = 0;
+  if(myverbose!=0) printf("Don't save RANLUX state at end of trajectory.\n");
+}
 <SEED>{DIGIT}+               {
   random_seed=atoi(yytext);
   if(myverbose!=0) printf("seed=%s \n", yytext);
 }
+<SEED>statefile {
+  random_seed=-1;
+  save_ranlux_state=1;
+  if(myverbose!=0) printf("seed=%s; Trying to read ranlux state from file! Saving RANLUX state at end of trajectory!\n", yytext);
+}
 <RLXDLEVEL>[12] {
   rlxd_level = atoi(yytext);
   if(myverbose!=0) printf("RanluxdLevel set to %d \n", rlxd_level);
@@ -2113,6 +2130,7 @@ int read_input(char * conf_file){
   g_N_s = _default_g_N_s;
   g_dflgcr_flag = _default_g_dflgcr_flag;
   random_seed = _default_random_seed;
+  save_ranlux_state = _default_save_ranlux_state;
   rlxd_level = _default_rlxd_level;
   startoption = _default_startoption;
   Ntherm = _default_Ntherm;
@@ -2122,6 +2140,7 @@ int read_input(char * conf_file){
   cp_interval = _default_cp_interval;
   nstore = _default_nstore;
   strcpy(rlxd_input_filename, _default_rlxd_input_filename);
+  strcpy(rlxs_input_filename, _default_rlxs_input_filename);
   strcpy(gauge_input_filename, _default_gauge_input_filename);
   g_stdio_proc = _default_g_stdio_proc;
   index_start = _default_index_start;

start.c

@@ -77,6 +77,7 @@
 #include "ranlxd.h"
 #include "ranlxs.h"
 #include "start.h"
+#include "fatal_error.h"
 
 static void gauss_vector(double v[],int n)
 {
@@ -856,6 +857,69 @@ void start_ranlux(int level, int seed)
    rlxd_init(level, loc_seed);
 }
 
+/* read warning in start.h before using this function! */
+void start_ranlux_from_file(char * const rlxd_state_filename, char * const rlxs_state_filename) {
+  FILE * rlxd_state_file;
+  FILE * rlxs_state_file;
+
+  char error_message[1000];
+
+  rlxd_state_file = fopen(rlxd_state_filename,"r");
+  rlxs_state_file = fopen(rlxs_state_filename,"r");
+
+  if(rlxd_state_file != NULL) {
+    int rlxd_state[105];
+    fread(rlxd_state, sizeof(rlxd_state), 1, rlxd_state_file);
+    fclose(rlxd_state_file);
+    rlxd_reset(rlxd_state);
+  } else {
+    snprintf(error_message,1000,"Problem reading RLXD state file \"%s\", aborting!",rlxd_state_filename);
+    fatal_error(error_message,"start_ranlux_from_file");
+  }
+
+  if(rlxs_state_file != NULL) {
+    int rlxs_state[105];
+    fread(rlxs_state, sizeof(rlxs_state), 1, rlxs_state_file);
+    fclose(rlxs_state_file);
+    rlxs_reset(rlxs_state);
+  } else {
+    snprintf(error_message,1000,"Problem reading RLXS state file \"%s\", aborting!",rlxs_state_filename);
+    fatal_error(error_message,"start_ranlux_from_file");
+  }
+
+}
+
+void store_ranlux_state(char * const rlxd_state_filename, char * const rlxs_state_filename) {
+  FILE * rlxd_state_file;
+  FILE * rlxs_state_file;
+
+  char error_message[1000];
+
+  rlxd_state_file = fopen(rlxd_state_filename,"w");
+  rlxs_state_file = fopen(rlxs_state_filename,"w");
+
+  if(rlxd_state_file != NULL) {
+    int rlxd_state[105];
+    rlxd_get(rlxd_state);
+    fwrite(rlxd_state, sizeof(rlxd_state), 1, rlxd_state_file);
+    fclose(rlxd_state_file);
+  } else {
+    snprintf(error_message,1000,"Problem opening RLXD state file \"%s\" for writing, aborting!",rlxd_state_filename);
+    fatal_error(error_message,"store_ranlux_state");
+  }
+  
+  if(rlxs_state_file != NULL) {
+    int rlxs_state[105];
+    rlxs_get(rlxs_state);
+    fwrite(rlxs_state, sizeof(rlxs_state), 1, rlxs_state_file);
+    fclose(rlxs_state_file);
+  } else {
+    snprintf(error_message,1000,"Problem opening RLXS state file \"%s\" for writing, aborting!",rlxs_state_filename);
+    fatal_error(error_message,"store_ranlux_state");
+  }
+
+}
+
 void gen_test_spinor_field(spinor * const k, const int eoflag) {
 
   int ix,iy,effvol;

start.h

@@ -68,6 +68,18 @@ void source_spinor_field_point_from_file(spinor * const P, spinor * const Q, int
 
 void start_ranlux(int level,int seed);
 
+/* This function allows initializing RANLUX from a saved state. IMPORTANT NOTE BELOW:
+ * Because of the way that random numbers are used in tmLQCD, this function should only ever be
+ * used in "reproducible random numbers" mode with the full understanding that all processes
+ * will have exactly the same random number generators!!
+ * The main routines in start.c all accomodate this by making every process generate random
+ * numbers for the whole volume and only using those relevant for the local volume while throwing
+ * all others away.
+ * If some function requests random numbers via any of the utility functions or ranlx[d,s] directly
+ * without taking this fact into account, all processes will generate the same ones! You have been warned. */
+void start_ranlux_from_file(char * const rlxd_state_filename, char * const rlxs_state_filename);
+void store_ranlux_state(char * const rlxd_state_filename, char * const rlxs_state_filename);
+
 void gen_test_spinor_field(spinor * const k , const int eoflag);
 void write_test_spinor_field(spinor * const k , const int eoflag, char * postfix);
 #endif

update_tm.c

@@ -334,7 +334,7 @@ int update_tm(double *plaquette_energy, double *rectangle_energy,
 #endif
   etime=gettime();
 
-  /* printing data in the .data file */
+  /* printing data in the .data file and save ranlux state if the option is set */
   if(g_proc_id==0) {
     datafile = fopen(filename, "a");
     if (!bc_flag) { /* if Periodic Boundary Conditions */
@@ -361,6 +361,10 @@ int update_tm(double *plaquette_energy, double *rectangle_energy,
     fprintf(datafile, "\n");
     fflush(datafile);
     fclose(datafile);
+
+    if(save_ranlux_state) {
+      store_ranlux_state(rlxd_input_filename, rlxs_input_filename);
+    }
   }
   return(accept);
 }