CHROMA
inline_prop_distillution_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Compute propagators from distillution
3  *
4  * Propagator calculation in distillution
5  */
6 
7 #include "fermact.h"
10 #include "meas/glue/mesplq.h"
13 #include "qdp_map_obj.h"
14 #include "qdp_map_obj_disk.h"
15 #include "qdp_disk_map_slice.h"
17 #include "util/ferm/transf.h"
18 #include "util/ferm/spin_rep.h"
19 #include "util/ferm/diractodr.h"
21 #include "util/ft/time_slice_set.h"
22 #include "util/info/proginfo.h"
26 
28 
29 namespace Chroma
30 {
31  //----------------------------------------------------------------------------
32  namespace InlinePropDistillutionEnv
33  {
34  //! Propagator input
35  void read(XMLReader& xml, const std::string& path, Params::NamedObject_t& input)
36  {
37  XMLReader inputtop(xml, path);
38 
39  read(inputtop, "save_solnP", input.save_solnP);
40  read(inputtop, "gauge_id", input.gauge_id);
41  read(inputtop, "distillution_id", input.distillution_id);
42  read(inputtop, "src_file", input.src_file);
43  read(inputtop, "soln_file", input.soln_file);
44  }
45 
46  //! Propagator output
47  void write(XMLWriter& xml, const std::string& path, const Params::NamedObject_t& input)
48  {
49  push(xml, path);
50 
51  write(xml, "save_solnP", input.save_solnP);
52  write(xml, "gauge_id", input.gauge_id);
53  write(xml, "distillution_id", input.distillution_id);
54  write(xml, "src_file", input.src_file);
55  write(xml, "soln_file", input.soln_file);
56 
57  pop(xml);
58  }
59 
60 
61  //! Propagator input
62  void read(XMLReader& xml, const std::string& path, Params::Param_t::Contract_t& input)
63  {
64  XMLReader inputtop(xml, path);
65 
66  read(inputtop, "quark_lines", input.quark_lines);
67  read(inputtop, "mass", input.mass);
68 
69  // Read quark-line parameters
70  input.quark_line_xml = readXMLGroup(inputtop, "QuarkLine", "QuarkLineType");
71  }
72 
73  //! Propagator output
74  void write(XMLWriter& xml, const std::string& path, const Params::Param_t::Contract_t& input)
75  {
76  push(xml, path);
77 
78  write(xml, "quark_lines", input.quark_lines);
79  write(xml, "mass", input.mass);
80  xml << input.quark_line_xml.xml;
81 
82  pop(xml);
83  }
84 
85 
86  //! Propagator input
87  void read(XMLReader& xml, const std::string& path, Params::Param_t& input)
88  {
89  XMLReader inputtop(xml, path);
90 
91  read(inputtop, "Propagator", input.prop);
92  read(inputtop, "Contractions", input.contract);
93  }
94 
95  //! Propagator output
96  void write(XMLWriter& xml, const std::string& path, const Params::Param_t& input)
97  {
98  push(xml, path);
99 
100  write(xml, "Propagator", input.prop);
101  write(xml, "Contractions", input.contract);
102 
103  pop(xml);
104  }
105 
106 
107  //! Propagator input
108  void read(XMLReader& xml, const std::string& path, Params& input)
109  {
110  Params tmp(xml, path);
111  input = tmp;
112  }
113 
114  //! Propagator output
115  void write(XMLWriter& xml, const std::string& path, const Params& input)
116  {
117  push(xml, path);
118 
119  write(xml, "Param", input.param);
120  write(xml, "NamedObject", input.named_obj);
121 
122  pop(xml);
123  }
124  } // namespace InlinePropDistillutionEnv
125 
126 
127  //----------------------------------------------------------------------------
128  namespace InlinePropDistillutionEnv
129  {
130  namespace
131  {
132  AbsInlineMeasurement* createMeasurement(XMLReader& xml_in,
133  const std::string& path)
134  {
135  return new InlineMeas(Params(xml_in, path));
136  }
137 
138  //! Local registration flag
139  bool registered = false;
140  }
141 
142  const std::string name = "PROP_DISTILLUTION";
143 
144  //! Register all the factories
145  bool registerAll()
146  {
147  bool success = true;
148  if (! registered)
149  {
152  success &= TheInlineMeasurementFactory::Instance().registerObject(name, createMeasurement);
153  registered = true;
154  }
155  return success;
156  }
157 
158 
159  //----------------------------------------------------------------------------
160  // Param stuff
162 
163  Params::Params(XMLReader& xml_in, const std::string& path)
164  {
165  try
166  {
167  XMLReader paramtop(xml_in, path);
168 
169  if (paramtop.count("Frequency") == 1)
170  read(paramtop, "Frequency", frequency);
171  else
172  frequency = 1;
173 
174  // Parameters for source construction
175  read(paramtop, "Param", param);
176 
177  // Read in the output propagator/source configuration info
178  read(paramtop, "NamedObject", named_obj);
179 
180  // Possible alternate XML file pattern
181  if (paramtop.count("xml_file") != 0)
182  {
183  read(paramtop, "xml_file", xml_file);
184  }
185  }
186  catch(const std::string& e)
187  {
188  QDPIO::cerr << __func__ << ": Caught Exception reading XML: " << e << std::endl;
189  QDP_abort(1);
190  }
191  }
192 
193 
194 
195  // Function call
196  void
197  InlineMeas::operator()(unsigned long update_no,
198  XMLWriter& xml_out)
199  {
200  // If xml file not empty, then use alternate
201  if (params.xml_file != "")
202  {
203  std::string xml_file = makeXMLFileName(params.xml_file, update_no);
204 
205  push(xml_out, "PropDistillution");
206  write(xml_out, "update_no", update_no);
207  write(xml_out, "xml_file", xml_file);
208  pop(xml_out);
209 
210  XMLFileWriter xml(xml_file);
211  func(update_no, xml);
212  }
213  else
214  {
215  func(update_no, xml_out);
216  }
217  }
218 
219 
220  // Real work done here
221  void
222  InlineMeas::func(unsigned long update_no,
223  XMLWriter& xml_out)
224  {
225  START_CODE();
226 
227  StopWatch snoop;
228  snoop.reset();
229  snoop.start();
230 
231  // Test and grab a reference to the gauge field
232  multi1d<LatticeColorMatrix> u;
233  XMLBufferWriter gauge_xml;
234  try
235  {
236  u = TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
237  TheNamedObjMap::Instance().get(params.named_obj.gauge_id).getRecordXML(gauge_xml);
238  }
239  catch( std::bad_cast )
240  {
241  QDPIO::cerr << name << ": caught dynamic cast error" << std::endl;
242  QDP_abort(1);
243  }
244  catch (const std::string& e)
245  {
246  QDPIO::cerr << name << ": std::map call failed: " << e << std::endl;
247  QDP_abort(1);
248  }
249 
250  push(xml_out, "PropDistillution");
251  write(xml_out, "update_no", update_no);
252 
253  QDPIO::cout << name << ": propagator calculation" << std::endl;
254 
255  proginfo(xml_out); // Print out basic program info
256 
257  // Write out the input
258  write(xml_out, "Input", params);
259 
260  // Write out the config header
261  write(xml_out, "Config_info", gauge_xml);
262 
263  push(xml_out, "Output_version");
264  write(xml_out, "out_version", 1);
265  pop(xml_out);
266 
267  // Calculate some gauge invariant observables just for info.
268  MesPlq(xml_out, "Observables", u);
269 
270  //
271  // Read in the source along with relevant information.
272  //
273  QDPIO::cout << "Snarf the distillution factory from a named buffer" << std::endl;
274  try
275  {
277  }
278  catch (std::bad_cast) {
279  QDPIO::cerr << name << ": caught dynamic cast error" << std::endl;
280  QDP_abort(1);
281  }
282  catch (const std::string& e) {
283  QDPIO::cerr << name << ": error extracting source_header: " << e << std::endl;
284  QDP_abort(1);
285  }
286  catch( const char* e) {
287  QDPIO::cerr << name <<": Caught some char* exception:" << std::endl;
288  QDPIO::cerr << e << std::endl;
289  QDPIO::cerr << "Rethrowing" << std::endl;
290  throw;
291  }
292 
293  // Cast should be valid now
294  const DistillutionNoise& dist_noise_obj =
296 
297  // Some diagnostics
298  QDPIO::cout << "Distillution factory: ensemble= XX" << dist_noise_obj.getEnsemble() << "XX "
299  << "sequence= XX" << dist_noise_obj.getSequence() << "XX\n"
300  << "t_origin= " << dist_noise_obj.getOrigin() << std::endl;
301 
302  // Will use TimeSliceSet-s a lot
303  const int decay_dir = dist_noise_obj.getDecayDir();
304  const int Lt = Layout::lattSize()[decay_dir];
305 
306  // A sanity check
307  if (decay_dir != Nd-1)
308  {
309  QDPIO::cerr << __func__ << ": TimeSliceIO only supports decay_dir= " << Nd-1 << "\n";
310  QDP_abort(1);
311  }
312 
313  // The time-slice set
314  TimeSliceSet time_slice_set(decay_dir);
315 
316 
317  //
318  // Map-object-disk storage of the source file
319  //
320  QDP::MapObjectDisk<KeyPropDistillution_t, TimeSliceIO<LatticeColorVectorF> > source_obj;
321  source_obj.setDebug(0);
322 
323  QDPIO::cout << "Open source file" << std::endl;
324 
325  if (! source_obj.fileExists(params.named_obj.src_file))
326  {
327  QDPIO::cerr << name << ": source file does not exist: src_file= " << params.named_obj.src_file << std::endl;
328  QDP_abort(1);
329  }
330  else
331  {
332  source_obj.open(params.named_obj.src_file, std::ios_base::in);
333  }
334 
335  QDPIO::cout << "Finished opening solution file" << std::endl;
336 
337 
338  //
339  // Map-object-disk storage
340  //
341  QDP::MapObjectDisk<KeyPropDistillution_t, TimeSliceIO<LatticeColorVectorF> > prop_obj;
342  prop_obj.setDebug(0);
343 
345  {
346  QDPIO::cout << "Open solution file" << std::endl;
347 
348  if (! prop_obj.fileExists(params.named_obj.soln_file))
349  {
350  XMLBufferWriter file_xml;
351 
352  push(file_xml, "MODMetaData");
353  write(file_xml, "id", std::string("propDist"));
354  write(file_xml, "lattSize", QDP::Layout::lattSize());
355  write(file_xml, "decay_dir", decay_dir);
356  write(file_xml, "ensemble", dist_noise_obj.getEnsemble());
357  write(file_xml, "sequence", dist_noise_obj.getSequence());
358  write(file_xml, "t_origin", dist_noise_obj.getOrigin());
359  file_xml << params.param.contract.quark_line_xml.xml;
360  write(file_xml, "quark_line", params.param.contract.quark_lines[0]);
361  proginfo(file_xml); // Print out basic program info
362  write(file_xml, "Params", params.param);
363  write(file_xml, "Config_info", gauge_xml);
364  pop(file_xml);
365 
366  std::string file_str(file_xml.str());
367 
368  prop_obj.insertUserdata(file_xml.str());
369  prop_obj.open(params.named_obj.soln_file, std::ios_base::in | std::ios_base::out | std::ios_base::trunc);
370  }
371  else
372  {
373  prop_obj.open(params.named_obj.soln_file);
374  }
375 
376  QDPIO::cout << "Finished opening solution file" << std::endl;
377  }
378 
379 
380  // Total number of iterations
381  int ncg_had = 0;
382 
383 
384  // Rotation from DR to DP
385  SpinMatrix diracToDRMat(DiracToDRMat());
386  std::vector<MatrixSpinRep_t> diracToDrMatPlus = convertTwoQuarkSpin(diracToDRMat);
387  std::vector<MatrixSpinRep_t> diracToDrMatMinus = convertTwoQuarkSpin(adj(diracToDRMat));
388 
389 
390  //
391  // Try the factories
392  //
393  try
394  {
395  StopWatch swatch;
396  swatch.reset();
397  QDPIO::cout << "Try the various factories" << std::endl;
398 
399  // Typedefs to save typing
400  typedef LatticeFermion T;
401  typedef multi1d<LatticeColorMatrix> P;
402  typedef multi1d<LatticeColorMatrix> Q;
403 
404  //
405  // Initialize fermion action
406  //
407  std::istringstream xml_s(params.param.prop.fermact.xml);
408  XMLReader fermacttop(xml_s);
409  QDPIO::cout << "FermAct = " << params.param.prop.fermact.id << std::endl;
410 
411  // Generic Wilson-Type stuff
414  fermacttop,
416 
417  Handle< FermState<T,P,Q> > state(S_f->createState(u));
418 
421 
422  QDPIO::cout << "Suitable factory found: compute all the quark props" << std::endl;
423  swatch.start();
424 
425 
426  //
427  // Loop over the source color and spin, creating the source
428  // and calling the relevant propagator routines.
429  //
430  // Loop over each quark-line
431  for(std::vector<int>::const_iterator quark_line= params.param.contract.quark_lines.begin();
432  quark_line != params.param.contract.quark_lines.end();
433  ++quark_line)
434  {
435  QDPIO::cout << "quark_line = " << *quark_line << " type= " << params.param.contract.quark_line_xml.id << std::endl;
436 
437  //
438  // Factory for quark line
439  //
440  Handle<AbsQuarkLine> quark_line_fact;
441 
442  try
443  {
444  std::istringstream xml_l(params.param.contract.quark_line_xml.xml);
445  XMLReader linktop(xml_l);
446  QDPIO::cout << "Quark line type = " << params.param.contract.quark_line_xml.id << std::endl;
447  QDPIO::cout << "Quark line xml = XX" << params.param.contract.quark_line_xml.xml << "XX" << std::endl;
448 
449  QDPIO::cout << "Create quark-line factory" << std::endl;
450 
451  quark_line_fact =
453  linktop,
455  dist_noise_obj,
456  source_obj,
458  *quark_line,
460 
461  QDPIO::cout << "Factory created" << std::endl;
462  }
463  catch(const std::string& e)
464  {
465  QDPIO::cerr << InlinePropDistillutionEnv::name << ": error creating quark-line factory: " << e << std::endl;
466  QDP_abort(1);
467  }
468  catch(...)
469  {
470  QDPIO::cerr << InlinePropDistillutionEnv::name << ": generic exception in creating quark-line factory" << std::endl;
471  QDP_abort(1);
472  }
473 
474 
475  // Loop over source locations
476  std::vector<int> t_sources(quark_line_fact->getTimeSources());
477 
478  for(int tt=0; tt < t_sources.size(); ++tt)
479  {
480  int t_source = t_sources[tt]; // This is the pretend time-slice. The actual value is shifted.
481  QDPIO::cout << "t_source = " << t_source << std::endl;
482 
483  // The space distillution loop
484  for(int dist_src=0; dist_src < quark_line_fact->getNumSpaceDils(); ++dist_src)
485  {
486  StopWatch sniss1;
487  sniss1.reset();
488  sniss1.start();
489  QDPIO::cout << "dist_src = " << dist_src << std::endl;
490 
491  // Prepare a distilluted source
492  LatticeColorVector vec_srce = quark_line_fact->getSrc(t_source, dist_src);
493 
494  //
495  // Loop over each spin source and invert.
496  // Use the same colorstd::vector source. No spin dilution will be used.
497  //
498  multi2d<LatticeColorVector> ferm_out(Ns,Ns);
499 
500  for(int spin_source=0; spin_source < Ns; ++spin_source)
501  {
502  QDPIO::cout << "spin_source = " << spin_source << std::endl;
503 
504  // Insert a ColorVector into spin index spin_source
505  // This only overwrites sections, so need to initialize first
506  LatticeFermion chi = zero;
507  CvToFerm(vec_srce, chi, spin_source);
508 
509  LatticeFermion quark_soln = zero;
510 
511  // Do the propagator inversion
512  SystemSolverResults_t res = (*PP)(quark_soln, chi);
513  ncg_had = res.n_count;
514 
515  // Extract into the temporary output array
516  for(int spin_sink=0; spin_sink < Ns; ++spin_sink)
517  {
518  ferm_out(spin_sink,spin_source) = peekSpin(quark_soln, spin_sink);
519  }
520  } // for spin_source
521 
522 
523  // Rotate from DeGrand-Rossi (DR) to Dirac-Pauli (DP)
524  {
525  multi2d<LatticeColorVector> ferm_tmp;
526 
527  multiplyRep(ferm_tmp, diracToDrMatMinus, ferm_out);
528  multiplyRep(ferm_out, ferm_tmp, diracToDrMatPlus);
529  }
530 
531  sniss1.stop();
532  QDPIO::cout << "Time to assemble and transmogrify propagators for dist_src= " << dist_src << " time = "
533  << sniss1.getTimeInSeconds()
534  << " secs" << std::endl;
535 
536 
537  // Write out each time-slice chunk of a lattice colorvec soln to disk
538  QDPIO::cout << "Potentially write propagator solutions to disk" << std::endl;
539  StopWatch sniss2;
540  sniss2.reset();
541  sniss2.start();
542 
543  // Write the solutions
545  {
546  QDPIO::cout << "Write propagator solution to disk" << std::endl;
547  std::list<KeyPropDistillution_t> snk_keys(quark_line_fact->getSnkKeys(t_source, dist_src));
548 
549  for(std::list<KeyPropDistillution_t>::const_iterator key= snk_keys.begin();
550  key != snk_keys.end();
551  ++key)
552  {
553  LatticeColorVectorF tmptmp = ferm_out(key->spin_snk,key->spin_src);
554 
555  prop_obj.insert(*key, TimeSliceIO<LatticeColorVectorF>(tmptmp, dist_noise_obj.getTime(key->t_slice)));
556  } // for key
557  }
558 
559  sniss2.stop();
560  QDPIO::cout << "Time to write propagators for dist_src= " << dist_src << " time = "
561  << sniss2.getTimeInSeconds()
562  << " secs" << std::endl;
563 
564  } // for dist_src
565  } // for tt
566  } // for quark_line
567 
568  swatch.stop();
569  QDPIO::cout << "Propagators computed: time= "
570  << swatch.getTimeInSeconds()
571  << " secs" << std::endl;
572  }
573  catch (const std::string& e)
574  {
575  QDPIO::cout << name << ": caught exception around qprop: " << e << std::endl;
576  QDP_abort(1);
577  }
578 
579  push(xml_out,"Relaxation_Iterations");
580  write(xml_out, "ncg_had", ncg_had);
581  pop(xml_out);
582 
583  pop(xml_out); // prop_dist
584 
585  snoop.stop();
586  QDPIO::cout << name << ": total time = "
587  << snoop.getTimeInSeconds()
588  << " secs" << std::endl;
589 
590  QDPIO::cout << name << ": ran successfully" << std::endl;
591 
592  END_CODE();
593  }
594 
595  }
596 
597 } // namespace Chroma
Inline measurement factory.
virtual std::string getSequence() const
Return the sequence.
virtual int getTime(int t_slice) const
Convenience - get shifted time.
virtual int getDecayDir() const
Return the decay direction.
virtual std::string getEnsemble() const
Return the ensemble.
virtual int getOrigin() const
Return the actual time origin.
Class for counted reference semantics.
Definition: handle.h:33
void func(const unsigned long update_no, XMLWriter &xml_out)
Do the measurement.
void operator()(const unsigned long update_no, XMLWriter &xml_out)
Do the measurement.
static T & Instance()
Definition: singleton.h:432
Builds time slice subsets.
Basis rotation matrix from Dirac to Degrand-Rossi (and reverse)
Distillution factory for producing keys * sources.
Support for distillution - random time-slices and quark line noises.
Class structure for fermion actions.
Fermion action factories.
All Wilson-type fermion actions.
SpinMatrixD DiracToDRMat()
The Dirac to Degrand-Rossi spin transformation matrix (and reverse)
Definition: diractodr.cc:22
void CvToFerm(const LatticeColorVectorF &a, LatticeFermionF &b, int spin_index)
Convert (insert) a LatticeColorVector into a LatticeFermion.
Definition: transf.cc:18
void proginfo(XMLWriter &xml)
Print out basic information about this program.
Definition: proginfo.cc:24
std::string makeXMLFileName(std::string xml_file, unsigned long update_no)
Return a xml file name for inline measurements.
GroupXML_t readXMLGroup(XMLReader &xml_in, const std::string &path, const std::string &type_name)
Read group and return as a std::string.
TimeSliceSet time_slice_set
Compute the propagator from distillution.
Key for distillution propagator sources and solutions.
Make xml file writer.
Nd
Definition: meslate.cc:74
Named object function std::map.
static bool registered
Local registration flag.
bool registerAll()
Register all the factories.
multi1d< LatticeColorMatrix > P
void write(XMLWriter &xml, const std::string &path, const Params::NamedObject_t &input)
Propagator output.
void read(XMLReader &xml, const std::string &path, Params::NamedObject_t &input)
Propagator input.
bool registerAll()
Register all the factories.
bool registerAll()
Register all the factories.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
static multi1d< LatticeColorMatrix > u
LinOpSysSolverMGProtoClover::Q Q
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
push(xml_out,"Condensates")
LinOpSysSolverMGProtoClover::T T
std::vector< MatrixSpinRep_t > convertTwoQuarkSpin(const SpinMatrix &in)
Convert a generic spin matrix into a sparse spin representation.
Definition: spin_rep.cc:67
multi1d< LatticeFermion > chi(Ncb)
pop(xml_out)
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
static QDP_ColorVector * out
Constructor.
START_CODE()
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
const WilsonTypeFermAct< multi1d< LatticeFermion > > & S_f
Definition: pbg5p_w.cc:27
Double zero
Definition: invbicg.cc:106
void multiplyRep(multi2d< LatticeColorVector > &dist_rep, const std::vector< MatrixSpinRep_t > &spin, const multi2d< LatticeColorVector > &prop_rep)
Dist(t2) = SpinMatrix*Prop(t2)
static QDP_ColorVector * in
::std::string string
Definition: gtest.h:1979
Print out basic info about this program.
Sparse matrix representation of spin matrices.
GroupXML_t invParam
Definition: qprop_io.h:84
GroupXML_t fermact
Definition: qprop_io.h:80
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Convenience for building time-slice subsets.
Contraction operators for two quarks.