CHROMA
const_hmc.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Main code for HMC with dynamical fermion generation
3  */
4 
5 #include "chroma.h"
6 #include <string>
7 
8 using namespace Chroma;
9 
10 namespace Chroma
11 {
12 
13  struct MCControl
14  {
16  QDP::Seed rng_seed;
17  unsigned long start_update_num;
18  unsigned long n_warm_up_updates;
19  unsigned long n_production_updates;
20  unsigned int n_updates_this_run;
21  unsigned int save_interval;
23  QDP_volfmt_t save_volfmt;
24  QDP_serialparallel_t save_pario;
28  bool rev_checkP;
31 
32  };
33 
34  void read(XMLReader& xml, const std::string& path, MCControl& p)
35  {
36  START_CODE();
37 
38  try {
39  XMLReader paramtop(xml, path);
40  p.cfg = readXMLGroup(paramtop, "Cfg", "cfg_type");
41  read(paramtop, "./RNG", p.rng_seed);
42  read(paramtop, "./StartUpdateNum", p.start_update_num);
43  read(paramtop, "./NWarmUpUpdates", p.n_warm_up_updates);
44  read(paramtop, "./NProductionUpdates", p.n_production_updates);
45  read(paramtop, "./NUpdatesThisRun", p.n_updates_this_run);
46  read(paramtop, "./SaveInterval", p.save_interval);
47  read(paramtop, "./SavePrefix", p.save_prefix);
48  read(paramtop, "./SaveVolfmt", p.save_volfmt);
49 
50  // -- Deal with parallel IO
51  p.save_pario = QDPIO_SERIAL; // Default
52 
53  // If there is a ParalelIO tag
54  if ( paramtop.count("./ParallelIO") > 0 ) {
55 
56  bool parioP=false;
57  read(paramtop, "./ParallelIO", parioP);
58  if ( parioP ) {
59  QDPIO::cout << "Setting parallel write mode" << std::endl;
60  p.save_pario = QDPIO_PARALLEL;
61  }
62  }
63 
64  // Default values: repro check is on, frequency is 10%
65  p.repro_checkP = true;
66  p.repro_check_frequency = 10;
67 
68 
69  // Now overwrite with user values
70  if( paramtop.count("./ReproCheckP") == 1 ) {
71  // Read user value if given
72  read(paramtop, "./ReproCheckP", p.repro_checkP);
73  }
74 
75  // If we do repro checking, read the frequency
76  if( p.repro_checkP ) {
77  if( paramtop.count("./ReproCheckFrequency") == 1 ) {
78  // Read user value if given
79  read(paramtop, "./ReproCheckFrequency", p.repro_check_frequency);
80  }
81  }
82 
83  // Reversibility checking enabled by default.
84  p.rev_checkP = true;
85  p.rev_check_frequency = 10;
86 
87  // Now overwrite with user values
88  if( paramtop.count("./ReverseCheckP") == 1 ) {
89  // Read user value if given
90  read(paramtop, "./ReverseCheckP", p.rev_checkP);
91  }
92 
93  // If we do repro checking, read the frequency
94  if( p.rev_checkP ) {
95  if( paramtop.count("./ReverseCheckFrequency") == 1 ) {
96  // Read user value if given
97  read(paramtop, "./ReverseCheckFrequency", p.rev_check_frequency);
98  }
99  }
100 
101  if( paramtop.count("./MonitorForces") == 1 ) {
102  read(paramtop, "./MonitorForces", p.monitorForcesP);
103  }
104  else {
105  p.monitorForcesP = true;
106  }
107 
108  if( paramtop.count("./InlineMeasurements") == 0 ) {
109  XMLBufferWriter dummy;
110  push(dummy, "InlineMeasurements");
111  pop(dummy); // InlineMeasurements
112  p.inline_measurement_xml = dummy.printCurrentContext();
113 
114  }
115  else {
116  XMLReader measurements_xml(paramtop, "./InlineMeasurements");
117  std::ostringstream inline_os;
118  measurements_xml.print(inline_os);
119  p.inline_measurement_xml = inline_os.str();
120  QDPIO::cout << "InlineMeasurements are: " << std::endl;
121  QDPIO::cout << p.inline_measurement_xml << std::endl;
122  }
123 
124 
125  }
126  catch(const std::string& e ) {
127  QDPIO::cerr << "Caught Exception: " << e << std::endl;
128  QDP_abort(1);
129  }
130 
131  END_CODE();
132  }
133 
134  void write(XMLWriter& xml, const std::string& path, const MCControl& p)
135  {
136  START_CODE();
137 
138  try {
139  push(xml, path);
140  xml << p.cfg.xml;
141  write(xml, "RNG", p.rng_seed);
142  write(xml, "StartUpdateNum", p.start_update_num);
143  write(xml, "NWarmUpUpdates", p.n_warm_up_updates);
144  write(xml, "NProductionUpdates", p.n_production_updates);
145  write(xml, "NUpdatesThisRun", p.n_updates_this_run);
146  write(xml, "SaveInterval", p.save_interval);
147  write(xml, "SavePrefix", p.save_prefix);
148  write(xml, "SaveVolfmt", p.save_volfmt);
149  {
150  bool pario = ( p.save_pario == QDPIO_PARALLEL );
151  write(xml, "ParallelIO", pario);
152  }
153  write(xml, "ReproCheckP", p.repro_checkP);
154  if( p.repro_checkP ) {
155  write(xml, "ReproCheckFrequency", p.repro_check_frequency);
156  }
157  write(xml, "ReverseCheckP", p.rev_checkP);
158  if( p.rev_checkP ) {
159  write(xml, "ReverseCheckFrequency", p.rev_check_frequency);
160  }
161  write(xml, "MonitorForces", p.monitorForcesP);
162 
163  xml << p.inline_measurement_xml;
164 
165  pop(xml);
166 
167  }
168  catch(const std::string& e ) {
169  QDPIO::cerr << "Caught Exception: " << e << std::endl;
170  QDP_abort(1);
171  }
172 
173  END_CODE();
174  }
175 
176 
177  struct HMCTrjParams
178  {
179  multi1d<int> nrow;
180 
181  // Polymorphic
182  std::string Monomials_xml; // XML for the monomials
183  std::string H_MC_xml; // XML for the Hamiltonian
184  std::string Integrator_xml; // XML for the Integrator
185  };
186 
187  void write(XMLWriter& xml, const std::string& path, const HMCTrjParams& p)
188  {
189  START_CODE();
190 
191  try {
192  push(xml, path);
193  write(xml, "nrow", p.nrow);
194  xml << p.Monomials_xml; // XML For the mon
195  xml << p.H_MC_xml;
196  xml << p.Integrator_xml;
197  pop(xml);
198  }
199  catch(const std::string& e ) {
200  QDPIO::cerr << "Caught Exception: " << e << std::endl;
201  QDP_abort(1);
202  }
203 
204  END_CODE();
205  }
206 
207 
208  void read(XMLReader& xml, const std::string& path, HMCTrjParams& p)
209  {
210  START_CODE();
211 
212  try {
213  XMLReader paramtop(xml, path);
214 
215  read(paramtop, "./nrow", p.nrow);
216  XMLReader Monomials_xml_reader(paramtop, "./Monomials");
217  std::ostringstream os_Monomials;
218  Monomials_xml_reader.print(os_Monomials);
219  p.Monomials_xml = os_Monomials.str();
220  QDPIO::cout << "Monomials xml is:" << std::endl;
221  QDPIO::cout << p.Monomials_xml << std::endl;
222 
223  // Now the XML for the Hamiltonians
224  XMLReader H_MC_xml(paramtop, "./Hamiltonian");
225  std::ostringstream os_H_MC;
226  H_MC_xml.print(os_H_MC);
227  p.H_MC_xml = os_H_MC.str();
228 
229  QDPIO::cout << "HMC_xml is: " << std::endl;
230  QDPIO::cout << p.H_MC_xml;
231 
232 
233  // Read the Integrator parameters
234  XMLReader MD_integrator_xml(paramtop, "./MDIntegrator");
235  std::ostringstream os_integrator;
236  MD_integrator_xml.print(os_integrator);
237  p.Integrator_xml = os_integrator.str();
238 
239  QDPIO::cout << "Integrator XML is: " << std::endl;
240  QDPIO::cout << p.Integrator_xml << std::endl;
241  }
242  catch( const std::string& e ) {
243  QDPIO::cerr << "Error reading XML : " << e << std::endl;
244  QDP_abort(1);
245  }
246 
247  END_CODE();
248  }
249 
250  template<typename UpdateParams>
251  void saveState(const UpdateParams& update_params,
252  MCControl& mc_control,
253  unsigned long update_no,
254  const multi1d<LatticeColorMatrix>& u) {
255  // Do nothing
256  }
257 
258  // Specialise
259  template<>
260  void saveState(const HMCTrjParams& update_params,
261  MCControl& mc_control,
262  unsigned long update_no,
263  const multi1d<LatticeColorMatrix>& u)
264  {
265  START_CODE();
266 
267  // File names
268  std::ostringstream restart_data_filename;
269  restart_data_filename << mc_control.save_prefix << "_restart_" << update_no << ".xml" ;
270 
271  std::ostringstream restart_config_filename;
272  restart_config_filename << mc_control.save_prefix << "_cfg_" << update_no << ".lime";
273 
274  XMLBufferWriter restart_data_buffer;
275 
276 
277  // Copy old params
278  MCControl p_new = mc_control;
279 
280  // Get Current RNG Seed
281  QDP::RNG::savern(p_new.rng_seed);
282 
283  // Set the current traj number
284  p_new.start_update_num = update_no;
285 
286  // Set the num_updates_this_run
287  unsigned long total = mc_control.n_warm_up_updates
288  + mc_control.n_production_updates ;
289 
290  if ( total < mc_control.n_updates_this_run + update_no ) {
291  p_new.n_updates_this_run = total - update_no;
292  }
293 
294  // Set the name and type of the config
295  {
296  // Parse the cfg XML including the parallel IO part
298 
299  // Reset the filename in it
300  cfg.cfg_file = restart_config_filename.str();
301  cfg.cfg_pario = mc_control.save_pario;
302 
303  // Prepare to write out
305  }
306 
307 
308  push(restart_data_buffer, "Params");
309  write(restart_data_buffer, "MCControl", p_new);
310  write(restart_data_buffer, "HMCTrj", update_params);
311  pop(restart_data_buffer);
312 
313 
314  // Save the config
315 
316  // some dummy header for the file
317  XMLBufferWriter file_xml;
318  push(file_xml, "HMC");
319  proginfo(file_xml);
320  pop(file_xml);
321 
322 
323  // Save the config
324  writeGauge(file_xml,
325  restart_data_buffer,
326  u,
327  restart_config_filename.str(),
328  p_new.save_volfmt,
329  p_new.save_pario);
330 
331 
332  // Write a restart DATA file from the buffer XML
333  // Do this after the config is written, so that if the cfg
334  // write fails, there is no restart file...
335  //
336  // production will then likely fall back to last good pair.
337 
338  XMLFileWriter restart_xml(restart_data_filename.str().c_str());
339  restart_xml << restart_data_buffer;
340  restart_xml.close();
341 
342  END_CODE();
343  }
344 
345 
346 
347  // Predeclare this
348  bool checkReproducability( const multi1d<LatticeColorMatrix>& P_new,
349  const multi1d<LatticeColorMatrix>& Q_new,
350  const QDP::Seed& seed_new,
351  const multi1d<LatticeColorMatrix>& P_old,
352  const multi1d<LatticeColorMatrix>& Q_old,
353  const QDP::Seed& seed_old );
354 
355  template<typename UpdateParams>
356  void doHMC(multi1d<LatticeColorMatrix>& u,
357  AbsHMCTrj<multi1d<LatticeColorMatrix>,
358  multi1d<LatticeColorMatrix> >& theHMCTrj,
359  MCControl& mc_control,
360  const UpdateParams& update_params,
361  multi1d< Handle<AbsInlineMeasurement> >& user_measurements)
362  {
363  START_CODE();
364 
365 
366  // Turn monitoring off/on
367  QDPIO::cout << "Setting Force monitoring to " << mc_control.monitorForcesP << std::endl;
368  setForceMonitoring(mc_control.monitorForcesP) ;
369  QDP::StopWatch swatch;
370 
371  XMLWriter& xml_out = TheXMLOutputWriter::Instance();
372  XMLWriter& xml_log = TheXMLLogWriter::Instance();
373 
374  push(xml_out, "doHMC");
375  push(xml_log, "doHMC");
376 
377  multi1d< Handle< AbsInlineMeasurement > > default_measurements(1);
378  InlinePlaquetteEnv::Params plaq_params;
379  plaq_params.frequency = 1;
380  // It is a handle
381  default_measurements[0] = new InlinePlaquetteEnv::InlineMeas(plaq_params);
382 
383  {
384  // Initialise the RNG
385  QDP::RNG::setrn(mc_control.rng_seed);
386 
387  // Fictitious momenta for now
388  multi1d<LatticeColorMatrix> p(Nd);
389 
390  // Create a field state
391  GaugeFieldState gauge_state(p,u);
392 
393  // Set the update number
394  unsigned long cur_update=mc_control.start_update_num;
395 
396  // Compute how many updates to do
397  unsigned long total_updates = mc_control.n_warm_up_updates
398  + mc_control.n_production_updates;
399 
400  unsigned long to_do = 0;
401  if ( total_updates >= mc_control.n_updates_this_run + cur_update +1 ) {
402  to_do = mc_control.n_updates_this_run;
403  }
404  else {
405  to_do = total_updates - cur_update ;
406  }
407 
408  QDPIO::cout << "MC Control: About to do " << to_do << " updates" << std::endl;
409 
410  // XML Output
411  push(xml_out, "MCUpdates");
412  push(xml_log, "MCUpdates");
413 
414  for(int i=0; i < to_do; i++)
415  {
416  push(xml_out, "elem"); // Caller writes elem rule
417  push(xml_log, "elem"); // Caller writes elem rule
418 
419  push(xml_out, "Update");
420  push(xml_log, "Update");
421  // Increase current update counter
422  cur_update++;
423 
424  // Decide if the next update is a warm up or not
425  bool warm_up_p = cur_update <= mc_control.n_warm_up_updates;
426  QDPIO::cout << "Doing Update: " << cur_update << " warm_up_p = " << warm_up_p << std::endl;
427 
428  // Log
429  write(xml_out, "update_no", cur_update);
430  write(xml_log, "update_no", cur_update);
431 
432  write(xml_out, "WarmUpP", warm_up_p);
433  write(xml_log, "WarmUpP", warm_up_p);
434 
435  bool do_reverse = false;
436  if( mc_control.rev_checkP
437  && ( cur_update % mc_control.rev_check_frequency == 0 )) {
438  do_reverse = true;
439  QDPIO::cout << "Doing Reversibility Test this traj" << std::endl;
440  }
441 
442 
443  // Check if I need to do any reproducibility testing
444  if( mc_control.repro_checkP
445  && (cur_update % mc_control.repro_check_frequency == 0 )
446  ) {
447 
448  // Yes - reproducibility trajectory
449  // Save fields and RNG at start of trajectory
450  QDPIO::cout << "Saving start config and RNG seed for reproducability test" << std::endl;
451 
452  GaugeFieldState repro_bkup_start( gauge_state.getP(), gauge_state.getQ());
453  QDP::Seed rng_seed_bkup_start;
454  QDP::RNG::savern(rng_seed_bkup_start);
455 
456  // DO the trajectory
457  QDPIO::cout << "Before HMC trajectory call" << std::endl;
458  swatch.reset();
459  swatch.start();
460 
461  // This may do a reversibility check
462  theHMCTrj( gauge_state, warm_up_p, do_reverse );
463  swatch.stop();
464 
465  QDPIO::cout << "After HMC trajectory call: time= "
466  << swatch.getTimeInSeconds()
467  << " secs" << std::endl;
468 
469  write(xml_out, "seconds_for_trajectory", swatch.getTimeInSeconds());
470  write(xml_log, "seconds_for_trajectory", swatch.getTimeInSeconds());
471 
472  // Save the fields and RNG at the end
473  QDPIO::cout << "Saving end config and RNG seed for reproducability test" << std::endl;
474  GaugeFieldState repro_bkup_end( gauge_state.getP(), gauge_state.getQ());
475  QDP::Seed rng_seed_bkup_end;
476  QDP::RNG::savern(rng_seed_bkup_end);
477 
478  // Restore the starting field and the starting RNG
479  QDPIO::cout << "Restoring start config and RNG for reproducability test" << std::endl;
480 
481  gauge_state.getP() = repro_bkup_start.getP();
482  gauge_state.getQ() = repro_bkup_start.getQ();
483  QDP::RNG::setrn(rng_seed_bkup_start);
484 
485  // Do the repro trajectory
486  QDPIO::cout << "Before HMC repro trajectory call" << std::endl;
487  swatch.reset();
488  swatch.start();
489  // Dont repeat the reversibility check in the repro test
490  theHMCTrj( gauge_state, warm_up_p, false );
491  swatch.stop();
492 
493  QDPIO::cout << "After HMC repro trajectory call: time= "
494  << swatch.getTimeInSeconds()
495  << " secs" << std::endl;
496 
497  write(xml_out, "seconds_for_repro_trajectory", swatch.getTimeInSeconds());
498  write(xml_log, "seconds_for_repro_trajectory", swatch.getTimeInSeconds());
499 
500  // Save seed at end of traj for comparison
501  QDP::Seed rng_seed_end2;
502  QDP::RNG::savern(rng_seed_end2);
503 
504 
505  // Check the reproducibility
506  bool pass = checkReproducability( gauge_state.getP(),
507  gauge_state.getQ(),
508  rng_seed_end2,
509  repro_bkup_end.getP(),
510  repro_bkup_end.getQ(),
511  rng_seed_bkup_end);
512 
513 
514  if( !pass ) {
515  QDPIO::cout << "Reproducability check failed on update " << cur_update << std::endl;
516  QDPIO::cout << "Aborting" << std::endl;
517  write(xml_out, "ReproCheck", pass);
518  write(xml_log, "ReproCheck", pass);
519  QDP_abort(1);
520  }
521  else {
522  QDPIO::cout << "Reproducability check passed on update " << cur_update << std::endl;
523  write(xml_out, "ReproCheck", pass);
524  write(xml_log, "ReproCheck", pass);
525  }
526 
527 
528  }
529  else {
530 
531  // Do the trajectory without accepting
532  QDPIO::cout << "Before HMC trajectory call" << std::endl;
533  swatch.reset();
534  swatch.start();
535  theHMCTrj( gauge_state, warm_up_p, do_reverse );
536  swatch.stop();
537 
538  QDPIO::cout << "After HMC trajectory call: time= "
539  << swatch.getTimeInSeconds()
540  << " secs" << std::endl;
541 
542  write(xml_out, "seconds_for_trajectory", swatch.getTimeInSeconds());
543  write(xml_log, "seconds_for_trajectory", swatch.getTimeInSeconds());
544 
545  }
546  swatch.reset();
547  swatch.start();
548 
549  // Create a gauge header for inline measurements.
550  // Since there are defaults always measured, we must always
551  // create a header.
552  //
553  // NOTE: THIS HEADER STUFF NEEDS A LOT MORE THOUGHT
554  //
555  QDPIO::cout << "HMC: start inline measurements" << std::endl;
556  {
557  XMLBufferWriter gauge_xml;
558  push(gauge_xml, "ChromaHMC");
559  write(gauge_xml, "update_no", cur_update);
560  write(gauge_xml, "HMCTrj", update_params);
561  pop(gauge_xml);
562 
563  // Reset and set the default gauge field
564  QDPIO::cout << "HMC: initial resetting default gauge field" << std::endl;
566  QDPIO::cout << "HMC: set default gauge field" << std::endl;
567  InlineDefaultGaugeField::set(gauge_state.getQ(), gauge_xml);
568  QDPIO::cout << "HMC: finished setting default gauge field" << std::endl;
569 
570  // Measure inline observables
571  push(xml_out, "InlineObservables");
572 
573  // Always measure defaults
574  for(int m=0; m < default_measurements.size(); m++)
575  {
576  QDPIO::cout << "HMC: do default measurement = " << m << std::endl;
577  QDPIO::cout << "HMC: dump named objects" << std::endl;
578  TheNamedObjMap::Instance().dump();
579 
580  // Caller writes elem rule
581  AbsInlineMeasurement& the_meas = *(default_measurements[m]);
582  push(xml_out, "elem");
583  the_meas(cur_update, xml_out);
584  pop(xml_out);
585 
586  QDPIO::cout << "HMC: finished default measurement = " << m << std::endl;
587  }
588 
589  // Only measure user measurements after warm up
590  if( ! warm_up_p )
591  {
592  QDPIO::cout << "Doing " << user_measurements.size()
593  <<" user measurements" << std::endl;
594  for(int m=0; m < user_measurements.size(); m++)
595  {
596  QDPIO::cout << "HMC: considering user measurement number = " << m << std::endl;
597  AbsInlineMeasurement& the_meas = *(user_measurements[m]);
598  if( cur_update % the_meas.getFrequency() == 0 )
599  {
600  // Caller writes elem rule
601  push(xml_out, "elem");
602  QDPIO::cout << "HMC: calling user measurement number = " << m << std::endl;
603  the_meas(cur_update, xml_out);
604  QDPIO::cout << "HMC: finished user measurement number = " << m << std::endl;
605  pop(xml_out);
606  }
607  }
608  QDPIO::cout << "HMC: finished user measurements" << std::endl;
609  }
610  pop(xml_out); // pop("InlineObservables");
611 
612  // Reset the default gauge field
613  QDPIO::cout << "HMC: final resetting default gauge field" << std::endl;
615  QDPIO::cout << "HMC: finished final resetting default gauge field" << std::endl;
616  }
617 
618  swatch.stop();
619  QDPIO::cout << "After all measurements: time= "
620  << swatch.getTimeInSeconds()
621  << " secs" << std::endl;
622 
623  write(xml_out, "seconds_for_measurements", swatch.getTimeInSeconds());
624  write(xml_log, "seconds_for_measurements", swatch.getTimeInSeconds());
625 
626  if( cur_update % mc_control.save_interval == 0 )
627  {
628  swatch.reset();
629  swatch.start();
630 
631  // Save state
632  saveState<UpdateParams>(update_params, mc_control, cur_update, gauge_state.getQ());
633 
634  swatch.stop();
635  QDPIO::cout << "After saving state: time= "
636  << swatch.getTimeInSeconds()
637  << " secs" << std::endl;
638  }
639 
640  pop(xml_log); // pop("Update");
641  pop(xml_out); // pop("Update");
642 
643  pop(xml_log); // pop("elem");
644  pop(xml_out); // pop("elem");
645  }
646 
647  // Save state
648  saveState<UpdateParams>(update_params, mc_control, cur_update, gauge_state.getQ());
649 
650  pop(xml_log); // pop("MCUpdates")
651  pop(xml_out); // pop("MCUpdates")
652  }
653 
654  pop(xml_log); // pop("doHMC")
655  pop(xml_out); // pop("doHMC")
656 
657  END_CODE();
658  }
659 
660  bool linkageHack(void)
661  {
662  bool foo = true;
663 
664  // Gauge Monomials
667 
668  // Ferm Monomials
670 
671  // MD Integrators
673 
674  // Chrono predictor
676 
677  // Inline Measurements
679 
680  // Gauge initialization
681  foo &= GaugeInitEnv::registerAll();
682 
683  return foo;
684  }
685 }
686 
687 using namespace Chroma;
688 
689 //! Hybrid Monte Carlo
690 /*! \defgroup hmcmain Hybrid Monte Carlo
691  * \ingroup main
692  *
693  * Main program for dynamical fermion generation
694 xml */
695 
696 int main(int argc, char *argv[])
697 {
698  Chroma::initialize(&argc, &argv);
699 
700  START_CODE();
701 
702  // Chroma Init stuff -- Open DATA and XMLDAT
703  QDPIO::cout << "Linkage = " << linkageHack() << std::endl;
704 
705  StopWatch snoop;
706  snoop.reset();
707  snoop.start();
708 
709  XMLFileWriter& xml_out = Chroma::getXMLOutputInstance();
710  XMLFileWriter& xml_log = Chroma::getXMLLogInstance();
711 
712  push(xml_out, "hmc");
713  push(xml_log, "hmc");
714 
715  HMCTrjParams trj_params;
716  MCControl mc_control;
717 
718  try
719  {
720  XMLReader xml_in(Chroma::getXMLInputFileName());
721 
722  XMLReader paramtop(xml_in, "/Params");
723  read( paramtop, "./HMCTrj", trj_params);
724  read( paramtop, "./MCControl", mc_control);
725 
726  // Write out the input
727  write(xml_out, "Input", xml_in);
728  write(xml_log, "Input", xml_in);
729  }
730  catch(const std::string& e) {
731  QDPIO::cerr << "hmc: Caught Exception while reading file: " << e << std::endl;
732  QDP_abort(1);
733  }
734 
735  if (mc_control.start_update_num >= mc_control.n_production_updates)
736  {
737  QDPIO::cout << "hmc: run is finished" << std::endl;
738  pop(xml_log);
739  pop(xml_out);
740  exit(0);
741  }
742 
743  QDPIO::cout << "Call QDP create layout" << std::endl;
744  Layout::setLattSize(trj_params.nrow);
745  Layout::create();
746  QDPIO::cout << "Finished with QDP create layout" << std::endl;
747 
748  proginfo(xml_out); // Print out basic program info
749  proginfo(xml_log); // Print out basic program info
750 
751  // Start up the config
752  multi1d<LatticeColorMatrix> u(Nd);
753  try
754  {
755  XMLReader file_xml;
756  XMLReader config_xml;
757 
758  QDPIO::cout << "Initialize gauge field" << std::endl;
759  StopWatch swatch;
760  swatch.reset();
761  swatch.start();
762  {
763  std::istringstream xml_c(mc_control.cfg.xml);
764  XMLReader cfgtop(xml_c);
765  QDPIO::cout << "Gauge initialization: cfg_type = " << mc_control.cfg.id << std::endl;
766 
768  gaugeInit(TheGaugeInitFactory::Instance().createObject(mc_control.cfg.id,
769  cfgtop,
770  mc_control.cfg.path));
771  (*gaugeInit)(file_xml, config_xml, u);
772  }
773  swatch.stop();
774  QDPIO::cout << "Gauge field successfully initialized: time= "
775  << swatch.getTimeInSeconds()
776  << " secs" << std::endl;
777 
778  swatch.reset();
779  swatch.start();
780  {
781  for(int mu=0; mu < Nd; mu++) {
782  reunit(u[mu]);
783  }
784  }
785  swatch.stop();
786  QDPIO::cout << "Gauge field reunitarized: time="
787  << swatch.getTimeInSeconds()
788  << " secs" << std::endl;
789 
790  // Write out the config header
791  write(xml_out, "Config_info", config_xml);
792  write(xml_log, "Config_info", config_xml);
793  }
794  catch(std::bad_cast)
795  {
796  QDPIO::cerr << "hmc: caught cast error" << std::endl;
797  QDP_abort(1);
798  }
799  catch(const std::string& e)
800  {
801  QDPIO::cerr << "hmc: Caught Exception: " << e << std::endl;
802  QDP_abort(1);
803  }
804  catch(std::exception& e)
805  {
806  QDPIO::cerr << "hmc: Caught standard library exception: " << e.what() << std::endl;
807  QDP_abort(1);
808  }
809  catch(...)
810  {
811  QDPIO::cerr << "hmc: caught generic exception during measurement" << std::endl;
812  QDP_abort(1);
813  }
814 
815 
816  // Set up the monomials
817  try {
818  std::istringstream Monomial_is(trj_params.Monomials_xml);
819  XMLReader monomial_reader(Monomial_is);
820  readNamedMonomialArray(monomial_reader, "/Monomials");
821  }
822  catch(const std::string& e) {
823  QDPIO::cout << "Caught Exception reading Monomials" << std::endl;
824  QDP_abort(1);
825  }
826 
827  std::istringstream H_MC_is(trj_params.H_MC_xml);
828  XMLReader H_MC_xml(H_MC_is);
829  ExactHamiltonianParams ham_params(H_MC_xml, "/Hamiltonian");
830 
832  multi1d<LatticeColorMatrix> > > H_MC(new ExactHamiltonian(ham_params));
833 
834 
835  std::istringstream MDInt_is(trj_params.Integrator_xml);
836  XMLReader MDInt_xml(MDInt_is);
837  LCMToplevelIntegratorParams int_par(MDInt_xml, "/MDIntegrator");
839  multi1d<LatticeColorMatrix> > > Integrator(new LCMToplevelIntegrator(int_par));
840 
841 
842  ConstLatColMatHMCTrj theHMCTrj( H_MC, Integrator );
843 
844 
845  multi1d < Handle< AbsInlineMeasurement > > the_measurements;
846 
847  // Get the measurements
848  try
849  {
850  std::istringstream Measurements_is(mc_control.inline_measurement_xml);
851 
852  XMLReader MeasXML(Measurements_is);
853 
854  std::ostringstream os;
855  MeasXML.print(os);
856  QDPIO::cout << os.str() << std::endl << std::flush;
857 
858  read(MeasXML, "/InlineMeasurements", the_measurements);
859  }
860  catch(const std::string& e) {
861  QDPIO::cerr << "hmc: Caught exception while reading measurements: " << e << std::endl
862  << std::flush;
863 
864  QDP_abort(1);
865  }
866 
867  QDPIO::cout << "There are " << the_measurements.size() << " user measurements " << std::endl;
868 
869 
870  // Run
871  try {
872  doHMC<HMCTrjParams>(u, theHMCTrj, mc_control, trj_params, the_measurements);
873  }
874  catch(std::bad_cast)
875  {
876  QDPIO::cerr << "HMC: caught cast error" << std::endl;
877  QDP_abort(1);
878  }
879  catch(std::bad_alloc)
880  {
881  // This might happen on any node, so report it
882  std::cerr << "HMC: caught bad memory allocation" << std::endl;
883  QDP_abort(1);
884  }
885  catch(const std::string& e)
886  {
887  QDPIO::cerr << "HMC: Caught std::string exception: " << e << std::endl;
888  QDP_abort(1);
889  }
890  catch(std::exception& e)
891  {
892  QDPIO::cerr << "HMC: Caught standard library exception: " << e.what() << std::endl;
893  QDP_abort(1);
894  }
895  catch(...)
896  {
897  QDPIO::cerr << "HMC: Caught generic/unknown exception" << std::endl;
898  QDP_abort(1);
899  }
900 
901  pop(xml_log); // hmc
902  pop(xml_out); // hmc
903 
904  snoop.stop();
905  QDPIO::cout << "HMC: total time = "
906  << snoop.getTimeInSeconds()
907  << " secs" << std::endl;
908 
909  END_CODE();
910 
912  exit(0);
913 }
914 
915 
916 // Repro check
917 namespace Chroma {
918 
919  bool
920  checkReproducability( const multi1d<LatticeColorMatrix>& P_new,
921  const multi1d<LatticeColorMatrix>& Q_new,
922  const QDP::Seed& seed_new,
923  const multi1d<LatticeColorMatrix>& P_old,
924  const multi1d<LatticeColorMatrix>& Q_old,
925  const QDP::Seed& seed_old )
926  {
927 
928  // Compare local contributions of P
929  int diffs_found = 0;
930  if ( P_new.size() != P_old.size() ) {
931  // Something bad happened if P_old and P_new are not the same
932  return false;
933  }
934 
935  if ( Q_new.size() != Q_old.size() ) {
936  // Something bad happened if Q_old and Q_new are not the same
937  return false;
938  }
939 
940  // Count the number of bytes to compare
941  int bytes=
942  2*Nc*Nc*Layout::sitesOnNode()*sizeof(WordType< LatticeColorMatrix >::Type_t);
943 
944  // Compare P_s
945  for(int mu=0; mu < P_new.size(); mu++) {
946  const unsigned char *p1 = (const unsigned char *)P_new[mu].getF();
947  const unsigned char *p2 = (const unsigned char *)P_old[mu].getF();
948  for(int b=0; b < bytes; b++) {
949  unsigned char diff = *p1 - *p2;
950  if( diff != 0 ) diffs_found++;
951  p1++; p2++;
952  }
953  }
954 
955  // Sum up the number of diffs found globally
956  QDPInternal::globalSum(diffs_found);
957 
958  if( diffs_found != 0 ) {
959  QDPIO::cout << "Found " << diffs_found << " different bytes in momentum repro check" << std::endl;
960  return false;
961  }
962 
963  diffs_found = 0;
964  for(int mu=0; mu < P_new.size(); mu++) {
965  const unsigned char *p1 = (const unsigned char *)Q_new[mu].getF();
966  const unsigned char *p2 = (const unsigned char *)Q_old[mu].getF();
967  for(int b=0; b < bytes; b++) {
968  unsigned char diff = *p1 - *p2;
969  if( diff != 0 ) diffs_found++;
970  p1++; p2++;
971  }
972  }
973 
974  // Sum up the number of diffs found globally
975  QDPInternal::globalSum(diffs_found);
976 
977  if( diffs_found != 0 ) {
978  QDPIO::cout << "Found " << diffs_found << " different bytes in gauge repro check" << std::endl;
979  return false;
980  }
981 
982  if( ! toBool( seed_new == seed_old ) ) {
983  QDPIO::cout << "New and old RNG seeds do not match " << std::endl;
984  return false;
985  }
986 
987  return true;
988  }
989 
990 }
Primary include file for CHROMA in application codes.
Abstract HMC trajectory.
Definition: abs_hmc.h:25
virtual unsigned long getFrequency(void) const =0
Tell me how often I should measure this beastie.
The Exact Hamiltonian Class - supplies internal field refreshment and energy calculations.
Pure gauge field state.
Definition: field_state.h:83
const multi1d< LatticeColorMatrix > & getP(void) const
Accessors.
Definition: field_state.h:117
const multi1d< LatticeColorMatrix > & getQ(void) const
Definition: field_state.h:118
Class for counted reference semantics.
Definition: handle.h:33
static T & Instance()
Definition: singleton.h:432
int main(int argc, char *argv[])
Definition: const_hmc.cc:696
int mu
Definition: cool.cc:24
void read(XMLReader &xml, const std::string &path, AsqtadFermActParams &param)
Read parameters.
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
void proginfo(XMLWriter &xml)
Print out basic information about this program.
Definition: proginfo.cc:24
void reset()
Reset the default gauge field state.
void set(const multi1d< LatticeColorMatrix > &u, XMLBufferWriter &record_xml)
Set the default gauge field.
void readNamedMonomialArray(XMLReader &xml, const std::string &path)
Read an array of named monomials from an XML reader. use factory to create the monomials and put them...
Definition: monomial_io.cc:106
void writeGauge(XMLBufferWriter &file_xml, XMLBufferWriter &record_xml, const multi1d< LatticeColorMatrix > &u, const std::string &file, QDP_volfmt_t volfmt, QDP_serialparallel_t serpar)
Write a gauge config in QIO format.
Definition: gauge_io.cc:60
GroupXML_t readXMLGroup(XMLReader &xml_in, const std::string &path, const std::string &type_name)
Read group and return as a std::string.
void setrn(int iseed[4])
Definition: make_seeds.cc:38
static int m[4]
Definition: make_seeds.cc:16
void savern(int iseed[4])
Definition: make_seeds.cc:46
Nd
Definition: meslate.cc:74
bool registerAll()
Register all the factories.
bool registerAll()
Register all the factories.
bool registerAll()
Register all the factories.
bool registerAll()
Register all the factories.
bool registerAll()
Register all the factories.
GroupXML_t createXMLGroup(const Params &p)
Returns a link smearing group with these params.
bool registerAll()
Register all the factories.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
XMLFileWriter & getXMLLogInstance()
Get xml log instance.
Definition: chroma_init.cc:378
static multi1d< LatticeColorMatrix > u
push(xml_out,"Condensates")
int i
Definition: pbg5p_w.cc:55
void saveState(const UpdateParams &update_params, MCControl &mc_control, unsigned long update_no, const multi1d< LatticeColorMatrix > &u)
Definition: const_hmc.cc:251
void initialize(int *argc, char ***argv)
Chroma initialisation routine.
Definition: chroma_init.cc:114
void finalize(void)
Chroma finalization routine.
Definition: chroma_init.cc:308
bool linkageHack(void)
Definition: const_hmc.cc:660
void reunit(LatticeColorMatrixF3 &xa)
Definition: reunit.cc:467
void doHMC(multi1d< LatticeColorMatrix > &u, AbsHMCTrj< multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > &theHMCTrj, MCControl &mc_control, const UpdateParams &update_params, multi1d< Handle< AbsInlineMeasurement > > &user_measurements)
Definition: const_hmc.cc:356
void setForceMonitoring(bool monitorP)
std::string getXMLInputFileName()
Get input file name.
Definition: chroma_init.cc:88
pop(xml_out)
START_CODE()
Complex b
Definition: invbicg.cc:96
bool checkReproducability(const multi1d< LatticeColorMatrix > &P_new, const multi1d< LatticeColorMatrix > &Q_new, const QDP::Seed &seed_new, const multi1d< LatticeColorMatrix > &P_old, const multi1d< LatticeColorMatrix > &Q_old, const QDP::Seed &seed_old)
Definition: const_hmc.cc:920
XMLFileWriter & getXMLOutputInstance()
Get xml output instance.
Definition: chroma_init.cc:359
::std::string string
Definition: gtest.h:1979
Real dummy
Definition: qtopcor.cc:36
Parameter structure for new Hamiltonian.
Hold group xml and type id.
std::string Integrator_xml
Definition: const_hmc.cc:184
multi1d< int > nrow
Definition: const_hmc.cc:179
std::string H_MC_xml
Definition: const_hmc.cc:183
std::string Monomials_xml
Definition: const_hmc.cc:182
A Structure to hold the top level parameters.
Params controlling running of monte carlo.
Definition: const_hmc.cc:14
unsigned long n_warm_up_updates
Definition: const_hmc.cc:18
std::string inline_measurement_xml
Definition: const_hmc.cc:25
unsigned long n_production_updates
Definition: const_hmc.cc:19
QDP_serialparallel_t save_pario
Definition: const_hmc.cc:24
QDP::Seed rng_seed
Definition: const_hmc.cc:16
unsigned int save_interval
Definition: const_hmc.cc:21
unsigned int n_updates_this_run
Definition: const_hmc.cc:20
QDP_volfmt_t save_volfmt
Definition: const_hmc.cc:23
GroupXML_t cfg
Definition: const_hmc.cc:15
int repro_check_frequency
Definition: const_hmc.cc:27
std::string save_prefix
Definition: const_hmc.cc:22
unsigned long start_update_num
Definition: const_hmc.cc:17
Params for initializing config.