CHROMA
gbaryonall2all.cc
Go to the documentation of this file.
1 /*
2  * \brief Group Theory based Baryon Operator Construction with All2All
3  */
4 #include "chroma.h"
5 #include "chromabase.h"
11 using namespace Chroma;
12 using namespace GroupBaryonOperatorEnv;
13 
14  void write( XMLWriter& xml, const std::string& path,
16  {
17  push( xml, path );
18  write( xml, "BaryonOpElement", input );
19  pop( xml );
20  }
21  void write( XMLWriter& xml, const std::string& path,
23  {
24  push( xml, path );
25  write( xml, "BaryonOpIndex", input );
26  pop( xml );
27  }
28  void write( XMLWriter& xml, const std::string& path,
30  {
31  push( xml, path );
32  write( xml, "BaryonOpIns", input );
33  pop( xml );
34  }
35  void write( XMLWriter& xml, const std::string& path,
37  {
38  push( xml, path );
39  write( xml, "gauge_id", input.gauge_id );
40  //write( xml, "BaryonOp", input.Bops );
41  pop( xml );
42  }
43  //! BaryonOperator header writer
44  void write( XMLWriter& xml, const std::string& path, const BaryonOperator_t& param )
45  {
46  if( path != "." ) push( xml, path );
47  int version = 1;
48  write( xml, "version", version );
49  write( xml, "mom2_max", param.mom2_max );
50  write( xml, "j_decay", param.j_decay );
51  write( xml, "seed_l", param.seed_l );
52  write( xml, "seed_m", param.seed_m );
53  write( xml, "seed_r", param.seed_r );
54  write( xml, "perms", param.quarkOrderings );
55  if( path != "." ) pop( xml );
56  }
57  bool linkageHack(void)
58  {
59  bool foo = true;
60  // Inline Measurements
62  return foo;
63  }
64 //
65 //! Make a group theory based baryon operator with all-to-all quark
66 //! propagators using the dilution method of TrinLat (2004)
67 //
68 /*! \defgroup GroupBaryonOperator
69  * \ingroup main
70  *
71  * Main program for making a group baryon operator
72  */
73 int main( int argc, char *argv[] )
74 {
75  // ==================================
76  // Put the machine into a known state
77  // ==================================
78  Chroma::initialize( &argc, &argv );
79  START_CODE();
80  QDPIO::cout << "GroupBaryonAll2All: Group Baryon Operators with All2All" << std::endl;
81  QDPIO::cout << "Linkage = " << linkageHack() << std::endl;
82  StopWatch swatch;
83  // ===============
84  // Read input data
85  // ===============
86  // Instantiate xml reader
87  XMLReader xml_in( Chroma::getXMLInputFileName() );
88  const std::string path = "/GroupBaryonAll2All";
89  Params params( xml_in, path );
90  // ==================
91  // Setup lattice size
92  // ==================
93  Layout::setLattSize( params.nrow );
94  Layout::create();
95  // ============
96  // XML business
97  // ============
98  XMLFileWriter& xml_out = Chroma::getXMLOutputInstance(); // XMLDAT by default
99  push( xml_out, "GroupBaryonAll2AllRun" );
100  proginfo( xml_out ); // Print out basic program info
101  write( xml_out, "InputCopy", xml_in ); // save a copy of the input
102  xml_out.flush();
103  // ============
104  // Startup gauge
105  // =============
106  XMLReader gauge_file_xml, gauge_xml;
107  gaugeStartup(gauge_file_xml, gauge_xml, params.gaugestuff.u, params.gaugestuff.cfg);
108  // test the config
111  QDPIO::cout<< "Gauge field measurements " << w_plaq <<" "<< s_plaq <<" "<< t_plaq <<" "<< link <<" "<<std::endl;
112 
113  // Smear the gauge field if needed
114  swatch.start();
115  std::istringstream xml_l( params.gaugestuff.link_smearing.xml );
116  XMLReader linktop( xml_l );
117  const std::string link_path = "/LinkSmearing";
118  QDPIO::cout << "Link smearing type = " << params.gaugestuff.link_smearing.id << std::endl;
120  linkSmearing( TheLinkSmearingFactory::Instance().createObject(
122  linktop,
123  link_path ) );
124  ( *linkSmearing ) ( params.gaugestuff.u );
125  swatch.stop();
126  QDPIO::cout << "Gauge links smeared : time = " << swatch.getTimeInSeconds() << " secs" << std::endl;
128  QDPIO::cout<< "gauge meas " << w_plaq <<" "<< s_plaq <<" "<< t_plaq <<" "<< link <<" "<<std::endl;
129  // ====================
130  // qqq's and baryonop's
131  // ====================
132  multi1d<GroupBaryonQQQ> AQQQ; // annihilation
133  multi1d<GroupBaryonQQQ> CQQQ; // creation
134  multi1d<GroupBaryonOp> AOp;
135  multi1d<GroupBaryonOp> COp;
136  // ====================================
137  // Read in the text file produced by
138  // Gen_Input_For_All2All_Baryons2.pl
139  // called "Chroma_Input.txt"
140  // ====================================
141  ReadTextInput( params, AOp, COp, AQQQ, CQQQ );
142  // ===============
143  // Noise Solutions
144  // ===============
145  multi1d<Params::QuarkSourceSolutions_t> quarks( params.qprop.solns.size() );
146  try
147  {
148  swatch.start();
149  for(int n=0; n < quarks.size(); ++n)
150  {
151  QDPIO::cout << "Attempt to read solutions for source number=" << n << std::endl;
152  quarks[n].dilutions.resize( params.qprop.solns[n].soln_file_names.size() );
153  QDPIO::cout << "dilutions.size= " << quarks[n].dilutions.size() << std::endl;
154  for(int i=0; i < quarks[n].dilutions.size(); ++i)
155  {
156  XMLReader file_xml, record_xml;
157  QDPFileReader from(file_xml, params.qprop.solns[n].soln_file_names[i], QDPIO_SERIAL);
158  read(from, record_xml, quarks[n].dilutions[i].soln);
159  close(from);
160  // need the seed
161  read(record_xml, "/Propagator/PropSource/Source/ran_seed", quarks[n].seed);
162  params.dilution[ n ].ran_seed = quarks[ n ].seed;
163  // read headers
164  read(record_xml, "/Propagator/PropSource", quarks[n].dilutions[i].source_header);
165  read(record_xml, "/Propagator/ForwardProp", quarks[n].dilutions[i].prop_header);
166  }
167  }
168  swatch.stop();
169  }
170  catch (const std::string& e)
171  {
172  QDPIO::cerr << "Error extracting headers: " << e << std::endl;
173  QDP_abort(1);
174  }
175  QDPIO::cout << "Read SOLUTION vectors from files : time = " << swatch.getTimeInSeconds() << " secs" << std::endl;
176  // ====================================
177  //
178  // Main Part of Program
179  //
180  // Run through the QQQ combinations and
181  // multiply by appropriate coefficients
182  // to make the various operators
183  //
184  // ====================================
185  multi1d<LatticeComplex> vresult(1);
186  LatticeComplex result;
187 #ifdef REDUCETOTIMEDILUTION
188  multi1d<Params::QuarkSourceSolutions_t> quarks2( params.qprop.solns.size() );
189 #endif
190  // don't average over equivalent momenta
191  SftMom phases(params.mom2_max, true, params.j_decay);
192  multi2d<DComplex> elem( params.Nmomenta, params.nrow[3] );
193 
194 #ifdef MAKE_SINK_OPERATORS
195 
196 #ifdef REDUCETOTIMEDILUTION
197  // =========================================
198  //
199  // Change the dilution scheme to a lower one
200  //
201  // =========================================
202  // From TfSxCfGf ==> TfSxCxGx hard-wired for test
203  try
204  {
205  swatch.start();
206  for(int n=0; n < quarks.size(); ++n)
207  {
208  QDPIO::cout << "Reducing spin-colour dilution levels for quark " << n << std::endl;
209  int indx, indx2;
210  //
211  // Fix these for now ... temporary fix
212  //
213  int NTimeDilution = params.nrow[3];
214  int NSpinDilution = 4;
215  int NColourDilution = 3;
216  int NSpaceDilution = 1;
217  quarks2[n].dilutions.resize( quarks[n].dilutions.size()/(params.NdilReduce) );
218  quarks2[n].seed = quarks[n].seed;
219  for(int i=0; i < (NTimeDilution); ++i)
220  {
221  // fix this
222  quarks2[n].dilutions[i].source_header = quarks[n].dilutions[i*(params.NdilReduce)].source_header;
223  quarks2[n].dilutions[i].prop_header = quarks[n].dilutions[i*(params.NdilReduce)].prop_header;
224  }
225  // compute the indices needed assuming they are
226  // in the order as below, time--colour--spin--space
227  // index = x + Nx * ( s + Ns * ( c + Nc * ( t ) ) )
228  // = t * Nc * Ns * Nx + c * Ns * Nx + s * Nx + x
229  for(int t=0; t < NTimeDilution; ++t)
230  {
231  for(int x=0; x < NSpaceDilution; ++x)
232  {
233  indx2 = x + NSpaceDilution * ( t );
234  LatticeFermion Q=zero;
235  for(int c=0; c < NColourDilution; ++c)
236  {
237  for(int s=0; s < NSpinDilution; ++s)
238  {
239  indx = x + NSpaceDilution * ( s + NSpinDilution * ( c + NColourDilution * ( t ) ) );
240  //indx = x + NSpaceDilution * s + NSpaceDilution * NSpinDilution * c + NSpaceDilution * NSpinDilution * NColourDilution * t;
241  Q += quarks[n].dilutions[ indx ].soln;
242  }
243  }
244  quarks2[n].dilutions[ indx2 ].soln = Q;
245  }
246  }
247  }
248  } // try
249  catch ( const std::string& e )
250  {
251  QDPIO::cerr << ": Caught Exception removing COLOUR and SPIN-dilution: " << e << std::endl;
252  QDP_abort( 1 );
253  }
254  swatch.stop();
255  QDPIO::cout << "COLOUR and SPIN dilution removed: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
256 #endif
257  try
258  { //
259  // Annihilation Operator (plus) first (only one term)
260  //
261  QDPIO::cout << "Making the SINK operators" << std::endl;
262  int ordering = 0;
263  for(int n=0; n < params.Noperators; ++n)
264  {
265 #ifdef REDUCETOTIMEDILUTION
266  AOp[ n ].baryonoperator.seed_l = quarks2[0].seed;
267  AOp[ n ].baryonoperator.seed_m = quarks2[1].seed;
268  AOp[ n ].baryonoperator.seed_r = quarks2[2].seed;
269 #else
270  AOp[ n ].baryonoperator.seed_l = quarks[0].seed;
271  AOp[ n ].baryonoperator.seed_m = quarks[1].seed;
272  AOp[ n ].baryonoperator.seed_r = quarks[2].seed;
273 #endif
274  }
275  swatch.start();
276  // Sink smear all the quarks up-front
277  std::istringstream xml_s( params.sink_smearing.sink.xml );
278  XMLReader sinktop( xml_s );
280  sinkSmearing(
282  params.sink_smearing.sink.id, sinktop,
284 #ifdef REDUCETOTIMEDILUTION
285  for(int n=0; n < quarks2.size(); ++n)
286  {
287  for(int i=0; i < (params.NH[ordering][n]/params.NdilReduce); ++i)
288  {
289  ( *sinkSmearing ) ( quarks2[ n ].dilutions[ i ].soln );
290  }
291  }
292 #else
293  for(int n=0; n < quarks.size(); ++n)
294  {
295  for(int i=0; i < (params.NH[ordering][n]/params.NdilReduce); ++i)
296  {
297  ( *sinkSmearing ) ( quarks[ n ].dilutions[ i ].soln );
298  }
299  }
300 #endif
301  swatch.stop();
302  QDPIO::cout << "SINK smearings done: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
303 
304  swatch.start();
305  for(int qqq=0; qqq < params.NQQQs; ++qqq)
306  {
307  for(int i=0; i < (params.NH[ordering][0]/params.NdilReduce); ++i)
308  for(int j=0; j < (params.NH[ordering][1]/params.NdilReduce); ++j)
309  for(int k=0; k < (params.NH[ordering][2]/params.NdilReduce); ++k)
310  {
311 #ifdef REDUCETOTIMEDILUTION
312  vresult = AQQQ[ qqq ](
313  quarks2[ 0 ].dilutions[ i ].soln,
314  quarks2[ 1 ].dilutions[ j ].soln,
315  quarks2[ 2 ].dilutions[ k ].soln,
316  PLUS
317  );
318 #else
319  vresult = AQQQ[ qqq ](
320  quarks[ 0 ].dilutions[ i ].soln,
321  quarks[ 1 ].dilutions[ j ].soln,
322  quarks[ 2 ].dilutions[ k ].soln,
323  PLUS
324  );
325 #endif
326  elem = phases.sft( vresult[ 0 ] );
327  for(int b=0; b < AQQQ[ qqq ].NBaryonOps; ++b)
328  {
329  for(int p=0; p < params.Nmomenta; ++p)
330  {
331  for(int t=0; t < params.nrow[ 3 ]; ++t)
332  {
333  AQQQ[ qqq ].baryon[ b ]->baryonoperator.orderings[ 0 ].op( i, j, k ).ind[ 0 ].elem( p, t )
334  += ( AQQQ[ qqq ].coef[ b ] * elem( p, t ) );
335  } // t
336  } // p
337  } // Bop_index b
338  } // i,j,k
339  } // qqq
340  swatch.stop();
341  } // try
342  catch ( const std::string& e )
343  {
344  QDPIO::cerr << ": Caught Exception creating sink baryon operator: " << e << std::endl;
345  QDP_abort( 1 );
346  }
347  QDPIO::cout << "SINK operators done: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
348  // Save the operators
349  // ONLY SciDAC output format is supported!
350  swatch.start();
351  for(int b=0; b < params.Noperators; ++b)
352  {
353  XMLBufferWriter file_xml;
354  push( file_xml, params.Names[ b ] );
355  file_xml << params.param.baryon_operator;
356  write( file_xml, "Config_info", gauge_xml );
357  pop( file_xml );
358  std::string output=params.Names[ b ]+"_snk";
359 
360  QDPFileWriter to( file_xml, output,
361  QDPIO_SINGLEFILE, QDPIO_SERIAL, QDPIO_OPEN );
362  /*
363  // Write the scalar data
364  {
365  XMLBufferWriter record_xml;
366  write( record_xml, "SourceBaryonOperator", COp[ b ].baryonoperator );
367  write( to, record_xml, COp[ b ].baryonoperator.serialize() );
368  }
369  */
370  // Write the scalar data
371  {
372  XMLBufferWriter record_xml;
373  write( record_xml, "SinkBaryonOperator", AOp[ b ].baryonoperator );
374  write( to, record_xml, AOp[ b ].baryonoperator.serialize() );
375  }
376  close( to );
377  }
378  swatch.stop();
379  QDPIO::cout << "Operators written: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
380 #endif // end annihilation operators
381 
382 #ifdef MAKE_SOURCE_OPERATORS
383  // =============
384  // Noise Sources ... re-generate from the rng seeds
385  // =============
386  QDPIO::cout << "Re-generating the noise SOURCES" << std::endl;
387  try
388  {
389  swatch.start();
390  int N;
391  for(int n=0; n < quarks.size(); ++n)
392  {
393  bool first = true;
394  LatticeFermion quark_noise; // noisy source on entire lattice
395  for(int i=0; i < quarks[ n ].dilutions.size(); ++i)
396  {
397  std::istringstream xml_s( quarks[ n ].dilutions[ i ].source_header.source.xml );
398  XMLReader sourcetop( xml_s );
399  //QDPIO::cout << "Source = " << quarks[n].dilutions[i].source_header.source.id << std::endl;
400  if ( quarks[ n ].dilutions[ i ].source_header.source.id != DiluteZNQuarkSourceConstEnv::name )
401  {
402  QDPIO::cerr << "Expected source_type = " << DiluteZNQuarkSourceConstEnv::name << std::endl;
403  QDP_abort( 1 );
404  }
405  //QDPIO::cout << "Quark num= " << n << " dilution num= " << i << std::endl;
406  DiluteZNQuarkSourceConstEnv::Params srcParams( sourcetop,
407  quarks[ n ].dilutions[ i ].source_header.source.path );
409  if ( first )
410  {
411  first = false;
412  quarks[ 0 ].j_decay = srcParams.j_decay;
413  // Grab N
414  N = srcParams.N;
415  // Set the seed to desired value
416  quarks[ n ].seed = srcParams.ran_seed;
417  QDP::RNG::setrn( quarks[ n ].seed );
418  // Create the noisy quark source on the entire lattice
419  zN_src( quark_noise, N );
420  }
421  // The seeds must always agree - here the seed is the unique id of the source
422  if ( toBool( srcParams.ran_seed != quarks[ n ].seed ) )
423  {
424  QDPIO::cerr << "quark source=" << n << " dilution=" << i << " seed does not match" << std::endl;
425  QDP_abort( 1 );
426  }
427  // The N's must always agree
428  if ( toBool( srcParams.N != N ) )
429  {
430  QDPIO::cerr << "quark source=" << n << " dilution=" << i << " N does not match" << std::endl;
431  QDP_abort( 1 );
432  }
433  // Use a trick here, create the source and subtract it from the global noisy
434  // Check at the end that the global noisy is zero everywhere.
435  // NOTE: the seed will be set every call
436  quarks[ n ].dilutions[ i ].source = srcConst( params.gaugestuff.u );
437  quark_noise -= quarks[ n ].dilutions[ i ].source;
438  } // end for i
439  Double dcnt = norm2( quark_noise );
440  if ( toDouble( dcnt ) != 0.0 ) // problematic - seems to work with unnormalized sources
441  {
442  QDPIO::cerr << "Noise not saturated by all potential solutions: dcnt=" << dcnt << std::endl;
443  QDP_abort( 1 );
444  }
445  } // end for n
446  swatch.stop();
447  } // end try
448  catch ( const std::string& e )
449  {
450  QDPIO::cerr << ": Caught Exception creating source: " << e << std::endl;
451  QDP_abort( 1 );
452  }
453  QDPIO::cout << "SOURCE vectors reconstructed from Seeds : time = " << swatch.getTimeInSeconds() << " secs" << std::endl;
454 
455 #ifdef REDUCETOTIMEDILUTION
456  // =========================================
457  //
458  // Change the dilution scheme to a lower one
459  //
460  // =========================================
461  // From TfSxCfGf ==> TfSxCxGx hard-wired for test
462  try
463  {
464  swatch.start();
465  for(int n=0; n < quarks.size(); ++n)
466  {
467  QDPIO::cout << "REDUCING spin-colour dilution levels for quark " << n << std::endl;
468  int indx, indx2;
469  //
470  // Fix these for now ... temporary fix
471  //
472  int NTimeDilution = params.nrow[3];
473  int NSpinDilution = 4;
474  int NColourDilution = 3;
475  int NSpaceDilution = 1;
476  quarks2[n].dilutions.resize( quarks[n].dilutions.size()/(params.NdilReduce) );
477  quarks2[n].seed = quarks[n].seed;
478  for(int i=0; i < (NTimeDilution); ++i)
479  {
480  // fix this
481  quarks2[n].dilutions[i].source_header = quarks[n].dilutions[i*(params.NdilReduce)].source_header;
482  quarks2[n].dilutions[i].prop_header = quarks[n].dilutions[i*(params.NdilReduce)].prop_header;
483  }
484  //DiluteZNQuarkSourceConstEnv::Params newsrcParams;
485  //newsrcParams.j_decay = 3;
486  //newsrcParams.ran_seed = quarks[ n ].seed;
487  // compute the indices needed assuming they are
488  // in the order as below, time--colour--spin--space
489  // index = x + Nx * ( s + Ns * ( c + Nc * ( t ) ) )
490  // = t * Nc * Ns * Nx + c * Ns * Nx + s * Nx + x
491  for(int t=0; t < NTimeDilution; ++t)
492  {
493  //newsrcParams.t_source = t;
494  for(int x=0; x < NSpaceDilution; ++x)
495  {
496  indx2 = x + NSpaceDilution * ( t );
497  LatticeFermion Q=zero;
498  for(int c=0; c < NColourDilution; ++c)
499  {
500  for(int s=0; s < NSpinDilution; ++s)
501  {
502  indx = x + NSpaceDilution * ( s + NSpinDilution * ( c + NColourDilution * ( t ) ) );
503  //indx = x + NSpaceDilution * s + NSpaceDilution * NSpinDilution * c + NSpaceDilution * NSpinDilution * NColourDilution * t;
504  Q += quarks[n].dilutions[ indx ].source;
505  }
506  }
507  quarks2[n].dilutions[ indx2 ].source = Q;
508  }
509  }
510  }
511  } // try
512  catch ( const std::string& e )
513  {
514  QDPIO::cerr << ": Caught Exception removing COLOUR and SPIN-dilution: " << e << std::endl;
515  QDP_abort( 1 );
516  }
517  swatch.stop();
518  QDPIO::cout << "COLOUR and SPIN DILUTION REMOVED: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
519 #endif
520  //
521  // setting the seeds for output (baryon_operator file)
522  // should also put the mass of the quarks in here as well
523  //
524  for(int n=0; n < params.Noperators; ++n)
525  {
526 #ifdef REDUCETOTIMEDILUTION
527  COp[ n ].baryonoperator.seed_l = quarks2[0].seed;
528  COp[ n ].baryonoperator.seed_m = quarks2[1].seed;
529  COp[ n ].baryonoperator.seed_r = quarks2[2].seed;
530 #else
531  COp[ n ].baryonoperator.seed_l = quarks[0].seed;
532  COp[ n ].baryonoperator.seed_m = quarks[1].seed;
533  COp[ n ].baryonoperator.seed_r = quarks[2].seed;
534 #endif
535  }
536  try
537  { //
538  // Creation Operator (MINUS)
539  //
540  for(int ordering=0; ordering < params.NsrcOrderings; ++ordering)
541  {
542  QDPIO::cout << "SOURCE Operator: ordering = " << ordering << std::endl;
543  swatch.start();
544  Real signs;
545  int not_used;
546  const int L=0,M=1,R=2;
547  //
548  // t = quarks[n].dilutions[i].source_header.t_source = timeslice of source
549  //
550  int qL = DilSwap( ordering, 0,1,2, L );
551  int qM = DilSwap( ordering, 0,1,2, M );
552  int qR = DilSwap( ordering, 0,1,2, R );
553  #define FIRSTTYPE
554  switch ( ordering )
555  {
556  case 0: // [012]
557  signs = 2.0; break;
558  case 1: // [210]
559  #ifdef FIRSTTYPE
560  signs = 2.0; break;
561  #else
562  signs = -2.0; break;
563  #endif
564  case 2: // [021]
565  #ifdef FIRSTTYPE
566  signs = -1.0; break;
567  #else
568  signs = 1.0; break;
569  #endif
570  case 3: // [102]
571  #ifdef FIRSTTYPE
572  signs = -1.0; break;
573  #else
574  signs = 1.0; break;
575  #endif
576  case 4: // [120]
577  signs = -1.0; break;
578  case 5: // [201]
579  signs = -1.0; break;
580  }
581  for(int qqq=0; qqq < params.NQQQs; ++qqq)
582  {
583  for(int iL=0; iL < (params.NH[ordering][L]/params.NdilReduce); ++iL)
584  for(int iM=0; iM < (params.NH[ordering][M]/params.NdilReduce); ++iM)
585  for(int iR=0; iR < (params.NH[ordering][R]/params.NdilReduce); ++iR)
586  {
587 #ifdef REDUCETOTIMEDILUTION
588  if( (quarks2[qL].dilutions[iL].source_header.t_source
589  == quarks2[qM].dilutions[iM].source_header.t_source)
590  &&(quarks2[qM].dilutions[iM].source_header.t_source
591  == quarks2[qR].dilutions[iR].source_header.t_source) )
592 #else
593  if( (quarks[qL].dilutions[iL].source_header.t_source
594  == quarks[qM].dilutions[iM].source_header.t_source)
595  &&(quarks[qM].dilutions[iM].source_header.t_source
596  == quarks[qR].dilutions[iR].source_header.t_source) )
597 #endif
598  {
599  int i0 = DilSwapInv( ordering, iL, iM, iR, 0 );
600  int i1 = DilSwapInv( ordering, iL, iM, iR, 1 );
601  int i2 = DilSwapInv( ordering, iL, iM, iR, 2 );
602  for(int t=0; t < phases.getSet().numSubsets(); ++t)
603  {
604  result[ phases.getSet()[t] ]
605 #ifdef REDUCETOTIMEDILUTION
606  = CQQQ[ qqq ](
607  quarks2[ qL ].dilutions[ iL ].source,
608  quarks2[ qM ].dilutions[ iM ].source,
609  quarks2[ qR ].dilutions[ iR ].source,
610  not_used,
611  MINUS
612  );
613 #else
614  = CQQQ[ qqq ](
615  quarks[ qL ].dilutions[ iL ].source,
616  quarks[ qM ].dilutions[ iM ].source,
617  quarks[ qR ].dilutions[ iR ].source,
618  not_used,
619  MINUS
620  );
621 #endif
622  }
623  elem = phases.sft( result );
624  for(int b=0; b < CQQQ[ qqq ].NBaryonOps; ++b)
625  {
626  for(int p=0; p < params.Nmomenta; ++p)
627  {
628  for(int t=0; t < params.nrow[ 3 ]; ++t)
629  {
630  CQQQ[ qqq ].baryon[ b ]->baryonoperator.orderings[ 0 ].op( i0, i1, i2 ).ind[ 0 ].elem( p, t )
631  += ( CQQQ[ qqq ].coef[ b ] * signs * elem( p, t ) );
632  } // t
633  } // p
634  } // Bop_index
635  } // if t_sources are all the same
636  } // iL,iM,iR
637  } // qqq
638  swatch.stop();
639  QDPIO::cout << "SOURCE operator " << ordering << " done: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
640  } // ordering
641  } // try
642  catch ( const std::string& e )
643  {
644  QDPIO::cerr << ": Caught Exception creating source baryon operator: " << e << std::endl;
645  QDP_abort( 1 );
646  }
647  QDPIO::cout << "All SOURCE operators done" << std::endl;
648 
649  // Save the operators
650  // ONLY SciDAC output format is supported!
651  swatch.start();
652  for(int b=0; b < params.Noperators; ++b)
653  {
654  XMLBufferWriter file_xml;
655  push( file_xml, params.Names[ b ] );
656  file_xml << params.param.baryon_operator;
657  write( file_xml, "Config_info", gauge_xml );
658  pop( file_xml );
659  std::string output=params.Names[ b ]+"_src";
660 
661  QDPFileWriter to( file_xml, output,
662  QDPIO_SINGLEFILE, QDPIO_SERIAL, QDPIO_OPEN );
663  // Write the scalar data
664  {
665  XMLBufferWriter record_xml;
666  write( record_xml, "SourceBaryonOperator", COp[ b ].baryonoperator );
667  write( to, record_xml, COp[ b ].baryonoperator.serialize() );
668  }
669  /*
670  // Write the scalar data
671  {
672  XMLBufferWriter record_xml;
673  write( record_xml, "SinkBaryonOperator", AOp[ b ].baryonoperator );
674  write( to, record_xml, AOp[ b ].baryonoperator.serialize() );
675  }
676  */
677  close( to );
678  }
679  swatch.stop();
680  QDPIO::cout << "Operators written: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
681 #endif // end of creation operators
682 
683  QDPIO::cout<<"end of group baryon operator calculation"<<std::endl;
684 
685 #ifdef WRITE_ALL_OUT
686  // Now writing all the operators out
687  // Save the operators
688  // ONLY SciDAC output format is supported!
689  swatch.start();
690  for(int b=0; b < params.Noperators; ++b)
691  {
692  XMLBufferWriter file_xml;
693  push( file_xml, params.Names[ b ] );
694  file_xml << params.param.baryon_operator;
695  write( file_xml, "Config_info", gauge_xml );
696  pop( file_xml );
697  QDPFileWriter to( file_xml, params.Names[ b ],
698  QDPIO_SINGLEFILE, QDPIO_SERIAL, QDPIO_OPEN );
699  // Write the scalar data
700  {
701  XMLBufferWriter record_xml;
702  write( record_xml, "SourceBaryonOperator", COp[ b ].baryonoperator );
703  write( to, record_xml, COp[ b ].baryonoperator.serialize() );
704  }
705  // Write the scalar data
706  {
707  XMLBufferWriter record_xml;
708  write( record_xml, "SinkBaryonOperator", AOp[ b ].baryonoperator );
709  write( to, record_xml, AOp[ b ].baryonoperator.serialize() );
710  }
711  close( to );
712  }
713  swatch.stop();
714  QDPIO::cout << "Operators written: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
715 #endif
716 
717  // Close the namelist output file XMLDAT
718  pop( xml_out ); // GroupBaryonAll2All
719  END_CODE();
720  // Time to bolt
722  exit( 0 );
723 }
Primary include file for CHROMA in application codes.
Primary include file for CHROMA library code.
Random complex Z(N) sources using dilution.
Class for counted reference semantics.
Definition: handle.h:33
Fourier transform phase factor support.
Definition: sftmom.h:35
multi2d< DComplex > sft(const LatticeComplex &cf) const
Do a sumMulti(cf*phases,getSet())
Definition: sftmom.cc:524
const Set & getSet() const
The set to be used in sumMulti.
Definition: sftmom.h:57
static T & Instance()
Definition: singleton.h:432
int main(int argc, char *argv[])
Construct group baryon operators.
void gaugeStartup(XMLReader &gauge_file_xml, XMLReader &gauge_xml, multi1d< LatticeColorMatrix > &u, Cfg_t &cfg)
Initialize the gauge fields.
void read(XMLReader &xml, const std::string &path, GroupBaryonOperatorEnv::Params::Qprop_t::Solutions_t &input)
Reader.
void proginfo(XMLWriter &xml)
Print out basic information about this program.
Definition: proginfo.cc:24
void zN_src(LatticeFermion &a, int N)
Volume source of Z(N) noise.
Definition: zN_src.cc:41
ForwardProp_t prop_header
Params params
unsigned j
Definition: ldumul_w.cc:35
unsigned n
Definition: ldumul_w.cc:36
void setrn(int iseed[4])
Definition: make_seeds.cc:38
int x
Definition: meslate.cc:34
int t
Definition: meslate.cc:37
const std::string name
Name to be used.
int DilSwapInv(int ord, int i, int j, int k, int which)
void ReadTextInput(Params &params, multi1d< GroupBaryonOp > &AB, multi1d< GroupBaryonOp > &CB, multi1d< GroupBaryonQQQ > &AQQQ, multi1d< GroupBaryonQQQ > &CQQQ)
Reader.
void write(XMLWriter &xml, const std::string &path, const GroupBaryonOperatorEnv::Params &param)
int DilSwap(int ord, int i, int j, int k, int which)
bool registerAll()
Register all the factories.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
LinOpSysSolverMGProtoClover::Q Q
push(xml_out,"Condensates")
Double c
Definition: invbicg.cc:108
int i
Definition: pbg5p_w.cc:55
void initialize(int *argc, char ***argv)
Chroma initialisation routine.
Definition: chroma_init.cc:114
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
void finalize(void)
Chroma finalization routine.
Definition: chroma_init.cc:308
bool linkageHack(void)
Definition: const_hmc.cc:660
std::string getXMLInputFileName()
Get input file name.
Definition: chroma_init.cc:88
pop(xml_out)
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
START_CODE()
Complex b
Definition: invbicg.cc:96
XMLFileWriter & getXMLOutputInstance()
Get xml output instance.
Definition: chroma_init.cc:359
Double zero
Definition: invbicg.cc:106
int k
Definition: invbicg.cc:119
multi1d< LatticeFermion > s(Ncb)
FloatingPoint< double > Double
Definition: gtest.h:7351
::std::string string
Definition: gtest.h:1979
Double link
Definition: pade_trln_w.cc:146
Double t_plaq
Definition: pade_trln_w.cc:145
Double w_plaq
Definition: pade_trln_w.cc:143
Double s_plaq
Definition: pade_trln_w.cc:144
All make sink constructors.
Factory for producing quark prop sinks.
All source smearing.
Factory for producing quark smearing objects.
Random complex Z(N) sources using dilution.
struct Chroma::GroupBaryonOperatorEnv::Params::Qprop_t qprop
struct Chroma::GroupBaryonOperatorEnv::Params::param_t param
struct Chroma::GroupBaryonOperatorEnv::Params::gaugestuff_t gaugestuff