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