CHROMA
inline_prop_matelem_pt_colorvec_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Compute the matrix element of LatticeColorVector*M^-1*LatticeColorVector
3  *
4  * Propagator calculation on a colorstd::vector
5  */
6 
7 #include "fermact.h"
8 //#include "meas/inline/hadron/inline_prop_matelem_pt_colorvec_w.h"
11 #include "meas/glue/mesplq.h"
13 #include "qdp_map_obj.h"
16 #include "util/ferm/key_val_db.h"
17 #include "util/ft/sftmom.h"
18 #include "util/info/proginfo.h"
20 
22 #include "util/ferm/transf.h"
23 #include "io/qprop_io.h"
24 
25 namespace Chroma
26 {
27  namespace InlinePropMatElemPtColorVecEnv
28  {
29  //! Propagator input
31  {
32  XMLReader inputtop(xml, path);
33 
34  read(inputtop, "gauge_id", input.gauge_id);
35  read(inputtop, "colorvec_id", input.colorvec_id);
36  read(inputtop, "prop_id", input.prop_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 InlinePropMatElemPtColorVecEnv::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_id", input.prop_id);
48  write(xml, "prop_op_file", input.prop_op_file);
49 
50  pop(xml);
51  }
52 
53 
54  //! Propagator input
55  void read(XMLReader& xml, const std::string& path, InlinePropMatElemPtColorVecEnv::Params::Param_t& input)
56  {
57  XMLReader inputtop(xml, path);
58 
59  read(inputtop, "num_vecs", input.num_vecs);
60  read(inputtop, "mass_label", input.mass_label);
61  }
62 
63  //! Propagator output
64  void write(XMLWriter& xml, const std::string& path, const InlinePropMatElemPtColorVecEnv::Params::Param_t& input)
65  {
66  push(xml, path);
67 
68  write(xml, "num_vecs", input.num_vecs);
69  write(xml, "mass_label", input.mass_label);
70 
71  pop(xml);
72  }
73 
74 
75  //! Propagator input
76  void read(XMLReader& xml, const std::string& path, InlinePropMatElemPtColorVecEnv::Params& input)
77  {
79  input = tmp;
80  }
81 
82  //! Propagator output
83  void write(XMLWriter& xml, const std::string& path, const InlinePropMatElemPtColorVecEnv::Params& input)
84  {
85  push(xml, path);
86 
87  write(xml, "Param", input.param);
88  write(xml, "NamedObject", input.named_obj);
89 
90  pop(xml);
91  }
92  } // namespace InlinePropMatElemPtColorVecEnv
93 
94 
95  namespace InlinePropMatElemPtColorVecEnv
96  {
97  namespace
98  {
99  AbsInlineMeasurement* createMeasurement(XMLReader& xml_in,
100  const std::string& path)
101  {
102  return new InlineMeas(Params(xml_in, path));
103  }
104 
105  //! Local registration flag
106  bool registered = false;
107  }
108 
109  const std::string name = "PROP_MATELEM_PT_COLORVEC";
110 
111  //! Register all the factories
112  bool registerAll()
113  {
114  bool success = true;
115  if (! registered)
116  {
117  success &= TheInlineMeasurementFactory::Instance().registerObject(name, createMeasurement);
118  registered = true;
119  }
120  return success;
121  }
122 
123 
124  //----------------------------------------------------------------------------
125  // Param stuff
127 
128  Params::Params(XMLReader& xml_in, const std::string& path)
129  {
130  try
131  {
132  XMLReader paramtop(xml_in, path);
133 
134  if (paramtop.count("Frequency") == 1)
135  read(paramtop, "Frequency", frequency);
136  else
137  frequency = 1;
138 
139  // Parameters for source construction
140  read(paramtop, "Param", param);
141 
142  // Read in the output propagator/source configuration info
143  read(paramtop, "NamedObject", named_obj);
144 
145  // Possible alternate XML file pattern
146  if (paramtop.count("xml_file") != 0)
147  {
148  read(paramtop, "xml_file", xml_file);
149  }
150  }
151  catch(const std::string& e)
152  {
153  QDPIO::cerr << __func__ << ": Caught Exception reading XML: " << e << std::endl;
154  QDP_abort(1);
155  }
156  }
157 
158 
159 
160  // Function call
161  void
162  InlineMeas::operator()(unsigned long update_no,
163  XMLWriter& xml_out)
164  {
165  // If xml file not empty, then use alternate
166  if (params.xml_file != "")
167  {
168  std::string xml_file = makeXMLFileName(params.xml_file, update_no);
169 
170  push(xml_out, "PropMatElemPtColorVec");
171  write(xml_out, "update_no", update_no);
172  write(xml_out, "xml_file", xml_file);
173  pop(xml_out);
174 
175  XMLFileWriter xml(xml_file);
176  func(update_no, xml);
177  }
178  else
179  {
180  func(update_no, xml_out);
181  }
182  }
183 
184 
185  // Real work done here
186  void
187  InlineMeas::func(unsigned long update_no,
188  XMLWriter& xml_out)
189  {
190  START_CODE();
191 
192  StopWatch snoop;
193  snoop.reset();
194  snoop.start();
195 
196  // Test and grab a reference to the gauge field
197  multi1d<LatticeColorMatrix> u;
198  XMLBufferWriter gauge_xml;
199  try
200  {
201  u = TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
202  TheNamedObjMap::Instance().get(params.named_obj.gauge_id).getRecordXML(gauge_xml);
203  }
204  catch( std::bad_cast )
205  {
206  QDPIO::cerr << name << ": caught dynamic cast error" << std::endl;
207  QDP_abort(1);
208  }
209  catch (const std::string& e)
210  {
211  QDPIO::cerr << name << ": std::map call failed: " << e << std::endl;
212  QDP_abort(1);
213  }
214 
215  push(xml_out, "PropMatElemPtColorVec");
216  write(xml_out, "update_no", update_no);
217 
218  QDPIO::cout << name << ": propagator colorstd::vector matrix element calculation" << std::endl;
219 
220  proginfo(xml_out); // Print out basic program info
221 
222  // Write out the input
223  write(xml_out, "Input", params);
224 
225  // Write out the config header
226  write(xml_out, "Config_info", gauge_xml);
227 
228  push(xml_out, "Output_version");
229  write(xml_out, "out_version", 1);
230  pop(xml_out);
231 
232  // Calculate some gauge invariant observables just for info.
233  MesPlq(xml_out, "Observables", u);
234 
235  //
236  // Read in the source along with relevant information.
237  //
238  XMLReader source_file_xml, source_record_xml;
239  XMLReader prop_file_xml, prop_record_xml;
240 
241  QDPIO::cout << "Snarf the source from a named buffer. Check for the prop std::map" << std::endl;
242  try
243  {
244  // NB We are just checking this is here.
246  TheNamedObjMap::Instance().getData<LatticePropagator>(params.named_obj.prop_id);
247 
248  // Snarf the source info. This is will throw if the colorvec_id is not there
249  TheNamedObjMap::Instance().get(params.named_obj.colorvec_id).getFileXML(source_file_xml);
250  TheNamedObjMap::Instance().get(params.named_obj.colorvec_id).getRecordXML(source_record_xml);
251 
252  TheNamedObjMap::Instance().get(params.named_obj.prop_id).getFileXML(prop_file_xml);
253  TheNamedObjMap::Instance().get(params.named_obj.prop_id).getRecordXML(prop_record_xml);
254 
255 
256  // Write out the source header
257  write(xml_out, "Source_file_info", source_file_xml);
258  write(xml_out, "Source_record_info", source_record_xml);
259 
260  write(xml_out, "Propagator_file_info", prop_file_xml);
261  write(xml_out, "Propagator_record_info", prop_record_xml);
262  }
263  catch (std::bad_cast)
264  {
265  QDPIO::cerr << name << ": caught dynamic cast error" << std::endl;
266  QDP_abort(1);
267  }
268  catch (const std::string& e)
269  {
270  QDPIO::cerr << name << ": error extracting source_header or prop std::map: " << e << std::endl;
271  QDP_abort(1);
272  }
273 
274  // Cast should be valid now
275  QDP::MapObject<int,EVPair<LatticeColorVector> >& eigen_source =
277 
278  // Cast should be valid now
279  const LatticePropagator& prop =
280  TheNamedObjMap::Instance().getData<LatticePropagator>(params.named_obj.prop_id);
281 
282  PropSourceConst_t source_header;
283 
284  read(prop_record_xml, "/Propagator/PropSource", source_header);
285 
286 
287  QDPIO::cout << "Point Propagator successfully found and parsed" << std::endl;
288 
289  //multi1d<int> source = prop_header.source_header.getTSrce();
290  multi1d<int> source = source_header.getTSrce();
291 
292  int decay_dir = source_header.j_decay ;
293  int t_source = source_header.t_source;
294 
295  //Real Mass = getMass(source_header.fermact);
296  //QDPIO::cout << "Mass = " << Mass << std::endl;
297  QDPIO::cout << "source = "
298  << source[0]<<" "
299  << source[1]<<" "
300  << source[2]<<" "
301  << source[3]<< std::endl;
302 
303 
304  // Sanity check - write out the norm2 of the source in the Nd-1 direction
305  // Use this for any possible verification
306  {
307  // Initialize the slow Fourier transform phases
308  SftMom phases(0, true, Nd-1);
309 
310  EVPair<LatticeColorVector> tmpvec; eigen_source.get(0,tmpvec);
311  multi1d<Double> source_corrs = sumMulti(localNorm2(tmpvec.eigenVector), phases.getSet());
312 
313  push(xml_out, "Source_correlators");
314  write(xml_out, "source_corrs", source_corrs);
315  pop(xml_out);
316  }
317 
318  // Another sanity check
319  if (params.param.num_vecs > eigen_source.size())
320  {
321  QDPIO::cerr << __func__ << ": num_vecs= " << params.param.num_vecs
322  << " is greater than the number of available colorvectors= "
323  << eigen_source.size() << std::endl;
324  QDP_abort(1);
325  }
326 
327 
328  //
329  // DB storage
330  //
331  BinaryStoreDB< SerialDBKey<KeyPropElementalOperator_t>, SerialDBData<ValPropElementalOperator_t> > qdp_db;
332 
333  // Open the file, and write the meta-data and the binary for this operator
334  if (! qdp_db.fileExists(params.named_obj.prop_op_file))
335  {
336  XMLBufferWriter file_xml;
337 
338  push(file_xml, "DBMetaData");
339  write(file_xml, "id", std::string("propElemOp"));
340  write(file_xml, "lattSize", QDP::Layout::lattSize());
341  write(file_xml, "decay_dir", decay_dir);
342  write(file_xml, "source", source);
343  proginfo(file_xml); // Print out basic program info
344  write(file_xml, "Params", params.param);
345  write(file_xml, "Config_info", gauge_xml);
346  write(file_xml, "Weights", getEigenValues(eigen_source, params.param.num_vecs));
347  pop(file_xml);
348 
349  std::string file_str(file_xml.str());
350  qdp_db.setMaxUserInfoLen(file_str.size());
351 
352  qdp_db.open(params.named_obj.prop_op_file, O_RDWR | O_CREAT, 0664);
353 
354  qdp_db.insertUserdata(file_str);
355  }
356  else
357  {
358  qdp_db.open(params.named_obj.prop_op_file, O_RDWR, 0664);
359  }
360 
361 
362  //
363  // Try the factories
364  //
365  try
366  {
367  StopWatch swatch;
368  swatch.reset();
369  swatch.start();
370 
371  //
372  // Loop over the source color and spin, creating the source
373  // and calling the relevant propagator routines.
374  //
375  const int num_vecs = params.param.num_vecs;
376  const int decay_dir = decay_dir;
377 
378  // Initialize the slow Fourier transform phases
379  SftMom phases(0, true, decay_dir);
380 
381  // Binary output
382  push(xml_out, "ColorVecMatElems");
383 
384  // Loop over each operator
385  QDPIO::cout << "t_source = " << t_source << std::endl;
386 
387  //
388  // All the loops
389  //
390  // NOTE: I pull the spin source and sink loops to the outside intentionally.
391  // The idea is to create a colorstd::vector index (2d) array. These are not
392  // too big, but are big enough to make the IO efficient, and the DB efficient
393  // on reading. For N=32 and Lt=128, the mats are 2MB.
394  //
395  for(int spin_source=0; spin_source < Ns; ++spin_source){
396  QDPIO::cout << "spin_source = " << spin_source << std::endl;
397  for(int spin_sink=0; spin_sink < Ns; ++spin_sink){
398  QDPIO::cout << "spin_sink = " << spin_sink << std::endl;
399 
400  // Invert the time - make it an independent key
401  multi1d<KeyValPropElementalOperator_t> buf(phases.numSubsets());
402  for(int t=0; t < phases.numSubsets(); ++t){
403  buf[t].key.key().t_slice = t;
404  buf[t].key.key().t_source = t_source;
405  buf[t].key.key().spin_src = spin_source;
406  buf[t].key.key().spin_snk = spin_sink;
407  buf[t].key.key().mass_label = params.param.mass_label;
408  buf[t].val.data().mat.resize(num_vecs,Nc);
409  }
410 
411  for(int colorvec_source=0; colorvec_source < Nc; ++colorvec_source){
412  LatticeFermion q ;
413  PropToFerm(prop, q, colorvec_source, spin_source);
414  LatticeColorVector vec_source(peekSpin(q, spin_sink));
415 
416  for(int colorvec_sink=0; colorvec_sink < num_vecs; ++colorvec_sink){
417  EVPair<LatticeColorVector> vec_sink; eigen_source.get(colorvec_sink,vec_sink);
418 
419  multi1d<ComplexD> hsum(sumMulti(localInnerProduct(vec_sink.eigenVector, vec_source), phases.getSet()));
420 
421  for(int t=0; t < hsum.size(); ++t){
422  buf[t].val.data().mat(colorvec_sink,colorvec_source) = hsum[t];
423  }
424 
425  } // for colorvec_sink
426  } // for colorvec_source
427 
428  QDPIO::cout << "insert: spin_source= " << spin_source << " spin_sink= " << spin_sink << std::endl;
429  for(int t=0; t < phases.numSubsets(); ++t){
430  qdp_db.insert(buf[t].key, buf[t].val);
431  }
432 
433  } // for spin_sink
434  } // for spin_source
435 
436  pop(xml_out);
437 
438  swatch.stop();
439  QDPIO::cout << "Propagators computed: time= "
440  << swatch.getTimeInSeconds()
441  << " secs" << std::endl;
442  }
443  catch (const std::string& e)
444  {
445  QDPIO::cout << name << ": caught exception around contractions: " << e << std::endl;
446  QDP_abort(1);
447  }
448 
449  pop(xml_out); // prop_matelem_colorvec
450 
451  snoop.stop();
452  QDPIO::cout << name << ": total time = "
453  << snoop.getTimeInSeconds()
454  << " secs" << std::endl;
455 
456  QDPIO::cout << name << ": ran successfully" << std::endl;
457 
458  END_CODE();
459  }
460  }//namespace Inline
461 
462 } // 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
int numSubsets() const
Number of subsets - length in decay direction.
Definition: sftmom.h:63
const Set & getSet() const
The set to be used in sumMulti.
Definition: sftmom.h:57
static T & Instance()
Definition: singleton.h:432
Class structure for fermion actions.
void PropToFerm(const LatticePropagatorF &b, LatticeFermionF &a, int color_index, int spin_index)
Extract a LatticeFermion from a LatticePropagator.
Definition: transf.cc:226
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.
std::map< std::string, SinkPropContainer_t > prop
MODS_t & eigen_source
Eigenvectors.
Compute the matrix element of LatticeColorVector*M^-1*LatticeColorVector.
Key for propagator colorstd::vector sources.
Key for propagator colorstd::vector matrix elements.
Key and values for DB.
Make xml file writer.
int t
Definition: meslate.cc:37
Nd
Definition: meslate.cc:74
Double q
Definition: mesq.cc:17
Named object function std::map.
static bool registered
Local registration flag.
QDPSubTypeTrait< typename BinaryReturn< C1, C2, FnLocalInnerProduct >::Type_t >::Type_t localInnerProduct(const QDPSubType< T1, C1 > &l, const QDPType< T2, C2 > &r)
void write(XMLWriter &xml, const std::string &path, const InlinePropMatElemPtColorVecEnv::Params::NamedObject_t &input)
Propagator output.
void read(XMLReader &xml, const std::string &path, InlinePropMatElemPtColorVecEnv::Params::NamedObject_t &input)
Propagator input.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
static multi1d< LatticeColorMatrix > u
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
push(xml_out,"Condensates")
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()
::std::string string
Definition: gtest.h:1979
Print out basic info about this program.
Routines associated with Chroma propagator IO.
Fourier transform phase factor support.
A Pair type.
struct Chroma::InlinePropMatElemPtColorVecEnv::Params::NamedObject_t named_obj
struct Chroma::InlinePropMatElemPtColorVecEnv::Params::Param_t param
Propagator source construction parameters.
Definition: qprop_io.h:27
Holds of vectors and weights.