CHROMA
inline_disco_eo_eigcg_w.cc
Go to the documentation of this file.
1 
2 /*! \file
3  * \brief Inline measurement 3pt_prop
4  *
5  * This uses eo-preconditioning on the noise std::vector side as well as
6  * with the eig-cg vectors. NOTE THIS CODE DOES NOT WORK.
7  */
8 
9 
10 #include "handle.h"
19 #include "meas/sources/zN_src.h"
25 #include "meas/hadron/barQll_w.h"
28 #include "meas/glue/mesplq.h"
29 #include "util/ft/sftmom.h"
30 #include "util/info/proginfo.h"
32 #include "util/info/unique_id.h"
33 #include "util/ferm/transf.h"
35 
38 
41 
45 
46 #include "util/ferm/key_val_db.h"
47 
48 #include <qdp-lapack.h>
49 #include <qdp_config.h>
50 
51 #include <vector>
52 #include <map>
53 
54 namespace Chroma{
55  namespace InlineDiscoEoEigCGEnv{
56  namespace{
57  AbsInlineMeasurement* createMeasurement(XMLReader& xml_in,
58  const std::string& path)
59  {
60  return new InlineMeas(Params(xml_in, path));
61  }
62 
63  //! Local registration flag
64  bool registered = false;
65  }
66 
67  const std::string name = "DISCO_EO_EIGCG";
68 
69  //! Register all the factories
70  bool registerAll()
71  {
72  bool success = true;
73  if (! registered)
74  {
75  //success &= BaryonOperatorEnv::registerAll();
76  success &= TheInlineMeasurementFactory::Instance().registerObject(name, createMeasurement);
77  registered = true;
78  }
79  return success;
80  }
81 
82  // Reader for input parameters
83  void read(XMLReader& xml, const std::string& path, Params::Param_t& param){
84  XMLReader paramtop(xml, path);
85 
86  int version;
87  read(paramtop, "version", version);
88 
89  switch (version)
90  {
91  case 1:
92  /************************************************************/
93  read(paramtop,"max_path_length",param.max_path_length);
94  read(paramtop,"p2_max",param.p2_max);
95  read(paramtop,"mass_label",param.mass_label);
96  param.chi = readXMLArrayGroup(paramtop, "Quarks", "DilutionType");
97  param.action = readXMLGroup(paramtop, "FermionAction","FermAct");
98 
99  break;
100 
101  default :
102  /**************************************************************/
103 
104  QDPIO::cerr << "Input parameter version " << version << " unsupported." << std::endl;
105  QDP_abort(1);
106  }
107  }
108 
109  // Writer for input parameters
110  void write(XMLWriter& xml, const std::string& path, const Params::Param_t& param){
111  push(xml, path);
112 
113  int version = 1;
114 
115  write(xml, "version", version);
116 
117  write(xml,"max_path_length",param.max_path_length);
118  write(xml,"p2_max",param.p2_max);
119  write(xml,"mass_label",param.mass_label);
120  push(xml,"FermionAction");
121  xml<<param.action.xml ;
122  pop(xml);
123  push(xml,"Quarks");
124  for( int t(0);t<param.chi.size();t++){
125  push(xml,"elem");
126  xml<<param.chi[t].xml;
127  pop(xml);
128  }
129  pop(xml);
130 
131  pop(xml); // final pop
132  }
133 
134  //! Gauge field parameters
135  void read(XMLReader& xml, const std::string& path, Params::NamedObject_t& input)
136  {
137  XMLReader inputtop(xml, path);
138 
139  read(inputtop, "gauge_id" , input.gauge_id ) ;
140  // If no file is included, then we will turn off the eig_cg part (ie, use 0 vectors)
141  if ( inputtop.count("evecs_file")!=0 )
142  read(inputtop, "evecs_file" , input.evecs_file ) ;
143  else
144  input.evecs_file="";
145 
146  read(inputtop, "op_db_file" , input.op_db_file ) ;
147  }
148 
149  //! Gauge field parameters
150  void write(XMLWriter& xml, const std::string& path, const Params::NamedObject_t& input){
151  push(xml, path);
152 
153  write(xml, "gauge_id" , input.gauge_id );
154  write(xml, "evecs_file" , input.evecs_file );
155  write(xml, "op_db_file" , input.op_db_file );
156  pop(xml);
157  }
158 
159 
160  // Param stuff
162  frequency = 0;
163  }
164 
165  Params::Params(XMLReader& xml_in, const std::string& path)
166  {
167  try
168  {
169  XMLReader paramtop(xml_in, path);
170 
171  if (paramtop.count("Frequency") == 1)
172  read(paramtop, "Frequency", frequency);
173  else
174  frequency = 1;
175 
176  // Read program parameters
177  read(paramtop, "Param", param);
178 
179  // Read in the output propagator/source configuration info
180  read(paramtop, "NamedObject", named_obj);
181 
182  // Possible alternate XML file pattern
183  if (paramtop.count("xml_file") != 0)
184  {
185  read(paramtop, "xml_file", xml_file);
186  }
187  }
188  catch(const std::string& e)
189  {
190  QDPIO::cerr << __func__ << ": Caught Exception reading XML: " << e << std::endl;
191  QDP_abort(1);
192  }
193  }
194 
195 
196  void Params::write(XMLWriter& xml_out, const std::string& path)
197  {
198  push(xml_out, path);
199 
200  // Parameters for source construction
201  InlineDiscoEoEigCGEnv::write(xml_out, "Param", param);
202 
203  // Write out the output propagator/source configuration info
204  InlineDiscoEoEigCGEnv::write(xml_out, "NamedObject", named_obj);
205 
206  pop(xml_out);
207  }
208 
209 
210  //! Meson operator
212  {
213  multi1d<short int> mom ; /*!< D-1 momentum of this operator */
214  unsigned short int t_slice ; /*!< Meson operator time slice */
215  multi1d<short int> disp ; /*!< Displacement dirs of quark (right)*/
216 
218  mom.resize(Nd-1);
219  }
220  };
221 
222  bool operator<(const KeyOperator_t& a, const KeyOperator_t& b){
223  return ((a.t_slice<b.t_slice)||(a.mom<b.mom)||(a.disp<b.disp));
224  }
225 
226  std::ostream& operator<<(std::ostream& os, const KeyOperator_t& d)
227  {
228  os << "KeyOperator_t:"
229  << " t_slice = " << d.t_slice
230  << ", disp = ";
231  for (int i=0; i<d.disp.size();i++)
232  os << d.disp[i] << " " ;
233 
234  os << ", mom = ";
235  for (int i=0; i<d.mom.size();i++)
236  os << d.mom[i] << " " ;
237 
238  os << std::endl;
239 
240  return os;
241  }
242 
244  public:
245  multi1d<ComplexD> op ;
246  ValOperator_t(){op.resize(Ns*Ns);} // Here go the 16 gamma matrices
248  } ;
249 
250  //-------------------------------------------------------------------------
251  //! stream IO
252  std::ostream& operator<<(std::ostream& os, const ValOperator_t& d)
253  {
254  os << "ValOperator_t:\n";
255  for (int i=0; i<d.op.size();i++)
256  os <<" gamma["<<i<<"] = "<< d.op[i] << std::endl ;
257 
258  return os;
259  }
260 
261  struct KeyVal_t{
264  };
265 
266  //! KeyOperator reader
267  void read(BinaryReader& bin, KeyOperator_t& d){
268  read(bin,d.t_slice);
269  unsigned short int n ;
270  read(bin,n);
271  d.disp.resize(n);
272  read(bin,d.disp);
273  d.mom.resize(Nd-1) ;
274  read(bin,d.mom);
275  }
276  //! KeyOperator writer
277  void write(BinaryWriter& bin, const KeyOperator_t& d){
278  write(bin,d.t_slice);
279  unsigned short int n ;
280  n = d.disp.size();
281  write(bin,n);
282  write(bin,d.disp);
283  write(bin,d.mom);
284  }
285 
286  //! ValOperator reader
287  void read(BinaryReader& bin, ValOperator_t& d){
288  d.op.resize(Ns*Ns);
289  read(bin,d.op);
290  }
291  //! ValOperator writer
292  void write(BinaryWriter& bin, const ValOperator_t& d){
293  write(bin,d.op);
294  }
295 
296  namespace{
297  StandardOutputStream& operator<<(StandardOutputStream& os, const multi1d<short int>& d){
298  if (d.size() > 0){
299  os << d[0];
300  for(int i=1; i < d.size(); i++)
301  os << " " << d[i];
302  }
303  return os;
304  }
305  }
306 
307  typedef LatticeFermion T;
308  typedef multi1d<LatticeColorMatrix> P;
309  typedef multi1d<LatticeColorMatrix> Q;
310 
312  createOddOdd_Op( const Params::Param_t& param, const P& u){
313  //
314  // Initialize fermion action
315  //
316  std::istringstream xml_s(param.action.xml);
317  XMLReader fermacttop(xml_s);
318  QDPIO::cout << "FermAct = " << param.action.id << std::endl;
319  //
320  // Try the factories
321  //
323  try{
324  QDPIO::cout << "Try the various Wilson fermion factories" << std::endl;
325  // Generic Wilson-Type stuff
327  Sf(TheFermionActionFactory::Instance().createObject(param.action.id,
328  fermacttop,
329  param.action.path));
330  state = Sf->createState(u);//got the state
331  QDPIO::cout << "Suitable factory found: compute the trace quark prop"<<std::endl;
332  }
333  catch (const std::string& e){
334  QDPIO::cout << name
335  << ": caught exception instantiating the action: " << e << std::endl;
336  }
337 
338  // Now need to construct the operator
339  // this seems tricky and maybe there is a better way... (Robert-Balint help!)
340  // We only work with Wilson and Clover Wilson fermions so check the following 2
341  // For the moment we restrict to EO preconditioned only
342 
343  //this is the Odd-Odd piece
345  if(param.action.id == "WILSON"){
346  WilsonFermActParams wp(fermacttop,param.action.path);
347  return new EvenOddPrecWilsonLinOp(state,wp.Mass,wp.anisoParam ) ;
348  }
349  else if(param.action.id == "CLOVER"){
350  CloverFermActParams cp(fermacttop,param.action.path);
351  return new EvenOddPrecCloverLinOp(state,cp) ;
352  }
353  else{
354  QDPIO::cout<<name<<" : Tough luck dude! No code for you..."<<std::endl ;
355  QDP_abort(1);
356  }
357  return NULL ;
358  }
359 
361  multi1d<Real> evals ;
362  multi1d<Complex> H ;
363  multi1d<Complex> HU ;
364  int ldh ;
365  int Nvec ;
366  } ;
367 
368  void do_disco(std::map< KeyOperator_t, ValOperator_t >& db,
369  const LatticeFermion& qbar,
370  const LatticeFermion& q,
371  const SftMom& p,
372  const int& t,
373  const Subset& trb,
374  const multi1d<short int>& path,
375  const int& max_path_length ){
376  QDPIO::cout<<" Computing Operator with path length "<<path.size()
377  <<" on timeslice "<<t<<". Path: "<<path <<std::endl;
378  /**
379  This do_disco routine should not be called, I don't think it's ever needed...
380 
381  **/
382  ValOperator_t val ;
383  KeyOperator_t key ;
384  std::pair<KeyOperator_t, ValOperator_t> kv ;
385  kv.first.t_slice = t ;
386  if(path.size()==0){
387  kv.first.disp.resize(1);
388  kv.first.disp[0] = 0 ;
389  }
390  else
391  kv.first.disp = path ;
392 
393  multi1d< multi1d<ComplexD> > foo(p.numMom()) ;
394  for (int m(0); m < p.numMom(); m++)
395  foo[m].resize(Ns*Ns);
396  for(int g(0);g<Ns*Ns;g++){
397  LatticeComplex cc = localInnerProduct(qbar,Gamma(g)*q);
398  for (int m(0); m < p.numMom(); m++){
399  // trb is the set of even/odd sites on timeslice t
400  foo[m][g] = sum(p[m]*cc,trb);
401  }
402  }
403  for (int m(0); m < p.numMom(); m++){
404  for(int i(0);i<(Nd-1);i++)
405  kv.first.mom[i] = p.numToMom(m)[i] ;
406 
407  kv.second.op = foo[m];
408  std::pair<std::map< KeyOperator_t, ValOperator_t >::iterator, bool> itbo;
409 
410  itbo = db.insert(kv);
411  if( itbo.second ){
412  QDPIO::cout<<"Inserting new entry in std::map\n";
413  }
414  else{ // if insert fails, key already exists, so add result
415  for(int i(0);i<kv.second.op.size();i++)
416  itbo.first->second.op[i] += kv.second.op[i] ;
417  }
418  }//m
419 
420  if(path.size()<max_path_length){
421  QDPIO::cout<<" attempt to add new path. "
422  <<" current path length is : "<<path.size();
423  multi1d<short int> new_path(path.size()+1);
424  QDPIO::cout<<" new path length is : "<<new_path.size()<<std::endl;
425  for(int i(0);i<path.size();i++)
426  new_path[i] = path[i] ;
427  for(int sign(-1);sign<2;sign+=2)
428  for(int mu(0);mu<Nd;mu++){
429  new_path[path.size()]= sign*(mu+1) ;
430  //skip back tracking
431  bool back_track=false ;
432  if(path.size()>0)
433  if(path[path.size()-1] == -new_path[path.size()])
434  back_track=true;
435  if(!back_track){
436  QDPIO::cout<<" Added path: "<<new_path<<std::endl;
437  LatticeFermion q_mu ;
438  if(sign>0)
439  q_mu = shift(q, FORWARD, mu);
440  else
441  q_mu = shift(q, BACKWARD, mu);
442 
443  do_disco(db, qbar, q_mu, p, t, trb, new_path, max_path_length);
444  } // skip backtracking
445  } // mu
446  }// path.size loop
447 
448  }// do_disco
449 
450  void do_disco(std::map< KeyOperator_t, ValOperator_t >& db,
451  const LatticeFermion& qbar,
452  const LatticeFermion& q,
454  const SftMom& p,
455  const int& t,
456  const Subset& trb,
457  const multi1d<short int>& path,
458  const int& max_path_length ){
459  QDPIO::cout<<" Computing Operator with path length "<<path.size()
460  <<" on timeslice "<<t<<". Path: "<<path <<std::endl;
461 
462  /** This version of do_disco is the one we call
463 
464  Currently we have editted this so as to put the timeslice and phase factor
465  on the correct part for the "even" sum, so right now this code is not terribly
466  efficient. Work on optimizing this all later.
467 
468  **/
469 
470  ValOperator_t val ;
471  KeyOperator_t key ;
472  std::pair<KeyOperator_t, ValOperator_t> kv ;
473  kv.first.t_slice = t ;
474  if(path.size()==0){
475  kv.first.disp.resize(1);
476  kv.first.disp[0] = 0 ;
477  }
478  else
479  kv.first.disp = path ;
480 
481  multi1d< multi1d<ComplexD> > foo(p.numMom()) ;
482  for (int m(0); m < p.numMom(); m++)
483  foo[m].resize(Ns*Ns);
484 
485  // First do the odd sites, and put in foo.
486  for(int g(0);g<Ns*Ns;g++){
487  LatticeComplex cc = localInnerProduct(qbar,Gamma(g)*q);
488  for (int m(0); m < p.numMom(); m++){
489  // trb is the set of odd sites on timeslice t
490  foo[m][g] = sum(p[m]*cc,trb);
491  }
492  }
493 
494  Set t0set;
495  t0set.make(TimeSliceFunc(Nd-1));// 3 is decay dir....need to change this later.
496 
497  // Now we have to do the even sites
498  LatticeFermion qtmp = zero ;
499  LatticeFermion qbar2 = zero ;
500  Doo->evenOddLinOp(qtmp,q,PLUS);
501  Doo->evenEvenInvLinOp(qbar2,qtmp,PLUS);
502  for(int g(0);g<Ns*Ns;g++){
503  for (int m(0); m < p.numMom(); m++){
504  LatticeFermion q2 = zero ;
505  LatticeFermion chi = zero ;
506  qtmp = zero ;
507  q2[t0set[t]] = p[m]*(Gamma(g)*qbar2);//rhs should only set the t0 terms non-zero...
508  Doo->evenEvenInvLinOp(qtmp,q2,PLUS);
509  Doo->oddEvenLinOp(chi,qtmp,PLUS);
510  foo[m][g] += innerProduct(qbar,chi);
511  }
512  }
513 
514  for (int m(0); m < p.numMom(); m++){
515  for(int i(0);i<(Nd-1);i++)
516  kv.first.mom[i] = p.numToMom(m)[i] ;
517 
518  kv.second.op = foo[m];
519  std::pair<std::map< KeyOperator_t, ValOperator_t >::iterator, bool> itbo;
520 
521  itbo = db.insert(kv);
522  if( itbo.second ){
523  QDPIO::cout<<"Inserting new entry in std::map\n";
524  }
525  else{ // if insert fails, key already exists, so add result
526  for(int i(0);i<kv.second.op.size();i++)
527  itbo.first->second.op[i] += kv.second.op[i] ;
528  }
529  }//m
530 
531  if(path.size()<max_path_length){
532  QDPIO::cout<<" attempt to add new path. "
533  <<" current path length is : "<<path.size();
534  multi1d<short int> new_path(path.size()+1);
535  QDPIO::cout<<" new path length is : "<<new_path.size()<<std::endl;
536  for(int i(0);i<path.size();i++)
537  new_path[i] = path[i] ;
538  for(int sign(-1);sign<2;sign+=2)
539  for(int mu(0);mu<Nd;mu++){
540  new_path[path.size()]= sign*(mu+1) ;
541  //skip back tracking
542  bool back_track=false ;
543  if(path.size()>0)
544  if(path[path.size()-1] == -new_path[path.size()])
545  back_track=true;
546  if(!back_track){
547  QDPIO::cout<<" Added path: "<<new_path<<std::endl;
548  LatticeFermion q_mu ;
549  if(sign>0)
550  q_mu = shift(q, FORWARD, mu);
551  else
552  q_mu = shift(q, BACKWARD, mu);
553 
554  do_disco(db, qbar, q_mu, Doo, p, t, trb, new_path, max_path_length);
555  } // skip backtracking
556  } // mu
557  }// path.size loop
558 
559  }// do_disco
560 
561  void do_disco(std::map< KeyOperator_t, ValOperator_t >& db,
562  const Params::Param_t& param,
563  const P& u,
564  const SftMom& p,
565  const int& t,
566  const Subset& trb,
567  const multi1d<short int>& path){
568  QDPIO::cout<<" Computing Operator with path length "<<path.size()
569  <<" on timeslice "<<t<<". Path: "<<path <<std::endl;
570 
571  int max_path_length = param.max_path_length;
572  ValOperator_t val ;
573  KeyOperator_t key ;
574  std::pair<KeyOperator_t, ValOperator_t> kv ;
575  kv.first.t_slice = t ;
576  if(path.size()==0){
577  kv.first.disp.resize(1);
578  kv.first.disp[0] = 0 ;
579  }
580  else
581  kv.first.disp = path ;
582 
584 
585  multi1d< multi1d<ComplexD> > foo(p.numMom()) ;
586  for (int m(0); m < p.numMom(); m++){
587  foo[m].resize(Ns*Ns);
588  for (int g(0); g<Ns*Ns;g++){
589  foo[m][g] = 0.0;
590  }
591  }
592  for(int g(0);g<Ns*Ns;g++){
593  for(int col=0;col<3;col++){
594  for(int sp=0;sp<4;sp++){
595  Fermion tt = zero ;
596  ColorVector cv = zero ;
597  Complex z = cmplx(Real(1.0),0.0) ;
598  pokeColor(cv,z,col);
599  pokeSpin(tt,cv,sp);
600  LatticeFermion V = tt ;
601  for (int m(0); m < p.numMom(); m++){
602  //only on even sites
603  LatticeFermion DV = zero;
604  Doo->evenEvenInvLinOp(DV,V,PLUS);
605  foo[m][g] += sum(p[m]*localInnerProduct(V,Gamma(g)*DV),trb);
606  }//m
607  }//spin
608  }//col
609  }//g
610 
611  for (int m(0); m < p.numMom(); m++){
612  for(int i(0);i<(Nd-1);i++)
613  kv.first.mom[i] = p.numToMom(m)[i] ;
614 
615  kv.second.op = foo[m];
616  std::pair<std::map< KeyOperator_t, ValOperator_t >::iterator, bool> itbo;
617 
618  itbo = db.insert(kv);
619  if( itbo.second ){
620  QDPIO::cout<<"Inserting new entry in std::map\n";
621  }
622  else{ // if insert fails, key already exists, so add result
623  for(int i(0);i<kv.second.op.size();i++)
624  itbo.first->second.op[i] += kv.second.op[i] ;
625  }
626  }//m
627 
628  }// do_disco
629 
630  void do_disco(std::map< KeyOperator_t, ValOperator_t >& db,
631  CholeskyFactors Clsk ,
632  multi1d<LatticeFermion>& vec,
633  const Params::Param_t& param,
635  const P& u,
636  const int& t,
637  const Subset& trb,
638  const SftMom& p,
639  const multi1d<short int>& path){
640 
641  QDPIO::cout<<" Computing Operator with path length "<<path.size()
642  <<" on timeslice "<<t<<". Path: "<<path <<std::endl;
643  int max_path_length = param.max_path_length;
644 
645  ValOperator_t val ;
646  KeyOperator_t key ;
647  std::pair<KeyOperator_t, ValOperator_t> kv ;
648 
649  kv.first.t_slice = t ;
650  if(path.size()==0){
651  kv.first.disp.resize(1);
652  kv.first.disp[0] = 0 ;
653  }
654  else
655  kv.first.disp = path ;
656 
657  int ldb = vec.size();
658  int info;
659  char U = 'U';
660  int Nrhs = ldb; // Because we have no dilution vectors, but rhs's are made of EigCG vecs...
661  multi2d<Complex> B(Nrhs,ldb);
662  multi1d< multi1d<LatticeFermion> > DDvec(ldb);
663  // Zero out the fermions
664  for(int i(0); i<ldb;i++){
665  DDvec[i].resize(Ns*Ns);
666  DDvec[i] = zero;
667  }
668  std::cout<<"We've initialized DDvec!\n";
669 
670  for(int i(0); i<ldb;i++){
671  LatticeFermion qtmp = zero ;
672  LatticeFermion qtmp2 = zero ;
673  LatticeFermion qtmp3 = zero ;
674  Doo->evenOddLinOp(qtmp,vec[i],PLUS);
675  Doo->evenEvenInvLinOp(qtmp2,qtmp,PLUS);
676  for(int g(0);g<Ns*Ns;g++){
677  Doo->evenEvenInvLinOp(qtmp3,Gamma(g)*qtmp2,PLUS);
678  qtmp = zero;
679  qtmp2 = zero;
680  Doo->oddEvenLinOp(qtmp,qtmp3,PLUS);
681  Doo->oddOddLinOp(qtmp2,qtmp,MINUS);
682  qtmp = zero;
683  Doo->oddOddLinOp(qtmp,Gamma(g)*vec[i],MINUS);
684  // I think this could run into memory problems with too many vectors
685  DDvec[i][g] = qtmp + qtmp2;
686  }
687  }
688 
689  multi1d< multi1d<ComplexD> > foo(p.numMom()) ;
690  for (int m(0); m < p.numMom(); m++){
691  foo[m].resize(Ns*Ns);
692  foo[m] = zero;
693  }
694  // Okay, this is probably inefficient, b/c we are doing some things multiple times...
695  for(int g(0);g<Ns*Ns;g++){
696  for (int m(0); m < p.numMom(); m++){
697  for (int i(0); i<ldb;i++){
698  for (int j(0); j<ldb;j++){
699  // LatticeComplex cc1 = localInnerProduct(Svec[j],Gamma(g)*vec[i]);
700  LatticeComplex cc2 = localInnerProduct(vec[j],DDvec[i][g]);
701  B[i][j] = sum(p[m]*cc2,trb) ;
702  }//j
703  }//i
704 
705 #if BASE_PRECISION == 32
706  int r = QDPLapack::cpotrs(U, Clsk.Nvec, Nrhs, Clsk.HU, Clsk.ldh, B, ldb, info);
707 #else
708  int r = QDPLapack::zpotrs(U, Clsk.Nvec, Nrhs, Clsk.HU, Clsk.ldh, B, ldb, info);
709 #endif
710  // and at this point, B is H^-1 Vdag Sdag gamma V, so
711  // For debugging purposes, both r and info should be zero...
712  QDPIO::cout<<"do_disco cpotrs r = "<<r<<std::endl;
713  QDPIO::cout<<"do_disco cpotrs info = "<<info<<std::endl;
714 
715  // foo[m][g] = 0.0;
716  for (int i(0); i<ldb;i++){
717  foo[m][g] += ComplexD(B[i][i]);// Trace over last set of indices
718  }
719  }//m
720  }//g
721  for (int m(0); m < p.numMom(); m++){
722  for(int i(0);i<(Nd-1);i++)
723  kv.first.mom[i] = p.numToMom(m)[i] ;
724 
725  kv.second.op = foo[m];
726 
727  std::pair<std::map< KeyOperator_t, ValOperator_t >::iterator, bool> itbo;
728 
729  itbo = db.insert(kv);
730  if( !(itbo.second) ){ // if insert fails, key already exists, so add result
731  for(int i(0);i<kv.second.op.size();i++)
732  itbo.first->second.op[i] += kv.second.op[i] ;
733  }
734  else
735  QDPIO::cout<<"Inserting new entry in std::map\n";
736  }
737 
738  if(path.size()<max_path_length){
739  QDPIO::cout<<" attempt to add new path. "
740  <<" current path length is : "<<path.size();
741  multi1d<short int> new_path(path.size()+1);
742  QDPIO::cout<<" new path length is : "<<new_path.size()<<std::endl;
743  for(int i(0);i<path.size();i++)
744  new_path[i] = path[i] ;
745  for(int sign(-1);sign<2;sign+=2)
746  for(int mu(0);mu<Nd;mu++){
747  new_path[path.size()]= sign*(mu+1) ;
748  //skip back tracking
749  bool back_track=false ;
750  if(path.size()>0)
751  if(path[path.size()-1] == -new_path[path.size()])
752  back_track=true;
753  if(!back_track){
754  QDPIO::cout<<" Added path: "<<new_path<<std::endl;
755  multi1d<LatticeFermion> vec_mu(vec.size()) ;
756  for(int j(0);j<vec.size();j++){
757  if(sign>0)
758  vec_mu[j] = shift(vec[j], FORWARD, mu);
759  else
760  vec_mu[j] = shift(vec[j], BACKWARD, mu);
761  }
762  do_disco(db, Clsk , vec_mu, param, Doo, u, t, trb, p, new_path);
763  } // skip backtracking
764  } // mu
765  }// path.size loop
766  }// do_disco
767 
768 
769  void ReadOPTEigCGVecs(multi1d<LatticeFermion>& vec,
770  CholeskyFactors& Clsk,
771  const std::string& evecs_file)
772  {
773  QDPIO::cout<<name<<" : Reading vecs from "
774  << evecs_file <<std::endl ;
775  StopWatch swatch;
776  swatch.reset();
777  swatch.start();
778 
779  int Nvecs,ldh ;
780  // File XML
781  XMLReader file_xml;
782  // Open file
783  QDPFileReader to(file_xml,evecs_file,QDPIO_SERIAL);
784  read(file_xml, "/OptEigInfo/ncurEvals", Nvecs);
785  read(file_xml, "/OptEigInfo/ldh", ldh);
786  // Added Nvecs and ldh to the Cholesky structure for later...
787  Clsk.Nvec = Nvecs;
788  Clsk.ldh = ldh;
789  vec.resize(Nvecs);
790  Clsk.evals.resize(ldh);
791  Clsk.H.resize(ldh*ldh);
792  Clsk.HU.resize(ldh*ldh);
793 
794  for(int v(0);v<Nvecs;v++){
795  XMLReader record_xml;
796  read(to, record_xml, vec[v]);
797  }
798 
799  XMLReader record_xml;
800  read(to, record_xml, Clsk.evals);
801  read(to, record_xml, Clsk.H);
802  read(to, record_xml, Clsk.HU);
803  swatch.stop();
804  QDPIO::cout<<name<<" : Time to read vecs= "
805  << swatch.getTimeInSeconds() <<" secs "<<std::endl ;
806  }
807 
808  // PRchi returns chitilde = (1 - V Hinv Vdag Sdag S)chi given
809  // an input chi std::vector, and of course the vectors and H.
810  void PRchi(multi1d<multi1d< multi1d<LatticeFermion> > >& quarkstilde,
811  const multi1d< Handle< DilutionScheme<LatticeFermion> > >& quarks,
813  CholeskyFactors Clsk , multi1d<LatticeFermion>& vec,
814  const Params::Param_t& param, const P& u){
815  char U = 'U';
816  int info;
817 
818  int ldb = vec.size();//This is the offset that will for now be the size of each std::vector
819 
820  Set timerb;
821  timerb.make(TimeSliceRBFunc(3));
822 
823  // Now we have to create the Sdag * S * quarks object to put into B
824  for(int n(0);n<quarks.size();n++){
825  for (int it(0) ; it < quarks[n]->getNumTimeSlices() ; ++it){
826  int t = quarks[n]->getT0(it) ;
827  int Nrhs = quarks[n]->getDilSize(it) ;
828  multi2d<Complex> B(Nrhs, ldb);
829  for(int i(0); i<ldb;i++){
830  for(int j = 0 ; j < Nrhs ; j++){
831  /**
832  Here, we are calculating
833  Sdag S chi
834  where chi (the solution) is Sinv eta, and eta is the noise std::vector (the source)
835  Thus, we can save time by using the source, and calculate
836  Sdag eta,
837  and that's what is being done here
838  **/
839  LatticeFermion qsrc = quarks[n]->dilutedSource(it,j);
840  LatticeFermion SdagSchi = zero ;
841  Doo->oddOddLinOp(SdagSchi,qsrc,MINUS);
842  // Sum over (odd) lattice sites, so we return a complex number for B
843  // SHOULDN'T THIS BE ONLY SITES THAT INCLUDE THE TIMESLICE WE'RE ON?
844  // B[j][i] = sum(localInnerProduct(vec[i],SdagSchi),rb[1]);
845  B[j][i] = sum(localInnerProduct(vec[i],SdagSchi),timerb[2*t+1]);
846  }
847  }
848 
849  // This could be better done.
850 #if BASE_PRECISION == 32
851  int r = QDPLapack::cpotrs(U, Clsk.Nvec, Nrhs, Clsk.HU, Clsk.ldh, B, ldb, info);
852 #else
853  int r = QDPLapack::zpotrs(U, Clsk.Nvec, Nrhs, Clsk.HU, Clsk.ldh, B, ldb, info);
854 #endif
855  QDPIO::cout<<"PRchi cpotrs r = "<<r<<std::endl;
856  QDPIO::cout<<"PRchi cpotrs info = "<<info<<std::endl;
857 
858  for(int j = 0 ; j < Nrhs ; j++){
859  LatticeFermion vB = zero; //B[j][0]*vec[0];
860  LatticeFermion q = quarks[n]->dilutedSolution(it,j);
861  for(int i(0); i<ldb;i++)
862  vB += B[j][i]*vec[i];
863  quarkstilde[n][it][j] = q - vB;
864  std::cout<<"Norm of PRchi: "<<norm2(quarkstilde[n][it][j])<<std::endl;
865  }
866  }
867  }
868  }// End of PRchi call...
869 
870  /**
871  This version is called if we have no EigCG vectors, so PR = 1
872  **/
873  void PRchi(multi1d<multi1d< multi1d<LatticeFermion> > >& quarkstilde,
874  multi1d< Handle< DilutionScheme<LatticeFermion> > >& quarks){
875 
876  for(int n(0);n<quarks.size();n++){
877  for (int it(0) ; it < quarks[n]->getNumTimeSlices() ; it++){
878  for (int j(0) ; j < quarks[n]->getDilSize(it) ; j++){
879  quarkstilde[n][it][j] = quarks[n]->dilutedSolution(it,j);
880  }
881  }
882  }
883  }// End of PRchi call...
884 
885  //--------------------------------------------------------------
886  // Function call
887  // void
888  //InlineDiscoEigCG::operator()(unsigned long update_no,
889  // XMLWriter& xml_out)
890  void InlineMeas::operator()(unsigned long update_no,
891  XMLWriter& xml_out)
892  {
893  // If xml file not empty, then use alternate
894  if (params.xml_file != ""){
895  std::string xml_file = makeXMLFileName(params.xml_file, update_no);
896 
897  push(xml_out, "discoEigCG");
898  write(xml_out, "update_no", update_no);
899  write(xml_out, "xml_file", xml_file);
900  pop(xml_out);
901 
902  XMLFileWriter xml(xml_file);
903  func(update_no, xml);
904  }
905  else{
906  func(update_no, xml_out);
907  }
908  }
909 
910 
911  // Function call
912  //void
913  //InlineDiscoEigCG::func(unsigned long update_no,
914  // XMLWriter& xml_out)
915  void InlineMeas::func(unsigned long update_no,
916  XMLWriter& xml_out)
917  {
918  START_CODE();
919 
920  StopWatch snoop;
921  snoop.reset();
922  snoop.start();
923 
924  /*** Stuff to set up everything! ***/
925  // Test and grab a reference to the gauge field
926  XMLBufferWriter gauge_xml;
927  try
928  {
929  TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
930  TheNamedObjMap::Instance().get(params.named_obj.gauge_id).getRecordXML(gauge_xml);
931  }
932  catch( std::bad_cast )
933  {
934  QDPIO::cerr << InlineDiscoEoEigCGEnv::name << ": caught dynamic cast error"
935  << std::endl;
936  QDP_abort(1);
937  }
938  catch (const std::string& e)
939  {
940  QDPIO::cerr << name << ": std::map call failed: " << e
941  << std::endl;
942  QDP_abort(1);
943  }
944  const multi1d<LatticeColorMatrix>& u =
945  TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
946 
947  push(xml_out, "discoEigCG");
948  write(xml_out, "update_no", update_no);
949 
950  QDPIO::cout << name << ": Disconnected diagrams with eigCG vectors" << std::endl;
951 
952  proginfo(xml_out); // Print out basic program info
953 
954  // Write out the input
955  params.write(xml_out, "Input");
956 
957  // Write out the config info
958  write(xml_out, "Config_info", gauge_xml);
959 
960  push(xml_out, "Output_version");
961  write(xml_out, "out_version", 1);
962  pop(xml_out);
963 
964  // First calculate some gauge invariant observables just for info.
965  // This is really cheap.
966  MesPlq(xml_out, "Observables", u);
967 
968  // Save current seed
969  Seed ran_seed;
970  QDP::RNG::savern(ran_seed);
971 
972  //
973  // Read the source and solutions
974  //
975  StopWatch swatch;
976  swatch.reset();
977  swatch.start();
978 
979  int N_quarks = params.param.chi.size() ;
980 
981  multi1d< Handle< DilutionScheme<LatticeFermion> > > quarks(N_quarks);
982 
983  try{
984  // Loop over quark dilutions
985  for(int n(0); n < params.param.chi.size(); n++){
986  const GroupXML_t& dil_xml = params.param.chi[n];
987  std::istringstream xml_d(dil_xml.xml);
988  XMLReader diltop(xml_d);
989  QDPIO::cout << "Dilution type = " << dil_xml.id << std::endl;
990  quarks[n] =
991  TheFermDilutionSchemeFactory::Instance().createObject(dil_xml.id,
992  diltop,
993  dil_xml.path);
994  }
995  }
996  catch(const std::string& e){
997  QDPIO::cerr << name << ": Caught Exception constructing dilution scheme: " << e << std::endl;
998  QDP_abort(1);
999  }
1000 
1001 
1002  //-------------------------------------------------------------------
1003  //Sanity checks
1004 
1005  //Another Sanity check, the three quarks must all be
1006  //inverted on the same cfg
1007  for (int n = 1 ; n < N_quarks ; n++){
1008  if (quarks[0]->getCfgInfo() != quarks[n]->getCfgInfo()){
1009  QDPIO::cerr << name << " : Quarks do not contain the same cfg info";
1010  QDPIO::cerr << ", quark "<< n << std::endl;
1011  QDP_abort(1);
1012  }
1013  }
1014 
1015  //Also ensure that the cfg on which the inversions were performed
1016  //is the same as the cfg that we are using
1017  {
1018  std::string cfgInfo;
1019 
1020  //Really ugly way of doing this (Robert Heeeelp!!)
1021  XMLBufferWriter top;
1022  write(top, "Config_info", gauge_xml);
1023  XMLReader from(top);
1024  XMLReader from2(from, "/Config_info");
1025  std::ostringstream os;
1026  from2.print(os);
1027 
1028  cfgInfo = os.str();
1029 
1030  if (cfgInfo != quarks[0]->getCfgInfo()){
1031  QDPIO::cerr << name << " : Quarks do not contain the same";
1032  QDPIO::cerr << " cfg info as the gauge field." ;
1033  QDPIO::cerr << "gauge: XX"<<cfgInfo<<"XX quarks: XX" ;
1034  QDPIO::cerr << quarks[0]->getCfgInfo()<<"XX"<< std::endl;
1035  QDP_abort(1);
1036  }
1037  }
1038 
1039  int decay_dir = quarks[0]->getDecayDir();
1040  //
1041  // Initialize the slow Fourier transform phases
1042  //
1043  SftMom phases(params.param.p2_max, false, decay_dir);
1044 
1045  // Create even/odd subsets on each timeslice:
1046  Set timerb;
1047  timerb.make(TimeSliceRBFunc(decay_dir));
1048 
1049  // Another sanity check. The seeds of all the quarks must be different
1050  // and thier decay directions must be the same
1051  for(int n = 1 ; n < quarks.size(); n++){
1052  if(toBool(quarks[n]->getSeed()==quarks[0]->getSeed())){
1053  QDPIO::cerr << name << ": error, quark seeds are the same" << std::endl;
1054  QDP_abort(1);
1055  }
1056 
1057  if(toBool(quarks[n]->getDecayDir()!=quarks[0]->getDecayDir())){
1058  QDPIO::cerr<<name<< ": error, quark decay dirs do not match" <<std::endl;
1059  QDP_abort(1);
1060  }
1061  }
1062 
1063  // Instantiate the Dirac operator:
1065 
1066  //Now read evecs from disk, if file is specified.
1067  multi1d<LatticeFermion> vec; // the vectors
1068  CholeskyFactors Clsk; // the Cholesky Factors
1069  if (params.named_obj.evecs_file!="")
1071 
1072  std::map< KeyOperator_t, ValOperator_t > data ;
1073 
1074  multi1d<multi1d< multi1d<LatticeFermion> > > quarkstilde(quarks.size());
1075  for(int n(0);n<quarks.size();n++){
1076  quarkstilde[n].resize(quarks[n]->getNumTimeSlices());
1077  for (int it(0) ; it < quarks[n]->getNumTimeSlices() ; it++){
1078  quarkstilde[n][it].resize(quarks[n]->getDilSize(it));
1079  }
1080  }
1081 
1082  /*** Everything is initialized, now to the real code ***/
1083 
1084  // chitilde = P_R S^-1 \eta_o
1085  // If there is no evecs file, then PR = 1
1086  if (params.named_obj.evecs_file!="")
1087  PRchi(quarkstilde, quarks, Doo, Clsk, vec, params.param, u);
1088  else
1089  PRchi(quarkstilde, quarks);
1090 
1091  for(int n(0);n<quarks.size();n++){
1092  for (int it(0) ; it < quarks[n]->getNumTimeSlices() ; it++){
1093  int t = quarks[n]->getT0(it) ;
1094  QDPIO::cout<<" Doing quark: "<< n <<std::endl ;
1095  QDPIO::cout<<" quark: "<< n <<" has "<<quarks[n]->getDilSize(it);
1096  QDPIO::cout<<" dilutions on time slice "<< t <<std::endl ;
1097  for(int i = 0 ; i < quarks[n]->getDilSize(it) ; i++){
1098  QDPIO::cout<<" Doing dilution : "<<i<<std::endl ;
1099  multi1d<short int> d ;
1100  // Now, I want to do the two trace terms that have the chitilde. Thus
1101  // The q chosen should be from quarkstilde = PR*quarks->dilutedSolution()
1102  LatticeFermion qbar = quarks[n]->dilutedSource(it,i);
1103  LatticeFermion q = quarkstilde[n][it][i];
1104  QDPIO::cout<<" Starting recursion "<<std::endl ;
1105  // this does both terms with the noise vectors
1106  // the 1 means to do sum over odd sites...
1107  do_disco(data, qbar, q, Doo, phases, t, timerb[2*t+1], d, params.param.max_path_length);
1108 
1109  QDPIO::cout<<" done with recursion! "
1110  <<" The length of the path is: "<<d.size()<<std::endl ;
1111  }
1112  QDPIO::cout<<" Done with dilutions for quark: "<<n <<std::endl ;
1113  }
1114  }
1115 
1116  // Okay, first we are going to normalize all the pieces above by the number of
1117  // quarks...
1120  std::map< KeyOperator_t, ValOperator_t >::iterator it;
1121 
1122  for(it=data.begin();it!=data.end();it++){
1123  key.key() = it->first ;
1124  val.data().op.resize(it->second.op.size()) ;
1125  for(int i(0);i<it->second.op.size();i++){
1126  val.data().op[i] = it->second.op[i]/toDouble(quarks.size());
1127  }
1128  }
1129 
1130  // Now calculate the other pieces which need no normalization...
1131  // Note using just timeslices for quarks[0]...may be a problem
1132  // if quarks dilute on diff timeslices...
1133 
1134  for (int it(0) ; it < quarks[0]->getNumTimeSlices() ; it++){
1135  multi1d<short int> d ;
1136  int t = quarks[0]->getT0(it) ;
1137  // Now let's do the Tr[Dee^-1 gamma]
1138  do_disco(data,params.param, u, phases, t, timerb[2*t+0], d);
1139  }
1140 
1141  if (params.named_obj.evecs_file!=""){// Only if we have vectors...
1142  for (int it(0) ; it < quarks[0]->getNumTimeSlices() ; it++){
1143  multi1d<short int> d ;
1144  int t = quarks[0]->getT0(it) ;
1145  QDPIO::cout<<" Starting recursion again"<<std::endl ;
1146  // Part with deflation..
1147  do_disco(data, Clsk, vec, params.param, Doo, u, t, timerb[2*t+1], phases,d);
1148 
1149  QDPIO::cout<<" done with recursion! "
1150  <<" The length of the path is: "<<d.size()<<std::endl ;
1151  }
1152  }
1153 
1154  //After all the pieces are computed we write the final result to the database
1155  // DB storage
1156  BinaryStoreDB<SerialDBKey<KeyOperator_t>,SerialDBData<ValOperator_t> > qdp_db;
1157 
1158  // Open the file, and write the meta-data and the binary for this operator
1159  {
1160  XMLBufferWriter file_xml;
1161 
1162  push(file_xml, "DBMetaData");
1163  write(file_xml, "id", std::string("discoEigElemOp"));
1164  write(file_xml, "lattSize", QDP::Layout::lattSize());
1165  write(file_xml, "decay_dir", decay_dir);
1166  write(file_xml, "Params", params.param);
1167  write(file_xml, "Config_info", gauge_xml);
1168  pop(file_xml);
1169 
1170  std::string file_str(file_xml.str());
1171  qdp_db.setMaxUserInfoLen(file_str.size());
1172 
1173  qdp_db.open(params.named_obj.op_db_file, O_RDWR | O_CREAT, 0664);
1174 
1175  qdp_db.insertUserdata(file_str);
1176  }
1177 
1178  // Store all the data
1179  for(it=data.begin();it!=data.end();it++){
1180  key.key() = it->first ;
1181  val.data().op.resize(it->second.op.size()) ;
1182  // DON'T normalize to number of quarks here, because we only do it on a
1183  // certain number of the terms in the trace...done above
1184  for(int i(0);i<it->second.op.size();i++)
1185  val.data().op[i] = it->second.op[i];
1186  qdp_db.insert(key,val) ;
1187  }
1188 
1189  // Close the namelist output file XMLDAT
1190  pop(xml_out); // DiscoEigCG
1191 
1192  snoop.stop();
1193  QDPIO::cout << name << ": total time = "
1194  << snoop.getTimeInSeconds()
1195  << " secs" << std::endl;
1196 
1197  QDPIO::cout << name << ": ran successfully" << std::endl;
1198 
1199  END_CODE();
1200  }
1201  } // namespace InlineDiscoEoEigCGEnv
1202 }// namespace chroma
Inline measurement factory.
Heavy Baryon (Qll) 2-pt function : Orginos and Savage.
Baryon spin and projector matrices.
All baryon operators.
Factory for producing baryon operators.
Even-odd preconditioned Clover-Dirac operator.
Even-odd preconditioned linear operator.
Definition: eoprec_linop.h:92
Even-odd preconditioned Wilson-Dirac operator.
Class for counted reference semantics.
Definition: handle.h:33
Inline measurement of stochastic baryon operators.
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
D & data()
Setter.
Definition: key_val_db.h:78
Serializable key harness.
Definition: key_val_db.h:21
K & key()
Setter.
Definition: key_val_db.h:30
Fourier transform phase factor support.
Definition: sftmom.h:35
static T & Instance()
Definition: singleton.h:432
Function object used for constructing the time-slice set.
Definition: barQll_w.h:95
Parameters for Clover fermion action.
int mu
Definition: cool.cc:24
Random Z(N) source construction using dilution.
All dilution scheme factories.
Factory for dilution schemes.
Parallel transport a lattice field.
Even-odd const determinant Wilson-like fermact.
Fermion action factories.
All Wilson-type 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.
multi1d< GroupXML_t > readXMLArrayGroup(XMLReader &xml_in, const std::string &path, const std::string &type_name)
Read group and return as a std::string.
GroupXML_t readXMLGroup(XMLReader &xml_in, const std::string &path, const std::string &type_name)
Read group and return as a std::string.
Class for counted reference semantics.
Inline measurement of stochastic 3pt functions.
Key and values for DB.
unsigned j
Definition: ldumul_w.cc:35
unsigned n
Definition: ldumul_w.cc:36
Linear operators.
static int m[4]
Definition: make_seeds.cc:16
void savern(int iseed[4])
Definition: make_seeds.cc:46
Make xml file writer.
int z
Definition: meslate.cc:36
int it
Definition: meslate.cc:33
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.
void do_disco(std::map< KeyOperator_t, ValOperator_t > &db, const LatticeFermion &qbar, const LatticeFermion &q, const SftMom &p, const int &t, const Subset &trb, const multi1d< short int > &path, const int &max_path_length)
void PRchi(multi1d< multi1d< multi1d< LatticeFermion > > > &quarkstilde, const multi1d< Handle< DilutionScheme< LatticeFermion > > > &quarks, Handle< EvenOddPrecLinearOperator< T, P, Q > > &Doo, CholeskyFactors Clsk, multi1d< LatticeFermion > &vec, const Params::Param_t &param, const P &u)
Handle< EvenOddPrecLinearOperator< T, P, Q > > createOddOdd_Op(const Params::Param_t &param, const P &u)
void read(XMLReader &xml, const std::string &path, Params::Param_t &param)
void write(XMLWriter &xml, const std::string &path, const Params::Param_t &param)
void ReadOPTEigCGVecs(multi1d< LatticeFermion > &vec, CholeskyFactors &Clsk, const std::string &evecs_file)
multi1d< LatticeColorMatrix > P
std::ostream & operator<<(std::ostream &os, const KeyOperator_t &d)
bool registerAll()
Register all the factories.
multi1d< LatticeColorMatrix > Q
bool operator<(const KeyOperator_t &a, const KeyOperator_t &b)
QDPSubTypeTrait< typename BinaryReturn< C1, C2, FnLocalInnerProduct >::Type_t >::Type_t localInnerProduct(const QDPSubType< T1, C1 > &l, const QDPType< T2, C2 > &r)
BinaryReturn< C1, C2, FnInnerProduct >::Type_t innerProduct(const QDPSubType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
const int N_quarks
Number of quarks to be used in this construction.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
for(cb=0;cb< Ncb;++cb) cp+
Definition: pbg5p_w.cc:145
static multi1d< LatticeColorMatrix > u
push(xml_out,"Condensates")
int i
Definition: pbg5p_w.cc:55
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
Double cp
Definition: invbicg.cc:107
multi1d< LatticeFermion > chi(Ncb)
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
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
static QOP_info_t info
Double zero
Definition: invbicg.cc:106
::std::string string
Definition: gtest.h:1979
#define FORWARD
Definition: primitives.h:82
#define BACKWARD
Definition: primitives.h:83
Print out basic info about this program.
Double sum
Definition: qtopcor.cc:37
All quark smearing constructors.
Factory for producing quark smearing objects.
Fourier transform phase factor support.
SpinMatrix sp
All make sink constructors.
Factory for producing quark prop sinks.
All source smearing.
Factory for producing quark smearing objects.
Params for clover ferm acts.
Hold group xml and type id.
struct Chroma::InlineDiscoEoEigCGEnv::Params::NamedObject_t named_obj
void write(XMLWriter &xml_out, const std::string &path)
struct Chroma::InlineDiscoEoEigCGEnv::Params::Param_t param
Params for wilson ferm acts.
multi1d< LatticeColorMatrix > U
Generate a unique id.
Wilson fermion action parameters.
Volume source of Z(N) noise.