CHROMA
inline_annih_prop_matelem_colorvec_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Compute the annihilation diagram propagator elements M^-1 * multi1d<LatticeColorVector>
3  *
4  * Annihilation diagrams version of propagator calculation on a colorstd::vector
5  */
6 
7 #include "fermact.h"
10 #include "meas/glue/mesplq.h"
11 #include "meas/sources/zN_src.h"
13 #include "qdp_map_obj_memory.h"
14 #include "util/ferm/key_val_db.h"
17 #include "util/ferm/transf.h"
18 #include "util/ft/sftmom.h"
19 #include "util/info/proginfo.h"
23 
25 
26 namespace Chroma
27 {
28  namespace InlineAnnihPropMatElemColorVecEnv
29  {
30  //! Propagator input
31  void read(XMLReader& xml, const std::string& path, Params::NamedObject_t& input)
32  {
33  XMLReader inputtop(xml, path);
34 
35  read(inputtop, "gauge_id", input.gauge_id);
36  read(inputtop, "colorvec_id", input.colorvec_id);
37  read(inputtop, "prop_op_file", input.prop_op_file);
38  }
39 
40  //! Propagator output
41  void write(XMLWriter& xml, const std::string& path, const Params::NamedObject_t& input)
42  {
43  push(xml, path);
44 
45  write(xml, "gauge_id", input.gauge_id);
46  write(xml, "colorvec_id", input.colorvec_id);
47  write(xml, "prop_op_file", input.prop_op_file);
48 
49  pop(xml);
50  }
51 
52 
53  //! Propagator input
54  void read(XMLReader& xml, const std::string& path, Params::Param_t::Contract_t& input)
55  {
56  XMLReader inputtop(xml, path);
57 
58  read(inputtop, "num_vecs", input.num_vecs);
59  read(inputtop, "t_sources_start", input.t_sources_start);
60  read(inputtop, "dt", input.dt);
61  read(inputtop, "decay_dir", input.decay_dir);
62  read(inputtop, "N", input.N);
63  read(inputtop, "ran_seed", input.ran_seed);
64  read(inputtop, "mass_label", input.mass_label);
65  }
66 
67  //! Propagator output
68  void write(XMLWriter& xml, const std::string& path, const Params::Param_t::Contract_t& input)
69  {
70  push(xml, path);
71 
72  write(xml, "num_vecs", input.num_vecs);
73  write(xml, "t_sources_start", input.t_sources_start);
74  write(xml, "dt", input.dt);
75  write(xml, "decay_dir", input.decay_dir);
76  write(xml, "N", input.N);
77  write(xml, "ran_seed", input.ran_seed);
78  write(xml, "mass_label", input.mass_label);
79 
80  pop(xml);
81  }
82 
83 
84  //! Propagator input
85  void read(XMLReader& xml, const std::string& path, Params::Param_t& input)
86  {
87  XMLReader inputtop(xml, path);
88 
89  read(inputtop, "Propagator", input.prop);
90  read(inputtop, "Contractions", input.contract);
91  }
92 
93  //! Propagator output
94  void write(XMLWriter& xml, const std::string& path, const Params::Param_t& input)
95  {
96  push(xml, path);
97 
98  write(xml, "Propagator", input.prop);
99  write(xml, "Contractions", input.contract);
100 
101  pop(xml);
102  }
103 
104 
105  //! Propagator input
106  void read(XMLReader& xml, const std::string& path, Params& input)
107  {
108  Params tmp(xml, path);
109  input = tmp;
110  }
111 
112  //! Propagator output
113  void write(XMLWriter& xml, const std::string& path, const Params& input)
114  {
115  push(xml, path);
116 
117  write(xml, "Param", input.param);
118  write(xml, "NamedObject", input.named_obj);
119 
120  pop(xml);
121  }
122  } // namespace InlinePropColorVecEnv
123 
124 
125  namespace InlineAnnihPropMatElemColorVecEnv
126  {
127  namespace
128  {
129  AbsInlineMeasurement* createMeasurement(XMLReader& xml_in,
130  const std::string& path)
131  {
132  return new InlineMeas(Params(xml_in, path));
133  }
134 
135  //! Local registration flag
136  bool registered = false;
137  }
138 
139  const std::string name = "ANNIH_PROP_MATELEM_COLORVEC";
140 
141  //! Register all the factories
142  bool registerAll()
143  {
144  bool success = true;
145  if (! registered)
146  {
148  success &= TheInlineMeasurementFactory::Instance().registerObject(name, createMeasurement);
149  registered = true;
150  }
151  return success;
152  }
153 
154 
155  //----------------------------------------------------------------------------
156  // Param stuff
158 
159  Params::Params(XMLReader& xml_in, const std::string& path)
160  {
161  try
162  {
163  XMLReader paramtop(xml_in, path);
164 
165  if (paramtop.count("Frequency") == 1)
166  read(paramtop, "Frequency", frequency);
167  else
168  frequency = 1;
169 
170  // Parameters for source construction
171  read(paramtop, "Param", param);
172 
173  // Read in the output propagator/source configuration info
174  read(paramtop, "NamedObject", named_obj);
175 
176  // Possible alternate XML file pattern
177  if (paramtop.count("xml_file") != 0)
178  {
179  read(paramtop, "xml_file", xml_file);
180  }
181  }
182  catch(const std::string& e)
183  {
184  QDPIO::cerr << __func__ << ": Caught Exception reading XML: " << e << std::endl;
185  QDP_abort(1);
186  }
187  }
188 
189 
190 
191  // Function call
192  void
193  InlineMeas::operator()(unsigned long update_no,
194  XMLWriter& xml_out)
195  {
196  // If xml file not empty, then use alternate
197  if (params.xml_file != "")
198  {
199  std::string xml_file = makeXMLFileName(params.xml_file, update_no);
200 
201  push(xml_out, "AnnihPropMatElemColorVec");
202  write(xml_out, "update_no", update_no);
203  write(xml_out, "xml_file", xml_file);
204  pop(xml_out);
205 
206  XMLFileWriter xml(xml_file);
207  func(update_no, xml);
208  }
209  else
210  {
211  func(update_no, xml_out);
212  }
213  }
214 
215 
216  // Real work done here
217  void
218  InlineMeas::func(unsigned long update_no,
219  XMLWriter& xml_out)
220  {
221  START_CODE();
222 
223  StopWatch snoop;
224  snoop.reset();
225  snoop.start();
226 
227  // Test and grab a reference to the gauge field
228  multi1d<LatticeColorMatrix> u;
229  XMLBufferWriter gauge_xml;
230  try
231  {
232  u = TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
233  TheNamedObjMap::Instance().get(params.named_obj.gauge_id).getRecordXML(gauge_xml);
234  }
235  catch( std::bad_cast )
236  {
237  QDPIO::cerr << name << ": caught dynamic cast error" << std::endl;
238  QDP_abort(1);
239  }
240  catch (const std::string& e)
241  {
242  QDPIO::cerr << name << ": std::map call failed: " << e << std::endl;
243  QDP_abort(1);
244  }
245 
246  push(xml_out, "AnnihPropMatElemColorVec");
247  write(xml_out, "update_no", update_no);
248 
249  QDPIO::cout << name << ": annihilation propagator matrix element calculation" << std::endl;
250 
251  proginfo(xml_out); // Print out basic program info
252 
253  // Write out the input
254  write(xml_out, "Input", params);
255 
256  // Write out the config header
257  write(xml_out, "Config_info", gauge_xml);
258 
259  push(xml_out, "Output_version");
260  write(xml_out, "out_version", 1);
261  pop(xml_out);
262 
263  // Calculate some gauge invariant observables just for info.
264  MesPlq(xml_out, "Observables", u);
265 
266  //
267  // Read in the source along with relevant information.
268  //
269  XMLReader source_file_xml, source_record_xml;
270 
271  QDPIO::cout << "Snarf the source from a named buffer" << std::endl;
272  try
273  {
275 
276  // Snarf the source info. This is will throw if the colorvec_id is not there
277  TheNamedObjMap::Instance().get(params.named_obj.colorvec_id).getFileXML(source_file_xml);
278  TheNamedObjMap::Instance().get(params.named_obj.colorvec_id).getRecordXML(source_record_xml);
279 
280  // Write out the source header
281  write(xml_out, "Source_file_info", source_file_xml);
282  write(xml_out, "Source_record_info", source_record_xml);
283  }
284  catch (std::bad_cast)
285  {
286  QDPIO::cerr << name << ": caught dynamic cast error" << std::endl;
287  QDP_abort(1);
288  }
289  catch (const std::string& e)
290  {
291  QDPIO::cerr << name << ": error extracting source_header: " << e << std::endl;
292  QDP_abort(1);
293  }
294 
295  // Cast should be valid now
296  const QDP::MapObject<int,EVPair<LatticeColorVector> >& eigen_source =
298 
299  QDPIO::cout << "Source successfully read and parsed" << std::endl;
300 
301  // Another sanity check
303  {
304  QDPIO::cerr << __func__ << ": num_vecs= " << params.param.contract.num_vecs
305  << " is greater than the number of available colorvectors= "
306  << eigen_source.size() << std::endl;
307  QDP_abort(1);
308  }
309 
310  // Interval of sources should divide into lattice extent
311  if ((QDP::Layout::lattSize()[params.param.contract.decay_dir] % params.param.contract.dt) != 0)
312  {
313  QDPIO::cerr << __func__ << ": dt= " << params.param.contract.dt
314  << " does not divide into lattice extent" << std::endl;
315  QDP_abort(1);
316  }
317 
318 
319  //
320  // DB storage
321  //
322  BinaryStoreDB< SerialDBKey<KeyPropElementalOperator_t>, SerialDBData<ValPropElementalOperator_t> > qdp_db;
323 
324  // Open the file, and write the meta-data and the binary for this operator
325  if (! qdp_db.fileExists(params.named_obj.prop_op_file))
326  {
327  XMLBufferWriter file_xml;
328 
329  push(file_xml, "DBMetaData");
330  write(file_xml, "id", std::string("propElemOp"));
331  write(file_xml, "lattSize", QDP::Layout::lattSize());
332  write(file_xml, "decay_dir", params.param.contract.decay_dir);
333  proginfo(file_xml); // Print out basic program info
334  write(file_xml, "Params", params.param.contract);
335  write(file_xml, "Config_info", gauge_xml);
337  pop(file_xml);
338 
339  std::string file_str(file_xml.str());
340  qdp_db.setMaxUserInfoLen(file_str.size());
341 
342  qdp_db.open(params.named_obj.prop_op_file, O_RDWR | O_CREAT, 0664);
343 
344  qdp_db.insertUserdata(file_str);
345  }
346  else
347  {
348  qdp_db.open(params.named_obj.prop_op_file, O_RDWR, 0664);
349  }
350 
351 
352 
353  // Save current seed
354  Seed ran_seed;
355  QDP::RNG::savern(ran_seed);
356 
357  // Set the seed to desired value
359 
360  // Total number of iterations
361  int ncg_had = 0;
362 
363  //
364  // Try the factories
365  //
366  try
367  {
368  StopWatch swatch;
369  swatch.reset();
370  QDPIO::cout << "Try the various factories" << std::endl;
371 
372  // Typedefs to save typing
373  typedef LatticeFermion T;
374  typedef multi1d<LatticeColorMatrix> P;
375  typedef multi1d<LatticeColorMatrix> Q;
376 
377  //
378  // Initialize fermion action
379  //
380  std::istringstream xml_s(params.param.prop.fermact.xml);
381  XMLReader fermacttop(xml_s);
382  QDPIO::cout << "FermAct = " << params.param.prop.fermact.id << std::endl;
383 
384  // Generic Wilson-Type stuff
387  fermacttop,
389 
390  Handle< FermState<T,P,Q> > state(S_f->createState(u));
391 
394 
395  QDPIO::cout << "Suitable factory found: compute all the quark props" << std::endl;
396  swatch.start();
397 
398  //
399  // Loop over the source color and spin, creating the source
400  // and calling the relevant propagator routines.
401  //
402  const int num_vecs = params.param.contract.num_vecs;
403  const int decay_dir = params.param.contract.decay_dir;
404  const int dt = params.param.contract.dt;
405  const int Lt = QDP::Layout::lattSize()[decay_dir];
406  const int num_sources = Lt / params.param.contract.dt;
407  const multi1d<int>& t_sources_start = params.param.contract.t_sources_start;
408 
409  // Initialize the slow Fourier transform phases
410  SftMom phases(0, true, decay_dir);
411 
412  for(int tt=0; tt < t_sources_start.size(); ++tt)
413  {
414  int t_source_start = params.param.contract.t_sources_start[tt];
415  // Construct all the solution vectors for only the first time source
416  QDPIO::cout << "t_source_start = " << t_source_start << std::endl;
417 
418  //
419  // Have to hold temporarily the solution vectors
420  //
421  MapObjectMemory<KeyPropColorVec_t,LatticeFermion> map_obj;
422 
423  // The random numbers used for the stochastic combination of sources is held here
424  MapObjectMemory<int,Complex> rng_map_obj;
425 
426  // Build up list of time sources
427  multi1d<int> t_sources(num_sources);
428 
429  for(int nn=0; nn < t_sources.size(); ++nn)
430  {
431  int t = (t_source_start + dt*nn + Lt) % Lt;
432  t_sources[nn] = t;
433 
434  rng_map_obj.insert(t, zN_rng(params.param.contract.N));
435  }
436 
437  for(int colorvec_source=0; colorvec_source < num_vecs; ++colorvec_source) {
438 
439  QDPIO::cout << "colorvec_source = " << colorvec_source << std::endl;
440 
441  // Accumulate the source from several time-slices
442  LatticeColorVector vec_srce = zero;
443  for(int nn=0; nn < t_sources.size(); ++nn) {
444  int t = t_sources[nn];
445 
446  // Pull out a time-slice of the color std::vector source, give it a random weight
447 
448  Complex weight; rng_map_obj.get(t, weight);
449  EVPair<LatticeColorVector> tmpvec; eigen_source.get(colorvec_source, tmpvec);
450  vec_srce[phases.getSet()[t]] = weight * tmpvec.eigenVector;
451  }
452 
453  for(int spin_source=0; spin_source < Ns; ++spin_source) {
454 
455  QDPIO::cout << "spin_source = " << spin_source << std::endl;
456 
457  // Insert a ColorVector into spin index spin_source
458  // This only overwrites sections, so need to initialize first
459  LatticeFermion chi = zero;
460  CvToFerm(vec_srce, chi, spin_source);
461 
462  LatticeFermion quark_soln = zero;
463 
464  // Do the propagator inversion
465  SystemSolverResults_t res = (*PP)(quark_soln, chi);
466  ncg_had = res.n_count;
467 
468  KeyPropColorVec_t key;
469  key.t_source = t_source_start;
470  key.colorvec_src = colorvec_source;
471  key.spin_src = spin_source;
472 
473  map_obj.insert(key, quark_soln);
474  } // for spin_source
475  } // for colorvec_source
476 
477 
478  swatch.stop();
479  QDPIO::cout << "Propagators computed: time= "
480  << swatch.getTimeInSeconds()
481  << " secs" << std::endl;
482 
483 
484  //
485  // All the loops
486  //
487  // NOTE: I pull the spin source and sink loops to the outside intentionally.
488  // The idea is to create a colorstd::vector index (2d) array. These are not
489  // too big, but are big enough to make the IO efficient, and the DB efficient
490  // on reading. For N=32 and Lt=128, the mats are 2MB.
491  //
492  swatch.reset();
493  swatch.start();
494  QDPIO::cout << "Extract matrix elements" << std::endl;
495 
496  for(int spin_source=0; spin_source < Ns; ++spin_source) {
497 
498  QDPIO::cout << "spin_source = " << spin_source << std::endl;
499 
500  for(int spin_sink=0; spin_sink < Ns; ++spin_sink) {
501 
502  QDPIO::cout << "spin_sink = " << spin_sink << std::endl;
503 
504  // Construct the keys and values. Have to hold the buffers here given that
505  // the source and sink spin indices are pulled out as well as the time
506  multi1d<KeyValPropElementalOperator_t> buf(t_sources.size());
507  for(int nn=0; nn < t_sources.size(); ++nn)
508  {
509  int t = t_sources[nn];
510 
511  buf[nn].key.key().t_slice = t;
512  buf[nn].key.key().t_source = t;
513  buf[nn].key.key().spin_src = spin_source;
514  buf[nn].key.key().spin_snk = spin_sink;
515  buf[nn].key.key().mass_label = params.param.contract.mass_label;
516  buf[nn].val.data().mat.resize(num_vecs,num_vecs);
517  }
518 
519  for(int colorvec_source=0; colorvec_source < num_vecs; ++colorvec_source)
520  {
521  KeyPropColorVec_t key;
522  key.t_source = t_source_start;
523  key.colorvec_src = colorvec_source;
524  key.spin_src = spin_source;
525 
526  LatticeColorVector vec_source;
527  {
528  LatticeFermion tmp;
529  QDPIO::cout << "MAP_OBJ LOOKUP" << std::endl << std::flush;
530  map_obj.get(key,tmp);
531  vec_source=peekSpin(tmp, spin_sink);
532  }
533 
534  for(int colorvec_sink=0; colorvec_sink < num_vecs; ++colorvec_sink)
535  {
536  EVPair<LatticeColorVector> vec_sink; eigen_source.get(colorvec_sink, vec_sink);
537 
538  multi1d<ComplexD> hsum(sumMulti(localInnerProduct(vec_sink.eigenVector, vec_source), phases.getSet()));
539 
540  for(int nn=0; nn < t_sources.size(); ++nn)
541  {
542  int t = t_sources[nn];
543  Complex weight; rng_map_obj.get(t,weight);
544  buf[nn].val.data().mat(colorvec_sink,colorvec_source) = conj(weight) * hsum[t];
545  }
546 
547  } // for colorvec_sink
548  } // for colorvec_source
549 
550  QDPIO::cout << "insert: spin_source= " << spin_source << " spin_sink= " << spin_sink << std::endl;
551  for(int nn=0; nn < t_sources.size(); ++nn)
552  {
553  int t = t_sources[nn];
554  qdp_db.insert(buf[nn].key, buf[nn].val);
555  }
556 
557 
558  } // for spin_sink
559  } // for spin_source
560 
561  swatch.stop();
562  QDPIO::cout << "Matrix elements computed: time= "
563  << swatch.getTimeInSeconds()
564  << " secs" << std::endl;
565  } // tt
566 
567  }
568  catch (const std::string& e)
569  {
570  QDPIO::cout << name << ": caught exception around qprop: " << e << std::endl;
571  QDP_abort(1);
572  }
573 
574  push(xml_out,"Relaxation_Iterations");
575  write(xml_out, "ncg_had", ncg_had);
576  pop(xml_out);
577 
578  pop(xml_out); // prop_colorvec
579 
580  // Reset the seed
581  QDP::RNG::setrn(ran_seed);
582 
583  snoop.stop();
584  QDPIO::cout << name << ": total time = "
585  << snoop.getTimeInSeconds()
586  << " secs" << std::endl;
587 
588  QDPIO::cout << name << ": ran successfully" << std::endl;
589 
590  END_CODE();
591  }
592 
593  }
594 
595 } // namespace Chroma
Inline measurement factory.
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.
Serializable value harness.
Definition: key_val_db.h:69
Fourier transform phase factor support.
Definition: sftmom.h:35
const Set & getSet() const
The set to be used in sumMulti.
Definition: sftmom.h:57
static T & Instance()
Definition: singleton.h:432
EXTERN Real dt
Class structure for fermion actions.
Fermion action factories.
All Wilson-type fermion actions.
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.
Complex zN_rng(int N)
Z(N)-rng.
Definition: zN_src.cc:12
Compute the annihilation diagram propagator elements M^-1 * multi1d<LatticeColorVector>
MODS_t & eigen_source
Eigenvectors.
Key for propagator colorstd::vector sources.
Key for propagator colorstd::vector matrix elements.
Key and values for DB.
void setrn(int iseed[4])
Definition: make_seeds.cc:38
void savern(int iseed[4])
Definition: make_seeds.cc:46
Make xml file writer.
int t
Definition: meslate.cc:37
Named object function std::map.
static bool registered
Local registration flag.
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.
multi1d< LatticeColorMatrix > P
QDPSubTypeTrait< typename BinaryReturn< C1, C2, FnLocalInnerProduct >::Type_t >::Type_t localInnerProduct(const QDPSubType< T1, C1 > &l, const QDPType< T2, C2 > &r)
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
multi1d< LatticeFermion > chi(Ncb)
multi1d< SubsetVectorWeight_t > getEigenValues(const MapObject< int, EVPair< LatticeColorVector > > &eigen_source, int num_vecs)
pop(xml_out)
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
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
::std::string string
Definition: gtest.h:1979
Print out basic info about this program.
Fourier transform phase factor support.
GroupXML_t invParam
Definition: qprop_io.h:84
GroupXML_t fermact
Definition: qprop_io.h:80
A Pair type.
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Holds of vectors and weights.
Volume source of Z(N) noise.