CHROMA
inline_stoch_group_meson_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Inline measurement of stochastic group meson operator
3  *
4  */
5 
6 #include "handle.h"
17 #include "meas/glue/mesplq.h"
19 #include "util/ferm/diractodr.h"
20 #include "util/ft/sftmom.h"
21 #include "util/info/proginfo.h"
23 #include <sstream>
24 
26 
27 namespace Chroma
28 {
29  /*!
30  * \ingroup hadron
31  *
32  * @{
33  */
34  namespace InlineStochGroupMesonEnv
35  {
36  //! Number of quarks to be used in this construction
37  const int N_quarks = 2;
38 
39 
40  //
41  // The spin basis matrix to goto Dirac
42  //
43  SpinMatrix rotate_mat(adj(DiracToDRMat()));
44 
45  // Reader for input parameters
46  void read(XMLReader& xml, const std::string& path, InlineStochGroupMesonEnv::Params::Param_t& param)
47  {
48  XMLReader paramtop(xml, path);
49 
50  int version;
51  read(paramtop, "version", version);
52 
53  switch (version)
54  {
55  case 1:
56  /**************************************************************************/
57  break;
58 
59  default :
60  /**************************************************************************/
61 
62  QDPIO::cerr << "Input parameter version " << version << " unsupported." << std::endl;
63  QDP_abort(1);
64  }
65 
66  read(paramtop, "creationOpContractType", param.creat_op_contract_type);
67  read(paramtop, "annihilationOpContractType", param.annih_op_contract_type);
68  read(paramtop, "mom2_max", param.mom2_max);
69  read(paramtop, "displacement_length", param.displacement_length);
70 
71  param.quark_smearing = readXMLGroup(paramtop, "QuarkSmearing", "wvf_kind");
72  param.link_smearing = readXMLGroup(paramtop, "LinkSmearing", "LinkSmearingType");
73  param.quark_dils = readXMLArrayGroup(paramtop, "QuarkDilutions", "DilutionType");
74  }
75 
76 
77  // Writer for input parameters
78  void write(XMLWriter& xml, const std::string& path, const InlineStochGroupMesonEnv::Params::Param_t& param)
79  {
80  push(xml, path);
81 
82  int version = 1;
83 
84  write(xml, "version", version);
85  write(xml, "creationOpContractType", param.creat_op_contract_type);
86  write(xml, "annihilationOpContractType", param.annih_op_contract_type);
87  write(xml, "mom2_max", param.mom2_max);
88  write(xml, "displacement_length", param.displacement_length);
89  xml << param.quark_smearing.xml;
90  xml << param.link_smearing.xml;
91 
92  push(xml, "QuarkDilutions");
93 
94  for(int i=0; i < param.quark_dils.size(); ++i)
95  xml << param.quark_dils[i].xml;
96 
97  pop(xml);
98 
99  pop(xml);
100  }
101 
102 
103  //Reader for 2-quark operator file
105  {
106  XMLReader inputtop(xml, path);
107 
108  read(inputtop, "ops_file", input.ops_file);
109  read(inputtop, "id", input.id);
110  }
111 
112 
113  // Writer for 2-quark operator file
115  {
116  push(xml, path);
117  write(xml, "ops_file", input.ops_file);
118  write(xml, "id", input.id);
119  pop(xml);
120  }
121 
122 
123  //! Read named objects
124  void read(XMLReader& xml, const std::string& path, InlineStochGroupMesonEnv::Params::NamedObject_t& input)
125  {
126  XMLReader inputtop(xml, path);
127 
128  read(inputtop, "gauge_id", input.gauge_id);
129  read(inputtop, "operators_file", input.operators_file);
130  read(inputtop, "Quark_ids", input.quark_ids);
131  }
132 
133  //! Write named objects
134  void write(XMLWriter& xml, const std::string& path, const InlineStochGroupMesonEnv::Params::NamedObject_t& input)
135  {
136  push(xml, path);
137 
138  write(xml, "gauge_id", input.gauge_id);
139  write(xml, "operators_file", input.operators_file);
140  write(xml, "Quark_ids", input.quark_ids);
141 
142  pop(xml);
143  }
144  }
145 
146 
147  namespace InlineStochGroupMesonEnv
148  {
149  // Anonymous namespace for registration
150  namespace
151  {
152  AbsInlineMeasurement* createMeasurement(XMLReader& xml_in,
153  const std::string& path)
154  {
155  return new InlineMeas(Params(xml_in, path));
156  }
157 
158  //! Local registration flag
159  bool registered = false;
160  }
161 
162  const std::string name = "STOCH_GROUP_MESON";
163 
164  //! Register all the factories
165  bool registerAll()
166  {
167  bool success = true;
168  if (! registered)
169  {
172  success &= DilutionSchemeEnv::registerAll();
173  success &= TheInlineMeasurementFactory::Instance().registerObject(name, createMeasurement);
174  registered = true;
175  }
176  return success;
177  }
178 
179 
180  //! Anonymous namespace
181  /*! Diagnostic stuff */
182  namespace
183  {
184  StandardOutputStream& operator<<(StandardOutputStream& os, const multi1d<int>& d)
185  {
186  if (d.size() > 0)
187  {
188  os << d[0];
189 
190  for(int i=1; i < d.size(); ++i)
191  os << " " << d[i];
192  }
193 
194  return os;
195  }
196  }
197 
198 
199  //----------------------------------------------------------------------------
200  // Param stuff
202  {
203  frequency = 0;
204  param.mom2_max = 0;
205  }
206 
207  Params::Params(XMLReader& xml_in, const std::string& path)
208  {
209  try
210  {
211  XMLReader paramtop(xml_in, path);
212 
213  if (paramtop.count("Frequency") == 1)
214  read(paramtop, "Frequency", frequency);
215  else
216  frequency = 1;
217 
218  // Read program parameters
219  read(paramtop, "Param", param);
220 
221  // Read in the output propagator/source configuration info
222  read(paramtop, "NamedObject", named_obj);
223 
224  // Possible alternate XML file pattern
225  if (paramtop.count("xml_file") != 0)
226  {
227  read(paramtop, "xml_file", xml_file);
228  }
229  }
230  catch(const std::string& e)
231  {
232  QDPIO::cerr << __func__ << ": Caught Exception reading XML: " << e << std::endl;
233  QDP_abort(1);
234  }
235  }
236 
237 
238  void
239  Params::writeXML(XMLWriter& xml_out, const std::string& path)
240  {
241  push(xml_out, path);
242 
243  // Parameters for source construction
244  write(xml_out, "Param", param);
245 
246  // Write out the output propagator/source configuration info
247  write(xml_out, "NamedObject", named_obj);
248 
249  pop(xml_out);
250  }
251 
252 
253 
254  //--------------------------------------------------------------
255  //! 2-quark operator structure
257  {
259  {
260  struct QuarkInfo_t
261  {
262  multi1d<int> displacement; /*!< Orig plus/minus 1-based directional displacements */
263  int spin; /*!< 1-based spin index */
264  };
265 
266  multi1d<QuarkInfo_t> quark; /*!< Displacement and spin for each quark */
267  std::string name; /*!< Name of the 2-quark operator */
268  };
269 
270  multi1d<TwoQuarkOp_t> ops; /*!< 2-quark ops within a file */
271  };
272 
273  //! Write quark
274  void write(XMLWriter& xml, const std::string& path,
276  {
277  push(xml, path);
278 
279  write(xml, "Displacement", input.displacement);
280  write(xml, "Spin", input.spin);
281 
282  pop(xml);
283  }
284 
285  //! Write two quark op
286  void write(XMLWriter& xml, const std::string& path,
287  const TwoQuarkOps_t::TwoQuarkOp_t& input)
288  {
289  push(xml, path);
290 
291  write(xml, "Quarks", input.quark);
292  write(xml, "Name", input.name);
293 
294  pop(xml);
295  }
296 
297 
298 
299  //---------------------------------------------------------------------------
300 
301 
302  //! The key for smeared quarks
304  {
305  int t0; /*!< Time of source */
306  int dil; /*!< dilution component per timeslice */
307  };
308 
309 
310  //! Support for the keys of smeared quarks
312  {
313  multi1d<int> lga(2);
314  lga[0] = a.t0;
315  lga[1] = a.dil;
316 
317  multi1d<int> lgb(2);
318  lgb[0] = b.t0;
319  lgb[1] = b.dil;
320 
321  return (lga < lgb);
322  }
323 
324 
326  {
327  LatticeFermion quark;
328  };
329 
330  //----------------------------------------------------------------------------
331  //! The key for smeared and displaced color vectors
333  {
334  int t0; /*!< Time of source */
335  int dil; /*!< dilution component per timeslice */
336 
337  multi1d<int> displacement; /*!< Orig plus/minus 1-based directional displacements */
338  int spin; /*!< 1-based spin index */
339  };
340 
341 
342  //! Support for the keys of smeared and displaced color vectors
344  {
345  multi1d<int> lgaa(3);
346  lgaa[0] = a.spin;
347  lgaa[1] = a.t0;
348  lgaa[2] = a.dil;
349  multi1d<int> lga = concat(lgaa, a.displacement);
350 
351  multi1d<int> lgbb(3);
352  lgbb[0] = b.spin;
353  lgbb[1] = b.t0;
354  lgbb[2] = b.dil;
355  multi1d<int> lgb = concat(lgbb, b.displacement);
356 
357  return (lga < lgb);
358  }
359 
360 
361  //! The value of the std::map
363  {
364  multi1d<LatticeComplex> vec;
365  };
366 
367 
368  //----------------------------------------------------------------------------
369  //! The smeared and displaced objects
371  {
372  public:
373  //! Constructor from smeared std::map
374  SmearedDispObjects(int disp_length,
375  multi1d< Handle< DilutionScheme<LatticeFermion> > > dil_quarks,
377  const multi1d<LatticeColorMatrix> & u_smr);
378 
379  //! Destructor
381 
382  //! Accessor
383  virtual multi1d<LatticeComplex> getDispSource(int quark_num,
384  const KeySmearedDispColorVector_t& key);
385 
386  //! Accessor
387  virtual multi1d<LatticeComplex> getDispSolution(int quark_num,
388  const KeySmearedDispColorVector_t& key);
389 
390  protected:
391  //! Displace an object
392  virtual const multi1d<LatticeComplex>& displaceObject(
393  std::map<KeySmearedDispColorVector_t, SmearedDispColorVector_t>& disp_quark_map,
394  const KeySmearedDispColorVector_t& key,
395  const LatticeFermion& smrd_q);
396 
397  //! Smear sources and solutions
398  virtual const LatticeFermion&
399  smearSource(int qnum , const KeySmearedQuark_t & key);
400 
401  virtual const LatticeFermion&
402  smearSolution(int qnum , const KeySmearedQuark_t & key);
403 
404  private:
405  multi1d< Handle< DilutionScheme<LatticeFermion> > > diluted_quarks;
407 
408  //!Gauge field
409  const multi1d<LatticeColorMatrix> & u;
410 
411  //! Displacement length
413 
414 
415  //! Maps of smeared color vectors
416  multi1d< std::map<KeySmearedQuark_t, SmearedQuark_t> > smeared_src_maps;
417  multi1d< std::map<KeySmearedQuark_t, SmearedQuark_t> > smeared_soln_maps;
418 
419 
420  //!Maps of smeared displaced color vectors
421  multi1d< std::map<KeySmearedDispColorVector_t, SmearedDispColorVector_t> > disp_src_maps;
422  multi1d< std::map<KeySmearedDispColorVector_t, SmearedDispColorVector_t> > disp_soln_maps;
423  };
424 
425 
426  // Constructor from smeared std::map
428  multi1d< Handle< DilutionScheme<LatticeFermion> > > dil_quarks,
430  const multi1d<LatticeColorMatrix> & u_smr) :
431  displacement_length(disp_length),diluted_quarks(dil_quarks),
432  quarkSmearing(qsmr), u(u_smr)
433  {
434  if (diluted_quarks.size() != N_quarks)
435  {
436  QDPIO::cerr << __func__ << ": expected num_quarks=" << N_quarks << std::endl;
437  QDP_abort(1);
438  }
439 
440  smeared_src_maps.resize(N_quarks);
441  smeared_soln_maps.resize(N_quarks);
442 
443  disp_src_maps.resize(N_quarks);
444  disp_soln_maps.resize(N_quarks);
445  }
446 
447 
448  const LatticeFermion&
450  const KeySmearedQuark_t & key)
451  {
452  std::map<KeySmearedQuark_t, SmearedQuark_t> & qmap = smeared_src_maps[qnum];
453 
454  //If entry is not in std::map create it
455  if ( qmap.find(key) == qmap.end() )
456  {
457  // Insert an empty entry and then modify it. This saves on
458  // copying the data around
459  {
460  SmearedQuark_t smrd_empty;
461  qmap.insert(std::make_pair(key, smrd_empty));
462 
463  // Sanity check - the entry better be there
464  if ( qmap.find(key) == qmap.end() )
465  {
466  QDPIO::cerr << __func__
467  << ": internal error - could not insert empty key in std::map"
468  << std::endl;
469  QDP_abort(1);
470  }
471  }
472 
473  // Modify the previous empty entry
474  SmearedQuark_t& smrd_q = qmap.find(key)->second;
475 
476  StopWatch snoop;
477 
478  snoop.reset();
479  snoop.start();
480 
481  LatticeFermion f = diluted_quarks[qnum]->dilutedSource(key.t0, key.dil);
482 
483  (*quarkSmearing)(f, u);
484 
485  // The first set of quarks will get a gamma_5 to implement the oper we need
486  // op= eta_1^dag * [whatever gamma and U oper structure] * gamma_5 eta_0^dag
487  SpinMatrix mat;
488  if (qnum == 0)
489  {
490  mat = rotate_mat * Gamma(15);
491  }
492  else
493  {
494  mat = rotate_mat;
495  }
496 
497  smrd_q.quark = mat * f;
498 
499  snoop.stop();
500 
501  QDPIO::cout << " Smeared Sources: Quark = "<< qnum <<" t0 = "
502  << key.t0 <<" dil = "<< key.dil << " Time = "<< snoop.getTimeInSeconds() <<" sec"<<std::endl;
503 
504  // Insert
505  qmap.insert(std::make_pair(key, smrd_q));
506 
507  } // if find in std::map
508 
509  // The key now must exist in the std::map, so return the smeared quark field
510 
511  return (qmap.find(key)->second).quark ;
512  }
513 
514 
515 
516  const LatticeFermion&
518  const KeySmearedQuark_t & key)
519  {
520  std::map<KeySmearedQuark_t, SmearedQuark_t> & qmap = smeared_soln_maps[qnum];
521 
522  //If entry is not in std::map create it
523  if ( qmap.find(key) == qmap.end() )
524  {
525 
526  // Insert an empty entry and then modify it. This saves on
527  // copying the data around
528  {
529  SmearedQuark_t smrd_empty;
530  qmap.insert(std::make_pair(key, smrd_empty));
531 
532  // Sanity check - the entry better be there
533  if ( qmap.find(key) == qmap.end() )
534  {
535  QDPIO::cerr << __func__
536  << ": internal error - could not insert empty key in std::map"
537  << std::endl;
538  QDP_abort(1);
539  }
540  }
541 
542  // Modify the previous empty entry
543  SmearedQuark_t& smrd_q = qmap.find(key)->second;
544 
545  StopWatch snoop;
546  snoop.reset();
547  snoop.start();
548 
549  LatticeFermion f = diluted_quarks[qnum]->dilutedSolution(key.t0, key.dil);
550 
551  (*quarkSmearing)(smrd_q.quark, u);
552 
553  // The first set of quarks will get a gamma_5 to implement the oper. we need
554  // op= psi_0^dag * gamma_5 * [whatever gamma and U oper structure] * psi_1^dag
555  SpinMatrix mat;
556  if (qnum == 0)
557  {
558  mat = rotate_mat * Gamma(15);
559  }
560  else
561  {
562  mat = rotate_mat;
563  }
564 
565  smrd_q.quark = mat * f;
566 
567  snoop.stop();
568 
569  QDPIO::cout << " Smeared Sinks: Quark = "<< qnum <<" t0 = "
570  << key.t0 <<" dil = "<< key.dil << " Time = "<< snoop.getTimeInSeconds() <<" sec"<<std::endl;
571 
572  // Insert
573  qmap.insert(std::make_pair(key, smrd_q));
574 
575  } // if find in std::map
576 
577  // The key now must exist in the std::map, so return the smeared quark field
578 
579  return (qmap.find(key)->second).quark ;
580  }
581 
582  //! Accessor
583  multi1d<LatticeComplex>
585  const KeySmearedDispColorVector_t& key)
586  {
587  //Get Smeared quark
588  KeySmearedQuark_t smr_key;
589  smr_key.t0 = key.t0;
590  smr_key.dil = key.dil;
591 
592  const LatticeFermion& smrd_q = smearSource(quark_num, smr_key);
593 
594  multi1d<LatticeComplex> vec;
595 
596  //Check if any displacement is needed
597  if (displacement_length == 0)
598  {
599  vec.resize(Nc);
600  LatticeColorVector colvec = peekSpin( smrd_q, key.spin - 1);
601 
602  for (int c = 0 ; c < Nc ; ++c)
603  vec[c] = peekColor( colvec, c);
604 
605  }
606  else
607  {
608  vec = displaceObject( disp_src_maps[quark_num] , key , smrd_q);
609  }
610 
611  return vec;
612  }
613 
614  //! Accessor
615  multi1d<LatticeComplex>
617  const KeySmearedDispColorVector_t& key)
618  {
619  //Get Smeared quark
620  KeySmearedQuark_t smr_key;
621  smr_key.t0 = key.t0;
622  smr_key.dil = key.dil;
623 
624  const LatticeFermion& smrd_q = smearSolution(quark_num, smr_key);
625 
626  multi1d<LatticeComplex> vec;
627 
628  //Check if any displacement is needed
629  if (displacement_length == 0)
630  {
631  vec.resize(Nc);
632  LatticeColorVector colvec = peekSpin( smrd_q, key.spin - 1);
633 
634  for (int c = 0 ; c < Nc ; ++c)
635  vec[c] = peekColor( colvec, c);
636 
637  }
638  else
639  {
640  vec = displaceObject( disp_soln_maps[quark_num] , key , smrd_q);
641  }
642 
643  return vec;
644  }
645 
646 
647  //! Accessor
648  const multi1d<LatticeComplex> &
650  std::map<KeySmearedDispColorVector_t, SmearedDispColorVector_t>& disp_quark_map,
651  const KeySmearedDispColorVector_t& key,
652  const LatticeFermion& smrd_q)
653  {
654  StopWatch snoop;
655 
656  // If no entry, then create a displaced version of the quark
657  if (disp_quark_map.find(key) == disp_quark_map.end())
658  {
659  // std::cout << __func__
660  // << ": n=" << n
661  // << " l=" << l
662  // << " i=" << i
663  // << " disp=" << term.quark[i].displacement
664  // << " len=" << term.quark[i].disp_len
665  // << " dir=" << term.quark[i].disp_dir
666  // << std::endl;
667 
668  // Insert an empty entry and then modify it. This saves on
669  // copying the data around
670  {
671  SmearedDispColorVector_t disp_empty;
672 
673  snoop.reset();
674  snoop.start();
675 
676  disp_quark_map.insert(std::make_pair(key, disp_empty));
677 
678  snoop.stop();
679 
680  QDPIO::cout<<"Inserted key in std::map: time = "<< snoop.getTimeInSeconds() << "secs"<<std::endl;
681 
682  // Sanity check - the entry better be there
683  if (disp_quark_map.find(key) == disp_quark_map.end())
684  {
685  QDPIO::cerr << __func__
686  << ": internal error - could not insert empty key in std::map"
687  << std::endl;
688  QDP_abort(1);
689  }
690  }
691 
692  // Modify the previous empty entry
693  SmearedDispColorVector_t& disp_q = disp_quark_map.find(key)->second;
694 
695  snoop.reset();
696  snoop.start();
697 
698  //Chroma uses a zero-based spin convention
699  LatticeColorVector vec = peekSpin(smrd_q, key.spin - 1);
700 
701  for(int i=0; i < key.displacement.size(); ++i)
702  {
703  if (key.displacement[i] > 0)
704  {
705  int disp_dir = key.displacement[i] - 1;
706  int disp_len = displacement_length;
707  displacement(u, vec, disp_len, disp_dir);
708  }
709  else if (key.displacement[i] < 0)
710  {
711  int disp_dir = -key.displacement[i] - 1;
712  int disp_len = -displacement_length;
713  displacement(u, vec, disp_len, disp_dir);
714  }
715  }
716 
717  snoop.stop();
718 
719  QDPIO::cout << "Displaced Quarks: Spin = "<<key.spin<<" Disp = "
720  << key.displacement <<" Time = "<<snoop.getTimeInSeconds() <<" sec"<<std::endl;
721 
722  disp_q.vec.resize(Nc);
723 
724  for(int i = 0 ; i < Nc ; ++i )
725  {
726  disp_q.vec[i] = peekColor(vec, i);
727  }
728 
729  } // if find in std::map
730 
731  snoop.reset();
732  snoop.start();
733 
734  // The key now must exist in the std::map, so return the std::vector
735  SmearedDispColorVector_t& disp_q = disp_quark_map.find(key)->second;
736 
737  snoop.stop();
738 
739  QDPIO::cout << "Retrieved entry from std::map : time = "<< snoop.getTimeInSeconds() << "secs "<<std::endl;
740 
741  return disp_q.vec;
742  }
743 
744 
745  //----------------------------------------------------------------------------
746  // Color contract two fermions
747  void makeColorSinglet (LatticeComplex& singlet,
748  const multi1d<LatticeComplex>& q0,
749  const multi1d<LatticeComplex>& q1,
750  const Subset& subset)
751  {
752  singlet[subset] = q0[0] * q1[0];
753  for(int i=1; i < q0.size(); ++i)
754  singlet[subset] += q0[i] * q1[i];
755  }
756 
757 
758  //----------------------------------------------------------------------------
759  // Color contract two fermions
760  multi2d<DComplex> contractOp (SmearedDispObjects& smrd_disp_vecs,
761  int n0, const KeySmearedDispColorVector_t& k0,
762  int n1, const KeySmearedDispColorVector_t& k1,
763  MesonOpType contractType,
764  const SftMom& phases,
765  int t0)
766  {
767  multi2d<DComplex> op_sum;
768  LatticeComplex singlet;
769 
770  switch(contractType)
771  {
773  makeColorSinglet(singlet,
774  smrd_disp_vecs.getDispSource(n0, k0),
775  smrd_disp_vecs.getDispSource(n1, k1),
776  phases.getSet()[t0]);
777  op_sum = phases.sft(singlet, t0);
778  break;
779 
781  makeColorSinglet(singlet,
782  smrd_disp_vecs.getDispSource(n0, k0),
783  smrd_disp_vecs.getDispSolution(n1, k1),
784  phases.getSet()[t0]);
785  op_sum = phases.sft(singlet, t0);
786  break;
787 
789  makeColorSinglet(singlet,
790  smrd_disp_vecs.getDispSolution(n0, k0),
791  smrd_disp_vecs.getDispSource(n1, k1),
792  phases.getSet()[t0]);
793  op_sum = phases.sft(singlet, t0);
794  break;
795 
797  makeColorSinglet(singlet,
798  smrd_disp_vecs.getDispSolution(n0, k0),
799  smrd_disp_vecs.getDispSolution(n1, k1),
800  all);
801  op_sum = phases.sft(singlet);
802  break;
803  }
804 
805  return op_sum;
806  }
807 
808 
809  //----------------------------------------------------------------------------
810  //! Meson operator
812  {
813  //! Meson operator time slices corresponding to location of operator source
815  {
816  //! Meson operator dilutions
817  struct Dilutions_t
818  {
819  //! Momentum projected correlator
820  struct Mom_t
821  {
822  multi1d<int> mom; /*!< D-1 momentum of this operator */
823  multi1d<DComplex> op; /*!< Momentum projected operator */
824  };
825 
826  multi1d<Mom_t> mom_projs; /*!< Holds momentum projections of the operator */
827  };
828 
829  multi2d<Dilutions_t> dilutions; /*!< Hybrid list indices */
830  int t0; /*!< Actual time corresponding to this timeslice */
831  };
832 
833  GroupXML_t smearing; /*!< String holding quark smearing xml */
834 
835  Seed seed_l; /*!< Id of left quark */
836  Seed seed_r; /*!< Id of right quark */
837 
838  GroupXML_t dilution_l; /*!< Dilution scheme of left quark */
839  GroupXML_t dilution_r; /*!< Dilution scheme of right quark */
840 
841  std::string id; /*!< Tag/ID used in analysis codes */
842 
843  MesonOpType op_contract_type; /*<! Contraction type for creation op */
844  int mom2_max; /*!< |\vec{p}|^2 */
845  int decay_dir; /*!< Direction of decay */
846 
847  multi1d<TimeSlices_t> time_slices; /*!< Time slices of the lattice that are used */
848  };
849 
850 
851  //----------------------------------------------------------------------------
852  //! MesonOperator header writer
853  void write(XMLWriter& xml, const std::string& path, const MesonOperator_t& param)
854  {
855  push(xml, path);
856 
857  int version = 1;
858  write(xml, "version", version);
859  write(xml, "id", param.id);
860  write(xml, "opContractType", param.op_contract_type);
861  write(xml, "mom2_max", param.mom2_max);
862  write(xml, "decay_dir", param.decay_dir);
863  write(xml, "seed_l", param.seed_l);
864  write(xml, "seed_r", param.seed_r);
865 
866  push(xml, "dilution_l");
867  xml << param.dilution_l.xml;
868  pop(xml);
869 
870  push(xml, "dilution_r");
871  xml << param.dilution_r.xml;
872  pop(xml);
873 
874  xml << param.smearing.xml;
875 
876  pop(xml);
877  }
878 
879 
880 
881  //! MesonOperator binary writer
882  void write(BinaryWriter& bin, const MesonOperator_t::TimeSlices_t::Dilutions_t::Mom_t& param)
883  {
884  write(bin, param.mom);
885  write(bin, param.op);
886  }
887 
888  //! MesonOperator binary writer
889  void write(BinaryWriter& bin, const MesonOperator_t::TimeSlices_t::Dilutions_t& param)
890  {
891  write(bin, param.mom_projs);
892  }
893 
894  //! MesonOperator binary writer
895  void write(BinaryWriter& bin, const MesonOperator_t::TimeSlices_t& param)
896  {
897  write(bin, param.dilutions);
898  write(bin, param.t0);
899  }
900 
901  //! MesonOperator binary writer
902  void write(BinaryWriter& bin, const MesonOperator_t& param)
903  {
904  write(bin, param.seed_l);
905  write(bin, param.seed_r);
906  write(bin, param.mom2_max);
907  write(bin, param.decay_dir);
908  write(bin, param.time_slices);
909  }
910 
911 
912  //! Read 2-quark operators file, assign correct displacement length
913  void readOps(TwoQuarkOps_t& oplist,
914  const std::string& ops_file)
915  {
916  START_CODE();
917 
918  TextFileReader reader(ops_file);
919 
920  int num_ops;
921  reader >> num_ops;
922  oplist.ops.resize(num_ops);
923 
924  //Loop over ops within a file
925  for(int n=0; n < oplist.ops.size(); ++n)
926  {
927  std::string name;
928  reader >> name;
929  oplist.ops[n].name = name;
930 
931  TwoQuarkOps_t::TwoQuarkOp_t& qqq = oplist.ops[n];
932  qqq.quark.resize(N_quarks);
933 
934  // Read 1-based spin
935  reader >> qqq.quark[0].spin >> qqq.quark[1].spin;
936 
937  // Read 1-based displacement only for the right quark
938  qqq.quark[0].displacement.resize(1);
939  qqq.quark[0].displacement[0] = 0;
940 
941  int ndisp;
942  reader >> ndisp;
943  multi1d<int> displacement(ndisp);
944  for(int i=0; i < ndisp; ++i)
945  reader >> displacement[i];
946 
947  qqq.quark[1].displacement = displacement;
948  } //n
949 
950  reader.close();
951 
952  END_CODE();
953  } //void readOps
954 
955 
956  //-------------------------------------------------------------------------------
957  // Function call
958  void
959  InlineMeas::operator()(unsigned long update_no,
960  XMLWriter& xml_out)
961  {
962  // If xml file not empty, then use alternate
963  if (params.xml_file != "")
964  {
965  std::string xml_file = makeXMLFileName(params.xml_file, update_no);
966 
967  push(xml_out, "stoch_meson");
968  write(xml_out, "update_no", update_no);
969  write(xml_out, "xml_file", xml_file);
970  pop(xml_out);
971 
972  XMLFileWriter xml(xml_file);
973  func(update_no, xml);
974  }
975  else
976  {
977  func(update_no, xml_out);
978  }
979  }
980 
981 
982  // Function call
983  void
984  InlineMeas::func(unsigned long update_no,
985  XMLWriter& xml_out)
986  {
987  START_CODE();
988 
989  StopWatch snoop;
990  snoop.reset();
991  snoop.start();
992 
993 
994  StopWatch swiss;
995 
996  // Test and grab a reference to the gauge field
997  XMLBufferWriter gauge_xml;
998  try
999  {
1000  TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
1001  TheNamedObjMap::Instance().get(params.named_obj.gauge_id).getRecordXML(gauge_xml);
1002  }
1003  catch( std::bad_cast )
1004  {
1005  QDPIO::cerr << InlineStochGroupMesonEnv::name << ": caught dynamic cast error"
1006  << std::endl;
1007  QDP_abort(1);
1008  }
1009  catch (const std::string& e)
1010  {
1011  QDPIO::cerr << InlineStochGroupMesonEnv::name << ": std::map call failed: " << e
1012  << std::endl;
1013  QDP_abort(1);
1014  }
1015  const multi1d<LatticeColorMatrix>& u =
1016  TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
1017 
1018  push(xml_out, "StochGroupMeson");
1019  write(xml_out, "update_no", update_no);
1020 
1021  QDPIO::cout << InlineStochGroupMesonEnv::name << ": Stochastic Meson Operator" << std::endl;
1022 
1023  proginfo(xml_out); // Print out basic program info
1024 
1025  // Write out the input
1026  params.writeXML(xml_out, "Input");
1027 
1028  // Write out the config info
1029  write(xml_out, "Config_info", gauge_xml);
1030 
1031  push(xml_out, "Output_version");
1032  write(xml_out, "out_version", 1);
1033  pop(xml_out);
1034 
1035  //First calculate some gauge invariant observables just for info.
1036  //This is really cheap.
1037  MesPlq(xml_out, "Observables", u);
1038 
1039  // Save current seed
1040  Seed ran_seed;
1041  QDP::RNG::savern(ran_seed);
1042 
1043  //
1044  // Construct the dilution scheme for each of the 2 quarks
1045  //
1046  if (params.param.quark_dils.size() != N_quarks)
1047  {
1048  QDPIO::cerr << name << ": expecting 2 quark dilutions" << std::endl;
1049  QDP_abort(1);
1050  }
1051 
1052  multi1d< Handle< DilutionScheme<LatticeFermion> > > diluted_quarks(N_quarks); /*!< Here is the big (dil) pickle */
1053 
1054  try
1055  {
1056  // Loop over the 2 quark dilutions
1057  for(int n = 0; n < params.param.quark_dils.size(); ++n)
1058  {
1059  const GroupXML_t& dil_xml = params.param.quark_dils[n];
1060 
1061  std::istringstream xml_d(dil_xml.xml);
1062  XMLReader diltop(xml_d);
1063  QDPIO::cout << "Dilution type = " << dil_xml.id << std::endl;
1064 
1065  diluted_quarks[n] = TheFermDilutionSchemeFactory::Instance().createObject(
1066  dil_xml.id, diltop, dil_xml.path);
1067  }
1068  }
1069  catch(const std::string& e)
1070  {
1071  QDPIO::cerr << name << ": Caught Exception constructing dilution scheme: " << e << std::endl;
1072  QDP_abort(1);
1073  }
1074 //------------------------------------------------------------------------
1075 //Sanity checks
1076 
1077  //The participating timeslices must match for each quark
1078  //grab info from first quark to prime the work
1079 
1080  multi1d<int> participating_timeslices( diluted_quarks[0]->getNumTimeSlices() );
1081 
1082  for (int t0 = 0 ; t0 < participating_timeslices.size() ; ++t0)
1083  {
1084  participating_timeslices[t0] = diluted_quarks[0]->getT0(t0);
1085  }
1086 
1087  for (int n = 1 ; n < N_quarks ; ++n)
1088  {
1089  if ( diluted_quarks[n]->getNumTimeSlices() != participating_timeslices.size() )
1090  {
1091  QDPIO::cerr << name << " : Quarks do not contain the same number of dilution timeslices: Quark "
1092  << n << std::endl;
1093 
1094  QDP_abort(1);
1095  }
1096 
1097  for (int t0 = 0 ; t0 < participating_timeslices.size() ; ++t0)
1098  {
1099  if ( diluted_quarks[n]->getT0(t0) != participating_timeslices[t0] )
1100  {
1101  QDPIO::cerr << name << " : Quarks do not contain the same participating timeslices: Quark "<<
1102  n << " timeslice "<< t0 << std::endl;
1103 
1104  QDP_abort(1);
1105  }
1106  }
1107  }
1108 
1109  //Another Sanity check, the two quarks must all be
1110  //inverted on the same cfg
1111  for (int n = 1 ; n < N_quarks ; ++n)
1112  {
1113  if (diluted_quarks[0]->getCfgInfo() != diluted_quarks[n]->getCfgInfo())
1114  {
1115  QDPIO::cerr << name
1116  << " : Quarks do not contain the same cfg info, quark "<< n << std::endl;
1117 
1118  QDP_abort(1);
1119  }
1120 
1121  }
1122 
1123  //Also ensure that the cfg on which the inversions were performed
1124  //is the same as the cfg that we are using
1125  {
1126  std::string cfgInfo;
1127 
1128  //Really ugly way of doing this(Robert Heeeelp!!)
1129  XMLBufferWriter top;
1130  write(top, "Config_info", gauge_xml);
1131  XMLReader from(top);
1132  XMLReader from2(from, "/Config_info");
1133  std::ostringstream os;
1134  from2.print(os);
1135 
1136  cfgInfo = os.str();
1137 
1138  if (cfgInfo != diluted_quarks[0]->getCfgInfo())
1139  {
1140  QDPIO::cerr << name
1141  << " : Quarks do not contain the same cfg info as the gauge field."
1142  << "gauge: XX"<<cfgInfo<<"XX quarks: XX"<<diluted_quarks[0]->getCfgInfo()<<"XX"<< std::endl;
1143 
1144 
1145  QDP_abort(1);
1146  }
1147  }
1148 
1149  //
1150  // Initialize the slow Fourier transform phases
1151  //
1152  int decay_dir = diluted_quarks[0]->getDecayDir();
1153 
1154  SftMom phases(params.param.mom2_max, false, decay_dir);
1155 
1156  // Sanity check - if this doesn't work we have serious problems
1157  if (phases.numSubsets() != QDP::Layout::lattSize()[decay_dir])
1158  {
1159  QDPIO::cerr << name << ": number of time slices not equal to that in the decay direction: "
1160  << QDP::Layout::lattSize()[decay_dir]
1161  << std::endl;
1162  QDP_abort(1);
1163  }
1164 
1165 
1166  // Another sanity check. The seeds of all the quarks must be different
1167  // and thier decay directions must be the same
1168  for(int n = 1 ; n < diluted_quarks.size(); ++n)
1169  {
1170  if ( toBool( diluted_quarks[n]->getSeed() == diluted_quarks[0]->getSeed() ) )
1171  {
1172  QDPIO::cerr << name << ": error, quark seeds are the same" << std::endl;
1173  QDP_abort(1);
1174  }
1175 
1176  if ( toBool( diluted_quarks[n]->getDecayDir() != diluted_quarks[0]->getDecayDir() ) )
1177  {
1178  QDPIO::cerr << name << ": error, quark decay dirs do not match" <<std::endl;
1179  QDP_abort(1);
1180  }
1181 
1182  }
1183 
1184  //
1185  // Smear the gauge field if needed
1186  //
1187  multi1d<LatticeColorMatrix> u_smr = u;
1188 
1189  try
1190  {
1191  std::istringstream xml_l(params.param.link_smearing.xml);
1192  XMLReader linktop(xml_l);
1193  QDPIO::cout << "Link smearing type = " << params.param.link_smearing.id << std::endl;
1194 
1195 
1197  linkSmearing(TheLinkSmearingFactory::Instance().createObject(params.param.link_smearing.id,
1198  linktop, params.param.link_smearing.path));
1199 
1200  (*linkSmearing)(u_smr);
1201  }
1202  catch(const std::string& e)
1203  {
1204  QDPIO::cerr << name << ": Caught Exception link smearing: " << e << std::endl;
1205  QDP_abort(1);
1206  }
1207 
1208  MesPlq(xml_out, "Smeared_Observables", u_smr);
1209 
1210  //
1211  // Read operator coefficients
1212  //
1213  QDPIO::cout << "Reading 2-quark operators" << std::endl;
1214  TwoQuarkOps_t qqq_oplist;
1215 
1217 
1218  //
1219  // Create the quark smearing factory
1220  //
1221  Handle< QuarkSmearing<LatticeFermion> > quarkSmearing;
1222 
1223  try
1224  {
1225  QDPIO::cout << "Create quark smearing object" << std::endl;
1226 
1227  // Create the quark smearing object
1228  std::istringstream xml_s(params.param.quark_smearing.xml);
1229  XMLReader smeartop(xml_s);
1230 
1231  quarkSmearing =
1233  smeartop, params.param.quark_smearing.path);
1234 
1235  }
1236  catch(const std::string& e)
1237  {
1238  QDPIO::cerr << ": Caught Exception creating quark smearing object: " << e << std::endl;
1239  QDP_abort(1);
1240  }
1241  catch(...)
1242  {
1243  QDPIO::cerr << ": Caught generic exception creating smearing object" << std::endl;
1244  QDP_abort(1);
1245  }
1246 
1247  //
1248  // Meson operators
1249  //
1250 
1251  //We make all source operators before we make all sink operators to
1252  //save on memory.
1253 
1254  for(int t0 = 0; t0 < participating_timeslices.size() ; ++t0)
1255  {
1256  StopWatch watch;
1257 
1258  // Bummer. We have to hold both the smeared sources and the smeared
1259  // solutions simultaneously since the operator construction might
1260  // involve both the sources and the solutions in one operator
1261  // for either the creation and/or annihilation operator
1262 
1263  // The object holding the smeared and displaced color std::vector std::maps
1265  diluted_quarks, quarkSmearing, u_smr);
1266 
1267  //Make the source operators
1268  {
1269  // Creation operator
1270  MesonOperator_t creat_oper;
1272  creat_oper.mom2_max = params.param.mom2_max;
1273  creat_oper.decay_dir = decay_dir;
1274  creat_oper.seed_l = diluted_quarks[1]->getSeed();
1275  creat_oper.seed_r = diluted_quarks[0]->getSeed();
1276  creat_oper.dilution_l = params.param.quark_dils[1];
1277  creat_oper.dilution_r = params.param.quark_dils[0];
1278  creat_oper.smearing = params.param.quark_smearing;
1279  creat_oper.time_slices.resize( 1 ); //Only a single t0 per file
1280 
1281  // Construct creation operator
1282  // Loop over each operator
1283  for(int l=0; l < qqq_oplist.ops.size(); ++l)
1284  {
1285  QDPIO::cout << "Elemental operator: op = " << l << std::endl;
1286 
1287  push(xml_out, "MesonOperator");
1288 
1289  creat_oper.id = qqq_oplist.ops[l].name;
1290 
1291  write(xml_out, "Name", creat_oper.id);
1292 
1293  // Build the operator
1294  swiss.reset();
1295  swiss.start();
1296 
1297  // The keys for the spin and displacements for this particular elemental operator
1298  multi1d<KeySmearedDispColorVector_t> keySmearedDispColorVector(N_quarks);
1299 
1300  for(int n = 0 ; n < N_quarks ; ++n)
1301  {
1302  keySmearedDispColorVector[n].displacement = qqq_oplist.ops[l].quark[n].displacement;
1303  keySmearedDispColorVector[n].spin = qqq_oplist.ops[l].quark[n].spin;
1304  }
1305 
1306 
1307  creat_oper.time_slices[0].t0 = participating_timeslices[t0];
1308 
1309  const int n0 = 1;
1310  const int n1 = 0;
1311 
1312  // The operator must hold all the dilutions
1313 
1314  // Creation operator
1315  MesonOperator_t::TimeSlices_t& cop = creat_oper.time_slices[0];
1316 
1317  cop.dilutions.resize(diluted_quarks[n0]->getDilSize(t0), diluted_quarks[n1]->getDilSize(t0));
1318 
1319  for (int n = 0 ; n < N_quarks ; ++n)
1320  {
1321  keySmearedDispColorVector[n].t0 = t0;
1322  }
1323 
1324  for(int i = 0 ; i < diluted_quarks[n0]->getDilSize(t0) ; ++i)
1325  {
1326  for(int j = 0 ; j < diluted_quarks[n1]->getDilSize(t0) ; ++j)
1327  {
1328  keySmearedDispColorVector[0].dil = i;
1329  keySmearedDispColorVector[1].dil = j;
1330 
1331  LatticeComplex c_oper;
1332 
1333  watch.reset();
1334  watch.start();
1335 
1336  // Do the relevant quark contraction
1337  // Slow fourier-transform
1338  // We can restrict what the FT routine requires to a subset.
1339  multi2d<DComplex> c_sum(contractOp(smrd_disp_vecs,
1340  n0, keySmearedDispColorVector[0],
1341  n1, keySmearedDispColorVector[1],
1343  phases,
1344  participating_timeslices[t0]));
1345 
1346  watch.stop();
1347 
1348  /*QDPIO::cout << " Spatial sums completed: time = " <<
1349  watch.getTimeInSeconds() << "secs " << std::endl;
1350  */
1351  // Unpack into separate momentum and correlator
1352  cop.dilutions(i,j).mom_projs.resize(phases.numMom());
1353 
1354  for(int mom_num = 0 ; mom_num < phases.numMom() ; ++mom_num)
1355  {
1356  cop.dilutions(i,j).mom_projs[mom_num].mom = phases.numToMom(mom_num);
1357  cop.dilutions(i,j).mom_projs[mom_num].op.resize(1);
1358 
1359  cop.dilutions(i,j).mom_projs[mom_num].op[ 0 ] =
1360  c_sum[mom_num][ participating_timeslices[t0] ];
1361  }
1362  } // end for j
1363  } // end for i
1364 
1365  swiss.stop();
1366 
1367 
1368  QDPIO::cout << "Creation operator construction: operator= " << l
1369  << " time= "
1370  << swiss.getTimeInSeconds()
1371  << " secs" << std::endl;
1372 
1373  QDPIO::cout << "Creation operator testval(t0 = " <<
1374  participating_timeslices[t0] << ") = " <<
1375  creat_oper.time_slices[0].dilutions(0,0).mom_projs[0].op[0]
1376  << std::endl;
1377 
1378  //Hard code the elemental op name for now
1379  std::stringstream cnvrt;
1380  cnvrt << creat_oper.id << "_t" << participating_timeslices[t0] << "_src.lime";
1381 
1382  std::string filename;
1383 
1384  filename = cnvrt.str();
1385 
1386  // Write the meta-data and the binary for this operator
1387  swiss.reset();
1388  swiss.start();
1389  {
1390  XMLBufferWriter src_record_xml, file_xml;
1391  BinaryBufferWriter src_record_bin;
1392 
1393  push(file_xml, "SourceMesonOperator");
1394  write(file_xml, "Params", params.param);
1395  write(file_xml, "Config_info", gauge_xml);
1396  write(file_xml, "Op_Info",qqq_oplist.ops[l]);
1397 
1398  push(file_xml, "QuarkSources");
1399 
1400  const int n0 = 1;
1401  const int n1 = 0;
1402 
1403  push(file_xml, "Quark_l");
1404  push(file_xml, "TimeSlice");
1405  push(file_xml, "Dilutions");
1406  for (int dil = 0; dil < diluted_quarks[n0]->getDilSize(t0) ; ++dil)
1407  {
1408  write(file_xml, "elem", diluted_quarks[n0]->getSourceHeader(t0, dil));
1409 
1410  // QDPIO::cout<< "t0 = " << t0 << " dil = "<< dil <<
1411  // " srdhdr = XX"<<diluted_quarks[n0]->getSourceHeader(t0,dil) << std::endl;
1412  }
1413  pop(file_xml); //dilutions
1414  pop(file_xml); //TimeSlice
1415  pop(file_xml); //Quark_l
1416 
1417  push(file_xml, "Quark_r");
1418  push(file_xml, "TimeSlice");
1419  push(file_xml, "Dilutions");
1420  for (int dil = 0; dil < diluted_quarks[n1]->getDilSize(t0) ; ++dil)
1421  {
1422  write(file_xml, "elem", diluted_quarks[n1]->getSourceHeader(t0, dil));
1423  }
1424  pop(file_xml); //dilutions
1425  pop(file_xml); //TimeSlice
1426  pop(file_xml); //Quark_r
1427 
1428  pop(file_xml);//QuarkSources
1429  push(file_xml, "QuarkSinks");
1430 
1431  push(file_xml, "Quark_l");
1432  write(file_xml, "PropHeader", diluted_quarks[n0]->getPropHeader(0,0));
1433  pop(file_xml);
1434 
1435  push(file_xml, "Quark_r");
1436  write(file_xml, "PropHeader", diluted_quarks[n1]->getPropHeader(0,0));
1437  pop(file_xml);
1438 
1439  pop(file_xml); //Quark Sinks
1440  pop(file_xml);//SourceMesonOp
1441 
1442  QDPFileWriter qdp_file(file_xml, filename, // are there one or two files???
1443  QDPIO_SINGLEFILE, QDPIO_SERIAL, QDPIO_OPEN);
1444 
1445 
1446  write(src_record_xml, "MesonCreationOperator", creat_oper);
1447  write(src_record_bin, creat_oper);
1448 
1449  write(qdp_file, src_record_xml, src_record_bin);
1450 
1451  }
1452  swiss.stop();
1453 
1454  QDPIO::cout << "Source Operator writing: operator = " <<
1455  l << " time= "
1456  << swiss.getTimeInSeconds()
1457  << " secs" << std::endl;
1458 
1459  pop(xml_out); // MesonOperator
1460 
1461  } // end for l (operator )
1462 
1463  } //End Make creation operator
1464 
1465 
1466 
1467  //Make Annilation Operator
1468  {
1469  // Annihilation operator
1470  MesonOperator_t annih_oper;
1472  annih_oper.mom2_max = params.param.mom2_max;
1473  annih_oper.decay_dir = decay_dir;
1474  annih_oper.seed_l = diluted_quarks[0]->getSeed();
1475  annih_oper.seed_r = diluted_quarks[1]->getSeed();
1476  annih_oper.dilution_l = params.param.quark_dils[0];
1477  annih_oper.dilution_r = params.param.quark_dils[1];
1478  annih_oper.smearing = params.param.quark_smearing;
1479  annih_oper.time_slices.resize( 1 );
1480 
1481  // Construct annihilation operator
1482  QDPIO::cout << "Building Sink operators" << std::endl;
1483 
1484  // Loop over each operator
1485  for(int l=0; l < qqq_oplist.ops.size(); ++l)
1486  {
1487  QDPIO::cout << "Elemental operator: op = " << l << std::endl;
1488 
1489  annih_oper.id = qqq_oplist.ops[l].name;
1490 
1491  // Build the operator
1492  swiss.reset();
1493  swiss.start();
1494 
1495  // The keys for the spin and displacements for this particular elemental operator
1496  multi1d<KeySmearedDispColorVector_t> keySmearedDispColorVector(N_quarks);
1497 
1498  for(int n = 0 ; n < N_quarks ; ++n)
1499  {
1500  keySmearedDispColorVector[n].displacement = qqq_oplist.ops[l].quark[n].displacement;
1501  keySmearedDispColorVector[n].spin = qqq_oplist.ops[l].quark[n].spin;
1502  }
1503 
1504  annih_oper.time_slices[0].t0 = participating_timeslices[t0];
1505 
1506  const int n0 = 0;
1507  const int n1 = 1;
1508 
1509  // The operator must hold all the dilutions
1510  // We know that all time slices match. However, not all time slices of the
1511  // lattice maybe used
1512 
1513  // Annihilation operator
1514  MesonOperator_t::TimeSlices_t& aop = annih_oper.time_slices[0];
1515 
1516  aop.dilutions.resize(diluted_quarks[n0]->getDilSize(t0), diluted_quarks[n1]->getDilSize(t0));
1517 
1518  for (int n = 0 ; n < N_quarks ; ++n)
1519  {
1520  keySmearedDispColorVector[n].t0 = t0;
1521  }
1522 
1523  for(int i = 0 ; i < diluted_quarks[n0]->getDilSize(t0) ; ++i)
1524  {
1525  for(int j = 0 ; j < diluted_quarks[n1]->getDilSize(t0) ; ++j)
1526  {
1527  keySmearedDispColorVector[0].dil = i;
1528  keySmearedDispColorVector[1].dil = j;
1529 
1530  watch.reset();
1531  watch.start();
1532 
1533  // Contract over color indices
1534  // Do the relevant quark contraction
1535  // Slow fourier-transform
1536  // We can restrict what the FT routine requires to a subset.
1537  multi2d<DComplex> a_sum(contractOp(smrd_disp_vecs,
1538  n0, keySmearedDispColorVector[0],
1539  n1, keySmearedDispColorVector[1],
1541  phases,
1542  participating_timeslices[t0]));
1543 
1544  watch.stop();
1545  /*
1546  QDPIO::cout << "Spatial Sums completed: time " <<
1547  watch.getTimeInSeconds() << "secs" << std::endl;
1548  */
1549  // Unpack into separate momentum and correlator
1550  aop.dilutions(i,j).mom_projs.resize(phases.numMom());
1551 
1552  for(int mom_num = 0 ; mom_num < phases.numMom() ; ++mom_num)
1553  {
1554  aop.dilutions(i,j).mom_projs[mom_num].mom = phases.numToMom(mom_num);
1555 
1556  aop.dilutions(i,j).mom_projs[mom_num].op = a_sum[mom_num];
1557  }
1558  } // end for j
1559  } // end for i
1560  swiss.stop();
1561 
1562 
1563  QDPIO::cout << "Annihilation operator construction: operator= " << l
1564  << " time= "
1565  << swiss.getTimeInSeconds()
1566  << " secs" << std::endl;
1567 
1568  QDPIO::cout << "Annihilation op testval( t0 = " <<
1569  participating_timeslices[t0] << ") = " <<
1570  annih_oper.time_slices[0].dilutions(0,0).mom_projs[0].op[0]
1571  << std::endl;
1572 
1573  //Hard code the elemental op name for now
1574  std::stringstream cnvrt;
1575  cnvrt << annih_oper.id << "_t" << participating_timeslices[t0] << "_snk.lime";
1576 
1577  std::string filename;
1578 
1579  filename = cnvrt.str();
1580 
1581  // Write the meta-data and the binary for this operator
1582  swiss.reset();
1583  swiss.start();
1584  {
1585  XMLBufferWriter src_record_xml, file_xml;
1586  BinaryBufferWriter src_record_bin;
1587 
1588  push(file_xml, "SinkMesonOperator");
1589  write(file_xml, "Params", params.param);
1590  write(file_xml, "Config_info", gauge_xml);
1591  write(file_xml, "Op_Info",qqq_oplist.ops[l]);
1592  push(file_xml, "QuarkSources");
1593 
1594  const int n0 = 0;
1595  const int n1 = 1;
1596 
1597  push(file_xml, "Quark_l");
1598  push(file_xml, "TimeSlice");
1599  push(file_xml, "Dilutions");
1600  for (int dil = 0; dil < diluted_quarks[n0]->getDilSize(t0) ; ++dil)
1601  {
1602  write(file_xml, "elem", diluted_quarks[n0]->getSourceHeader(t0, dil));
1603  }
1604  pop(file_xml); //dilutions
1605  pop(file_xml); //TimeSlice
1606  pop(file_xml); //Quark_l
1607 
1608  push(file_xml, "Quark_r");
1609  push(file_xml, "TimeSlice");
1610  push(file_xml, "Dilutions");
1611  for (int dil = 0; dil < diluted_quarks[n1]->getDilSize(t0) ; ++dil)
1612  {
1613  write(file_xml, "elem", diluted_quarks[n1]->getSourceHeader(t0, dil));
1614  }
1615  pop(file_xml); //dilutions
1616  pop(file_xml); //TimeSlice
1617  pop(file_xml); //Quark_r
1618 
1619  pop(file_xml);//QuarkSources
1620  push(file_xml, "QuarkSinks");
1621 
1622  push(file_xml, "Quark_l");
1623  write(file_xml, "PropHeader", diluted_quarks[n0]->getPropHeader(0,0));
1624  pop(file_xml);
1625 
1626  push(file_xml, "Quark_r");
1627  write(file_xml, "PropHeader", diluted_quarks[n1]->getPropHeader(0,0));
1628  pop(file_xml);
1629 
1630  pop(file_xml);//QuarkSinks
1631  pop(file_xml);//SinkMesonOperator
1632 
1633  QDPFileWriter qdp_file(file_xml, filename, // are there one or two files???
1634  QDPIO_SINGLEFILE, QDPIO_SERIAL, QDPIO_OPEN);
1635 
1636 
1637  write(src_record_xml, "MesonAnnihilationOperator", annih_oper);
1638  write(src_record_bin, annih_oper);
1639 
1640  write(qdp_file, src_record_xml, src_record_bin);
1641 
1642  }
1643  swiss.stop();
1644 
1645  QDPIO::cout << "Annihilation Operator writing: operator = " << l
1646  << " time= " << swiss.getTimeInSeconds() << " secs" << std::endl;
1647 
1648  } // end for l (operator )
1649 
1650  } //End Make annihilation operator
1651 
1652  } //end t0
1653 
1654  // Close the namelist output file XMLDAT
1655  pop(xml_out); // StochMeson
1656 
1657  snoop.stop();
1658  QDPIO::cout << InlineStochGroupMesonEnv::name << ": total time = "
1659  << snoop.getTimeInSeconds()
1660  << " secs" << std::endl;
1661 
1662  QDPIO::cout << InlineStochGroupMesonEnv::name << ": ran successfully" << std::endl;
1663 
1664  END_CODE();
1665  } // func
1666 
1667  } // namespace InlineStochGroupMesonEnv
1668 
1669  /*! @} */ // end of group hadron
1670 
1671 } // namespace Chroma
Inline measurement factory.
Class for counted reference semantics.
Definition: handle.h:33
void operator()(const unsigned long update_no, XMLWriter &xml_out)
Do the measurement.
void func(const unsigned long update_no, XMLWriter &xml_out)
Do the measurement.
multi1d< std::map< KeySmearedDispColorVector_t, SmearedDispColorVector_t > > disp_soln_maps
virtual const LatticeFermion & smearSource(int qnum, const KeySmearedQuark_t &key)
Smear sources and solutions.
virtual const multi1d< LatticeComplex > & displaceObject(std::map< KeySmearedDispColorVector_t, SmearedDispColorVector_t > &disp_quark_map, const KeySmearedDispColorVector_t &key, const LatticeFermion &smrd_q)
Displace an object.
multi1d< std::map< KeySmearedQuark_t, SmearedQuark_t > > smeared_soln_maps
multi1d< std::map< KeySmearedDispColorVector_t, SmearedDispColorVector_t > > disp_src_maps
Maps of smeared displaced color vectors.
virtual multi1d< LatticeComplex > getDispSolution(int quark_num, const KeySmearedDispColorVector_t &key)
Accessor.
SmearedDispObjects(int disp_length, multi1d< Handle< DilutionScheme< LatticeFermion > > > dil_quarks, Handle< QuarkSmearing< LatticeFermion > > qsmr, const multi1d< LatticeColorMatrix > &u_smr)
Constructor from smeared std::map.
const multi1d< LatticeColorMatrix > & u
Gauge field.
multi1d< Handle< DilutionScheme< LatticeFermion > > > diluted_quarks
virtual const LatticeFermion & smearSolution(int qnum, const KeySmearedQuark_t &key)
virtual multi1d< LatticeComplex > getDispSource(int quark_num, const KeySmearedDispColorVector_t &key)
Accessor.
Handle< QuarkSmearing< LatticeFermion > > quarkSmearing
multi1d< std::map< KeySmearedQuark_t, SmearedQuark_t > > smeared_src_maps
Maps of smeared color vectors.
Base class for quark smearing.
Fourier transform phase factor support.
Definition: sftmom.h:35
int numSubsets() const
Number of subsets - length in decay direction.
Definition: sftmom.h:63
multi1d< int > numToMom(int mom_num) const
Convert momenta id to actual array of momenta.
Definition: sftmom.h:78
multi2d< DComplex > sft(const LatticeComplex &cf) const
Do a sumMulti(cf*phases,getSet())
Definition: sftmom.cc:524
int numMom() const
Number of momenta.
Definition: sftmom.h:60
const Set & getSet() const
The set to be used in sumMulti.
Definition: sftmom.h:57
static T & Instance()
Definition: singleton.h:432
All dilution scheme factories.
Factory for dilution schemes.
Basis rotation matrix from Dirac to Degrand-Rossi (and reverse)
Parallel transport a lattice field.
SpinMatrixD DiracToDRMat()
The Dirac to Degrand-Rossi spin transformation matrix (and reverse)
Definition: diractodr.cc:22
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.
multi1d< GroupXML_t > readXMLArrayGroup(XMLReader &xml_in, const std::string &path, const std::string &type_name)
Read group and return as a std::string.
MesonOpType
Meson operator contraction orderings.
GroupXML_t readXMLGroup(XMLReader &xml_in, const std::string &path, const std::string &type_name)
Read group and return as a std::string.
@ MESON_OP_TYPE_SOLUTION_SOURCE
@ MESON_OP_TYPE_SOURCE_SOLUTION
@ MESON_OP_TYPE_SOLUTION_SOLUTION
@ MESON_OP_TYPE_SOURCE_SOURCE
LatticePropagator displacement(const multi1d< LatticeColorMatrix > &u, const LatticePropagator &chi, int length, int dir)
Apply a displacement operator to a lattice field.
Class for counted reference semantics.
Inline measurement of stochastic group meson operators.
unsigned j
Definition: ldumul_w.cc:35
unsigned n
Definition: ldumul_w.cc:36
void savern(int iseed[4])
Definition: make_seeds.cc:46
Make xml file writer.
Named object function std::map.
static bool registered
Local registration flag.
bool registerAll()
Register all the factories.
void readOps(TwoQuarkOps_t &oplist, const std::string &ops_file)
Read 2-quark operators file, assign correct displacement length.
const int N_quarks
Number of quarks to be used in this construction.
bool registerAll()
Register all the factories.
void read(XMLReader &xml, const std::string &path, InlineStochGroupMesonEnv::Params::Param_t &param)
void makeColorSinglet(LatticeComplex &singlet, const multi1d< LatticeComplex > &q0, const multi1d< LatticeComplex > &q1, const Subset &subset)
bool operator<(const KeySmearedQuark_t &a, const KeySmearedQuark_t &b)
Support for the keys of smeared quarks.
SpinMatrix rotate_mat(adj(DiracToDRMat()))
void write(XMLWriter &xml, const std::string &path, const InlineStochGroupMesonEnv::Params::Param_t &param)
multi2d< DComplex > contractOp(SmearedDispObjects &smrd_disp_vecs, int n0, const KeySmearedDispColorVector_t &k0, int n1, const KeySmearedDispColorVector_t &k1, MesonOpType contractType, const SftMom &phases, int t0)
bool registerAll()
Register all the factories.
bool registerAll()
Register all the factories.
Handle< FermBC< LatticeStaggeredFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > > reader(XMLReader &xml_in, const std::string &path)
Helper function for the FermionAction readers.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
QDP::StandardOutputStream & operator<<(QDP::StandardOutputStream &s, const multi1d< int > &d)
Definition: npr_vertex_w.cc:12
static multi1d< LatticeColorMatrix > u
push(xml_out,"Condensates")
Double c
Definition: invbicg.cc:108
int i
Definition: pbg5p_w.cc:55
Complex a
Definition: invbicg.cc:95
pop(xml_out)
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
DComplex d
Definition: invbicg.cc:99
START_CODE()
Complex b
Definition: invbicg.cc:96
::std::string string
Definition: gtest.h:1979
int l
Definition: pade_trln_w.cc:111
Print out basic info about this program.
All quark smearing constructors.
Factory for producing quark smearing objects.
Fourier transform phase factor support.
All make sink constructors.
Factory for producing quark prop sinks.
All source smearing.
Factory for producing quark smearing objects.
Hold group xml and type id.
Meson operator time slices corresponding to location of operator source.
void writeXML(XMLWriter &xml_out, const std::string &path)