CHROMA
ks_local_loops.cc
Go to the documentation of this file.
1 /* + */
2 
3 #include "fermact.h"
5 #include "meas/hadron/hadron_s.h"
6 #include "meas/sources/z2_src.h"
7 #include "meas/sources/srcfil.h"
8 
11 
12 namespace Chroma {
13 
14 
15 
16  void write_out_source_type(XMLWriter & xml_out, VolSrc_type volume_source){
17 
18 
19  if( volume_source == Z2NOISE ){
20  write(xml_out, "Random_volume_source" , "Z2NOISE");
21 
22  }else if( volume_source == GAUSSIAN ){
23  write(xml_out, "Random_volume_source" , "GAUSSIAN");
24 
25  }else if( volume_source == T_DILUTE_GAUSS ){
26  write(xml_out, "Random_volume_source" , "T_DILUTE_GAUSS");
27 
28  }else if( volume_source == C_DILUTE_GAUSS ){
29  write(xml_out, "Random_volume_source" , "C_DILUTE_GAUSS");
30 
31  }else if( volume_source == P_DILUTE_GAUSS ){
32  write(xml_out, "Random_volume_source" , "P_DILUTE_GAUSS");
33 
34  }else if( volume_source == CT_DILUTE_GAUSS ){
35  write(xml_out, "Random_volume_source" , "CT_DILUTE_GAUSS");
36 
37  }else if( volume_source == CP_DILUTE_GAUSS ){
38  write(xml_out, "Random_volume_source" , "CP_DILUTE_GAUSS");
39 
40  }else if( volume_source == PT_DILUTE_GAUSS ){
41  write(xml_out, "Random_volume_source" , "PT_DILUTE_GAUSS");
42 
43  }else if( volume_source == MOD_T_DILUTE_GAUSS ){
44  write(xml_out, "Random_volume_source" , "MOD_T_DILUTE_GAUSS");
45 
46  }else if( volume_source == CORNER_DILUTE_GAUSS ){
47  write(xml_out, "Random_volume_source" , "CORNER_DILUTE_GAUSS");
48 
49  }else if( volume_source == COR_DBL_T_DILUTE_GAUSS ){
50  write(xml_out, "Random_volume_source" , "COR_DBL_T_DILUTE_GAUSS");
51 
52  }else if( volume_source == COR_MOD_DBL_T_DILUTE_GAUSS ){
53  write(xml_out, "Random_volume_source" , "COR_MOD_DBL_T_DILUTE_GAUSS");
54 
55  }else if( volume_source == C_MOD_T_DILUTE_GAUSS ){
56  write(xml_out, "Random_volume_source" , "C_MOD_T_DILUTE_GAUSS");
57  }
58 
59 
60  }
61 
62 
63  /**************************************************************************/
64 
65  Real fill_volume_source(LatticeStaggeredFermion & q_source,
66  VolSrc_type volume_source, int t_length,
67  int *p_src_tslice, int *p_src_color_ind,
68  int *p_src_parity_ind, int *p_src_corner_ind,
69  int src_seperation, int j_decay){
70 
71  // Fill the volume with random noise
72 
73  Real coverage_fraction;
74 
75  if( volume_source == GAUSSIAN ){
76  gaussian(q_source);
77  }else if( volume_source == Z2NOISE ){
78  z2_src(q_source);
79 
80  }else if( volume_source == T_DILUTE_GAUSS ){
81  gaussian_on_timeslice(q_source,(*p_src_tslice),j_decay);
82  (*p_src_tslice)++;
83  if((*p_src_tslice)>=t_length){
84  (*p_src_tslice)=0;
85  }
86  coverage_fraction=1.0/t_length;
87 
88  }else if( volume_source == C_DILUTE_GAUSS ){
89  gaussian_color_src(q_source,(*p_src_color_ind));
90  (*p_src_color_ind)++;
91  if((*p_src_color_ind) >= Nc){
92  (*p_src_color_ind)=0;
93  }
94  coverage_fraction=1.0/Nc;
95 
96  }else if( volume_source == P_DILUTE_GAUSS ){
97  gaussian_on_parity(q_source,(*p_src_parity_ind));
98  (*p_src_parity_ind)++;
99  if((*p_src_parity_ind) > 1){
100  (*p_src_parity_ind) = 0;
101  }
102  coverage_fraction=1.0/2;
103 
104  }else if( volume_source == CT_DILUTE_GAUSS ){
105  gaussian_color_src_on_slice(q_source,(*p_src_color_ind),
106  (*p_src_tslice), j_decay);
107  (*p_src_tslice)++;
108  if((*p_src_tslice)>=t_length){
109  (*p_src_tslice)=0;
110  (*p_src_color_ind)++;
111  if((*p_src_color_ind) >=Nc){
112  (*p_src_color_ind)=0;
113  }
114  }
115  coverage_fraction=1.0/(Nc*t_length);
116 
117  }else if( volume_source == CP_DILUTE_GAUSS ){
118  gaussian_color_src_on_parity(q_source, (*p_src_color_ind),
119  (*p_src_parity_ind));
120  (*p_src_parity_ind)++;
121  if((*p_src_parity_ind) > 1){
122  (*p_src_parity_ind)=0;
123  (*p_src_color_ind)++;
124  if((*p_src_color_ind)>=Nc){
125  (*p_src_color_ind)=0;
126  }
127  }
128  coverage_fraction=1.0/(2*Nc);
129 
130  }else if( volume_source == PT_DILUTE_GAUSS ){
131  gaussian_parity_src_on_slice(q_source, (*p_src_parity_ind),
132  (*p_src_tslice), j_decay);
133  (*p_src_parity_ind)++;
134  if((*p_src_parity_ind) > 1){
135  (*p_src_parity_ind)=0;
136  (*p_src_tslice)++;
137  if((*p_src_tslice) >= t_length){
138  (*p_src_tslice) = 0;
139  }
140  }
141  coverage_fraction=1.0/(2*t_length);
142 
143  }else if( volume_source == MOD_T_DILUTE_GAUSS ){
144  gaussian_on_mod_timeslice(q_source, (*p_src_tslice), j_decay,
145  src_seperation);
146  (*p_src_tslice)++;
147  if((*p_src_tslice)>=t_length){
148  (*p_src_tslice)=0;
149  }
150  if(t_length%src_seperation==0){
151  coverage_fraction=(1.0/(src_seperation));
152  }else{
153  coverage_fraction=9999999;
154  }
155 
156  }else if( volume_source == CORNER_DILUTE_GAUSS ){
157  gaussian_on_corner(q_source, (*p_src_corner_ind));
158  (*p_src_corner_ind)++;
159  if((*p_src_corner_ind)>=16){
160  (*p_src_corner_ind)=0;
161  }
162  coverage_fraction=1.0/16;
163 
164  }else if( volume_source == COR_DBL_T_DILUTE_GAUSS ){
165  gaussian_corner_on_dbl_slice(q_source, (*p_src_corner_ind),
166  (*p_src_tslice), j_decay);
167  (*p_src_corner_ind)++;
168  if((*p_src_corner_ind) >= 16){
169  (*p_src_corner_ind) = 0;
170  (*p_src_tslice)++;
171  if((*p_src_tslice)>=t_length){
172  (*p_src_tslice) = 0;
173  }
174  }
175  coverage_fraction=1.0/(16*t_length);
176 
177  }else if( volume_source == COR_MOD_DBL_T_DILUTE_GAUSS ){
178  gaussian_corner_on_mod_dbl_slice(q_source,(*p_src_corner_ind),
179  (*p_src_tslice), j_decay,
180  src_seperation);
181  (*p_src_corner_ind)++;
182  if((*p_src_corner_ind) >= 16){
183  (*p_src_corner_ind)=0;
184  (*p_src_tslice)++;
185  if((*p_src_tslice)>=t_length){
186  (*p_src_tslice)=0;
187  }
188  }
189  if(t_length%src_seperation==0){
190  coverage_fraction=1.0/(2*Nc);
191  }else{
192  coverage_fraction=9999999;
193  }
194 
195  }else if( volume_source == C_MOD_T_DILUTE_GAUSS ){
196  gaussian_color_src_on_mod_slice(q_source,(*p_src_color_ind),
197  (*p_src_tslice), j_decay, src_seperation);
198  (*p_src_tslice)++;
199  if((*p_src_tslice)>=t_length){
200  (*p_src_tslice)=0;
201  (*p_src_color_ind)++;
202  if((*p_src_color_ind)>=Nc){
203  (*p_src_color_ind) = 0;
204  }
205  }
206 
207  coverage_fraction=1.0/(Nc*t_length);
208 
209  }else{
210  QDP_error_exit("Wrong type of volume source");
211  }
212 
213  /**********************************************/
214  /* local src for testing */
215 
216  // q_source=zero;
217  // multi1d<int> coord(Nd);
218  // coord[0]=0; coord[1] = 0; coord[2] = 0; coord[3] = 0;
219  // srcfil(q_source, coord,0 ) ;
220  // srcfil(q_source, coord,1 ) ;
221  // srcfil(q_source, coord,2 ) ;
222 
223 
224 
225 
226 
227  /***********************************************/
228 
229 
230 
231  return coverage_fraction;
232 
233 
234  }
235 
236 
237  /**************************************************************************/
238 
239 
240 
243  LatticeStaggeredFermion & q_source,
244  LatticeStaggeredFermion & psi ,
245  const multi1d<LatticeColorMatrix> & u,
246  XMLWriter & xml_out,
247  bool gauge_shift,
248  bool sym_shift,
249  bool loop_checkpoint,
250  int t_length,
251  Real Mass,
252  int Nsamp,
253  Real RsdCG,
254  int CFGNO,
255  VolSrc_type volume_source,
256  int src_seperation,
257  int j_decay){
258 
259 
260  push(xml_out,"local_loops_s");
261 
262  // write common parameters
263  write(xml_out, "Mass" , Mass);
264 
265  write_out_source_type(xml_out, volume_source);
266 
267  write(xml_out, "Number_of_samples" , Nsamp);
268 
269 
270  Stag_shift_option type_of_shift ;
271  if( gauge_shift ){
272  if(sym_shift){
273  type_of_shift = SYM_GAUGE_INVAR ;
274  }else{
275  type_of_shift = GAUGE_INVAR ;
276  }
277  }else{
278  if(sym_shift){
279  type_of_shift = SYM_NON_GAUGE_INVAR ;
280  }else{
281  type_of_shift = NON_GAUGE_INVAR ;
282  }
283  }
284 
285 
286  // set up the loop code
287  local_scalar_loop scalar_one_loop(t_length,Nsamp,
288  u,type_of_shift) ;
289  local_scalar_kilcup_loop scalar_one_kilcup_loop(t_length,
290  Nsamp, u,
291  type_of_shift);
292  non_local_scalar_loop scalar_two_loop(t_length,Nsamp,
293  u,type_of_shift) ;
294  fourlink_scalar_loop scalar_four_loop(t_length, Nsamp,
295  u, type_of_shift);
296  fourlink_scalar_kilcup_loop scalar_four_kilcup_loop(t_length, Nsamp,
297  u, type_of_shift);
298  threelink_pseudoscalar_loop eta3_loop(t_length,Nsamp,
299  u,type_of_shift) ;
300  fourlink_pseudoscalar_loop eta4_loop(t_length,Nsamp,
301  u,type_of_shift) ;
302 
303  fourlink_pseudoscalar_kilcup_loop eta4_kilcup_loop(t_length,Nsamp,
304  u,type_of_shift) ;
305 
306  zerolink_pseudoscalar_loop eta0_loop(t_length, Nsamp,
307  u,type_of_shift) ;
308 
309 
310  // Seed the RNG with the cfg number for now
311  QDP::Seed seed;
312  seed = CFGNO;
313  RNG::setrn(seed);
314 
315  int src_tslice=0;
316  int src_color_ind = 0;
317  int src_parity_ind = 0;
318  int src_corner_ind =0;
319 
320  Real coverage_fraction;
321 
322  for(int i = 0; i < Nsamp; ++i){
323  psi = zero; // note this is ``zero'' and not 0
324  RNG::savern(seed);
325  QDPIO::cout << "SEED = " << seed << std::endl;
326 
327  QDPIO::cout << "Noise sample: " << i << std::endl;
328 
329  // Fill the volume with random noise
330  coverage_fraction = fill_volume_source(q_source, volume_source,
331  t_length, &src_tslice,
332  &src_color_ind, &src_parity_ind,
333  &src_corner_ind, src_seperation,
334  j_decay);
335 
336 
337  // Compute the solution std::vector for the particular source
338  SystemSolverResults_t res = (*qprop)(psi, q_source);
339 
340  push(xml_out,"Qprop_noise");
341  write(xml_out, "Noise_number" , i);
342  write(xml_out, "RsdCG" , RsdCG);
343  write(xml_out, "n_count", res.n_count);
344  write(xml_out, "Seed" , seed);
345  pop(xml_out);
346 
347 
348  scalar_one_loop.compute(q_source,psi,i) ;
349  scalar_one_kilcup_loop.compute(psi,i,Mass);
350  scalar_two_loop.compute(q_source,psi,i) ;
351  scalar_four_loop.compute(q_source,psi,i);
352  scalar_four_kilcup_loop.compute(psi,i,Mass);
353  eta3_loop.compute(q_source,psi,i) ;
354  eta4_loop.compute(q_source,psi,i) ;
355  eta4_kilcup_loop.compute(psi,i, Mass);
356  eta0_loop.compute(q_source,psi,i) ;
357 
358  if(loop_checkpoint){
359  //write each measurement to the XML file
360 
361  scalar_one_loop.dump(xml_out,i) ;
362  scalar_one_kilcup_loop.dump(xml_out,i);
363  scalar_two_loop.dump(xml_out,i) ;
364  scalar_four_loop.dump(xml_out,i) ;
365  scalar_four_kilcup_loop.dump(xml_out,i) ;
366  eta3_loop.dump(xml_out,i) ;
367  eta4_loop.dump(xml_out,i) ;
368  eta4_kilcup_loop.dump(xml_out,i) ;
369  eta0_loop.dump(xml_out,i) ;
370  }
371 
372  } // Nsamples
373 
374 
375  // write output from the loop calc
376  scalar_one_loop.dump(xml_out) ;
377  scalar_one_kilcup_loop.dump(xml_out);
378  scalar_two_loop.dump(xml_out) ;
379  scalar_four_loop.dump(xml_out) ;
380  scalar_four_kilcup_loop.dump(xml_out) ;
381  eta3_loop.dump(xml_out) ;
382  eta4_loop.dump(xml_out) ;
383  eta4_kilcup_loop.dump(xml_out) ;
384  eta0_loop.dump(xml_out) ;
385 
386  // end of this section
387  pop(xml_out);
388 
389  }
390 
391  /**********************************************************************/
392 
393 
394  //
395  // version used in test code.
396  //
397 
400  LatticeStaggeredFermion & q_source,
401  LatticeStaggeredFermion & psi ,
402  multi1d<LatticeColorMatrix> & u,
403  XMLFileWriter & xml_out,
404  XMLReader & xml_in ,
405  int t_length,
406  Real Mass,
407  int Nsamp,
408  Real RsdCG,
409  int CFGNO,
410  VolSrc_type volume_source,
411  int src_seperation,
412  int j_decay){
413 
414  int src_tslice=0;
415  int src_color_ind = 0;
416  int src_parity_ind = 0;
417  int src_corner_ind =0;
418  Real coverage_fraction;
419 
420  push(xml_out,"ks_local_loops");
421 
422  // write common parameters
423  write(xml_out, "Mass" , Mass);
424 
425  write_out_source_type(xml_out, volume_source);
426 
427  write(xml_out, "Number_of_samples" , Nsamp);
428 
429 
430  // parse input files
431 
432 
433  // the wrapped disconnected loops
434  bool gauge_shift ;
435  bool sym_shift ;
436  bool loop_checkpoint;
437 
438  try{
439  read(xml_in, "/propagator/param/use_gauge_invar_oper", gauge_shift ) ;
440  }catch (const std::string& e){
441  QDPIO::cerr << "Error reading data: " << e << std::endl;
442  throw;
443  }
444  try{
445  read(xml_in, "/propagator/param/use_sym_shift_oper", sym_shift ) ;
446  }catch (const std::string& e){
447  QDPIO::cerr << "Error reading data: " << e << std::endl;
448  throw;
449  }
450  try{
451  read(xml_in, "/propagator/param/loop_checkpoint", loop_checkpoint ) ;
452  }catch (const std::string& e){
453  QDPIO::cerr << "Error reading data: " << e << std::endl;
454  throw;
455  }
456 
457 
458  Stag_shift_option type_of_shift ;
459  if( gauge_shift ){
460  if(sym_shift){
461  type_of_shift = SYM_GAUGE_INVAR ;
462  }else{
463  type_of_shift = GAUGE_INVAR ;
464  }
465  }else{
466  if(sym_shift){
467  type_of_shift = SYM_NON_GAUGE_INVAR ;
468  }else{
469  type_of_shift = NON_GAUGE_INVAR ;
470  }
471  }
472 
473  // set up the loop code
474  local_scalar_loop scalar_one_loop(t_length,Nsamp,
475  u,type_of_shift) ;
476  local_scalar_kilcup_loop scalar_one_kilcup_loop(t_length, Nsamp,
477  u, type_of_shift);
478  non_local_scalar_loop scalar_two_loop(t_length,Nsamp,
479  u,type_of_shift) ;
480  fourlink_scalar_loop scalar_four_loop(t_length, Nsamp,
481  u, type_of_shift);
482  fourlink_scalar_kilcup_loop scalar_four_kilcup_loop(t_length, Nsamp,
483  u, type_of_shift);
484  threelink_pseudoscalar_loop eta3_loop(t_length,Nsamp,
485  u,type_of_shift) ;
486  fourlink_pseudoscalar_loop eta4_loop(t_length,Nsamp,
487  u,type_of_shift) ;
488  fourlink_pseudoscalar_kilcup_loop eta4_kilcup_loop(t_length, Nsamp,
489  u,type_of_shift) ;
490 
491  // for test purposes
492  zerolink_pseudoscalar_loop eta0_loop(t_length, Nsamp,
493  u,type_of_shift) ;
494 
495  // Seed the RNG with the cfg number for now
496  QDP::Seed seed;
497  seed = CFGNO;
498  RNG::setrn(seed);
499 
500  for(int i = 0; i < Nsamp; ++i){
501  psi = zero; // note this is ``zero'' and not 0
502  RNG::savern(seed);
503  QDPIO::cout << "SEED = " << seed << std::endl;
504 
505  QDPIO::cout << "Noise sample: " << i << std::endl;
506 
507  // Fill the volume with random noise
508 
509  coverage_fraction = fill_volume_source(q_source, volume_source,
510  t_length, &src_tslice,
511  &src_color_ind, &src_parity_ind,
512  &src_corner_ind, src_seperation,
513  j_decay);
514 
515  // Compute the solution std::vector for the particular source
516  SystemSolverResults_t res = (*qprop)(psi, q_source);
517 
518  push(xml_out,"Qprop_noise");
519  write(xml_out, "Noise_number" , i);
520  write(xml_out, "RsdCG" , RsdCG);
521  write(xml_out, "n_count", res.n_count);
522  write(xml_out, "Seed" , seed);
523  pop(xml_out);
524 
525 
526  scalar_one_loop.compute(q_source,psi,i) ;
527  scalar_one_kilcup_loop.compute(psi, i, Mass);
528  scalar_two_loop.compute(q_source,psi,i) ;
529  scalar_four_loop.compute(q_source,psi,i);
530  scalar_four_kilcup_loop.compute(psi,i,Mass);
531  eta3_loop.compute(q_source,psi,i) ;
532  eta4_loop.compute(q_source,psi,i) ;
533  eta4_kilcup_loop.compute(psi,i, Mass);
534  eta0_loop.compute(q_source,psi,i) ;
535 
536  if(loop_checkpoint){
537  //write each measurement to the XML file
538 
539  scalar_one_loop.dump(xml_out,i) ;
540  scalar_one_kilcup_loop.dump(xml_out,i);
541  scalar_two_loop.dump(xml_out,i) ;
542  scalar_four_loop.dump(xml_out,i) ;
543  scalar_four_kilcup_loop.dump(xml_out,i) ;
544  eta3_loop.dump(xml_out,i) ;
545  eta4_loop.dump(xml_out,i) ;
546  eta4_kilcup_loop.dump(xml_out,i) ;
547  eta0_loop.dump(xml_out,i) ;
548  }
549 
550  } // Nsamples
551 
552 
553  // write output from the loop calc
554  scalar_one_loop.dump(xml_out) ;
555  scalar_one_kilcup_loop.dump(xml_out);
556  scalar_two_loop.dump(xml_out) ;
557  scalar_four_loop.dump(xml_out) ;
558  scalar_four_kilcup_loop.dump(xml_out) ;
559  eta3_loop.dump(xml_out) ;
560  eta4_loop.dump(xml_out) ;
561  eta4_kilcup_loop.dump(xml_out) ;
562  eta0_loop.dump(xml_out) ;
563 
564  // end of this section
565  pop(xml_out);
566 
567  }
568 
569  /**********************************************************************/
570 
571  // fuzz the loops
572 
575  LatticeStaggeredFermion & q_source,
576  LatticeStaggeredFermion & psi ,
577  LatticeStaggeredFermion & psi_fuzz ,
578  const multi1d<LatticeColorMatrix> & u,
579  const multi1d<LatticeColorMatrix> & u_smr,
580  XMLWriter & xml_out,
581  bool gauge_shift,
582  bool sym_shift,
583  bool loop_checkpoint,
584  int t_length,
585  Real Mass,
586  int Nsamp,
587  Real RsdCG,
588  int CFGNO,
589  VolSrc_type volume_source,
590  int fuzz_width,
591  int src_seperation,
592  int j_decay){
593 
594  int src_tslice=0;
595  int src_color_ind = 0;
596  int src_parity_ind = 0;
597  int src_corner_ind =0;
598  Real coverage_fraction;
599 
600  push(xml_out,"fuzz_loops_s");
601 
602  // write common parameters
603  write(xml_out, "Mass" , Mass);
604 
605 
606  write_out_source_type(xml_out, volume_source);
607 
608 
609  write(xml_out, "Number_of_samples" , Nsamp);
610  write(xml_out, "fuzz_width" , fuzz_width);
611 
612 
613  Stag_shift_option type_of_shift ;
614  if( gauge_shift ){
615  if(sym_shift){
616  type_of_shift = SYM_GAUGE_INVAR ;
617  }else{
618  type_of_shift = GAUGE_INVAR ;
619  }
620  }else{
621  if(sym_shift){
622  type_of_shift = SYM_NON_GAUGE_INVAR ;
623  }else{
624  type_of_shift = NON_GAUGE_INVAR ;
625  }
626  }
627 
628 
629 
630  // set up the loop code
631  local_scalar_loop scalar_one_loop(t_length,Nsamp,
632  u,type_of_shift) ;
633  local_scalar_kilcup_loop scalar_one_kilcup_loop(t_length,Nsamp,
634  u, type_of_shift);
635  non_local_scalar_loop scalar_two_loop(t_length,Nsamp,
636  u,type_of_shift) ;
637  fourlink_scalar_loop scalar_four_loop(t_length, Nsamp,
638  u, type_of_shift);
639  fourlink_scalar_kilcup_loop scalar_four_kilcup_loop(t_length, Nsamp,
640  u, type_of_shift);
641  threelink_pseudoscalar_loop eta3_loop(t_length,Nsamp,
642  u,type_of_shift) ;
643  fourlink_pseudoscalar_loop eta4_loop(t_length,Nsamp,
644  u,type_of_shift) ;
645 
646  fourlink_pseudoscalar_kilcup_loop eta4_kilcup_loop(t_length, Nsamp,
647  u,type_of_shift) ;
648 
649  // for test purposes
650  zerolink_pseudoscalar_loop eta0_loop(t_length, Nsamp,
651  u,type_of_shift) ;
652 
653  // Seed the RNG with the cfg number for now
654  QDP::Seed seed;
655  seed = CFGNO;
656  RNG::setrn(seed);
657 
658 
659  for(int i = 0; i < Nsamp; ++i){
660  psi = zero; // note this is ``zero'' and not 0
661  RNG::savern(seed);
662  QDPIO::cout << "SEED = " << seed << std::endl;
663 
664  QDPIO::cout << "Noise sample: " << i << std::endl;
665 
666  // Fill the volume with random noise
667  coverage_fraction = fill_volume_source(q_source, volume_source,
668  t_length, &src_tslice,
669  &src_color_ind, &src_parity_ind,
670  &src_corner_ind, src_seperation,
671  j_decay);
672 
673  // Compute the solution std::vector for the particular source
674  SystemSolverResults_t res = (*qprop)(psi, q_source);
675 
676  push(xml_out,"Qprop_noise");
677  write(xml_out, "Noise_number" , i);
678  write(xml_out, "RsdCG" , RsdCG);
679  write(xml_out, "n_count", res.n_count);
680  write(xml_out, "Seed" , seed);
681  pop(xml_out);
682 
683 
684  fuzz_smear(u_smr, psi,psi_fuzz,
685  fuzz_width, j_decay) ;
686 
687 
688  scalar_one_loop.compute(q_source,psi_fuzz,i) ;
689  scalar_one_kilcup_loop.compute(psi, i, Mass);
690  scalar_two_loop.compute(q_source,psi_fuzz,i) ;
691  scalar_four_loop.compute(q_source,psi,i);
692  scalar_four_kilcup_loop.compute(psi,i,Mass);
693  eta3_loop.compute(q_source,psi_fuzz,i) ;
694  eta4_loop.compute(q_source,psi_fuzz,i) ;
695  eta4_kilcup_loop.compute(psi,i, Mass);
696  eta0_loop.compute(q_source,psi,i) ;
697 
698 
699  if(loop_checkpoint){
700  //write each measurement to the XML file
701 
702  scalar_one_loop.dump(xml_out,i) ;
703  scalar_one_kilcup_loop.dump(xml_out,i);
704  scalar_two_loop.dump(xml_out,i) ;
705  scalar_four_loop.dump(xml_out,i) ;
706  scalar_four_kilcup_loop.dump(xml_out,i) ;
707  eta3_loop.dump(xml_out,i) ;
708  eta4_loop.dump(xml_out,i) ;
709  eta4_kilcup_loop.dump(xml_out,i) ;
710  eta0_loop.dump(xml_out,i) ;
711  }
712 
713  } // Nsamples
714 
715 
716  // write output from the loop calc
717  scalar_one_loop.dump(xml_out) ;
718  scalar_one_kilcup_loop.dump(xml_out);
719  scalar_two_loop.dump(xml_out) ;
720  scalar_four_loop.dump(xml_out) ;
721  scalar_four_kilcup_loop.dump(xml_out) ;
722  eta3_loop.dump(xml_out) ;
723  eta4_loop.dump(xml_out) ;
724  eta4_kilcup_loop.dump(xml_out) ;
725  eta0_loop.dump(xml_out) ;
726 
727  // end of this section
728  pop(xml_out);
729 
730  }
731 
732  /**********************************************************************/
733 
734  // fuzz the loops
735  // HACK
736 
739  LatticeStaggeredFermion & q_source,
740  LatticeStaggeredFermion & psi ,
741  LatticeStaggeredFermion & psi_fuzz ,
742  const multi1d<LatticeColorMatrix> & u,
743  const multi1d<LatticeColorMatrix> & u_smr,
744  XMLWriter & xml_out,
745  bool gauge_shift,
746  bool sym_shift,
747  bool loop_checkpoint,
748  int t_length,
749  Real Mass,
750  int Nsamp,
751  Real RsdCG,
752  int CFGNO,
753  VolSrc_type volume_source,
754  int fuzz_width,
755  int src_seperation,
756  int j_decay, bool binary_loop_checkpoint,
757  std::string binary_name){
758 
759  int src_tslice=0;
760  int src_color_ind = 0;
761  int src_parity_ind = 0;
762  int src_corner_ind =0;
763  Real coverage_fraction;
764  int j;
765  int fuzz_index;
766 
767  push(xml_out,"fuzz_loops_s");
768 
769  // write common parameters
770  write(xml_out, "Mass" , Mass);
771 
772 
773  write_out_source_type(xml_out, volume_source);
774 
775 
776  write(xml_out, "Number_of_samples" , Nsamp);
777  write(xml_out, "fuzz_width" , fuzz_width);
778 
779 
780  Stag_shift_option type_of_shift ;
781  if( gauge_shift ){
782  if(sym_shift){
783  type_of_shift = SYM_GAUGE_INVAR ;
784  }else{
785  type_of_shift = GAUGE_INVAR ;
786  }
787  }else{
788  if(sym_shift){
789  type_of_shift = SYM_NON_GAUGE_INVAR ;
790  }else{
791  type_of_shift = NON_GAUGE_INVAR ;
792  }
793  }
794 
795  // set up the loop code
796  local_scalar_loop scalar_one_loop(t_length,Nsamp,
797  u,type_of_shift) ;
798  local_scalar_kilcup_loop scalar_one_kilcup_loop(t_length, Nsamp,
799  u, type_of_shift);
800  non_local_scalar_loop scalar_two_loop(t_length,Nsamp,
801  u,type_of_shift) ;
802  fourlink_scalar_loop scalar_four_loop(t_length, Nsamp,
803  u, type_of_shift);
804  fourlink_scalar_kilcup_loop scalar_four_kilcup_loop(t_length, Nsamp,
805  u, type_of_shift);
806  threelink_pseudoscalar_loop eta3_loop(t_length,Nsamp,
807  u,type_of_shift) ;
808  fourlink_pseudoscalar_loop eta4_loop(t_length,Nsamp,
809  u,type_of_shift) ;
810 
811  fourlink_pseudoscalar_kilcup_loop eta4_kilcup_loop(t_length, Nsamp,
812  u,type_of_shift) ;
813 
814  // for test purposes
815  zerolink_pseudoscalar_loop eta0_loop(t_length, Nsamp,
816  u,type_of_shift);
817 
818 
819 
820  // fuzzed loops here
821  local_scalar_loop_fuzz scalar_one_loop_fuzz(t_length, Nsamp, u,
822  type_of_shift);
823  local_scalar_kilcup_loop_fuzz scalar_one_kilcup_loop_fuzz(t_length, Nsamp,
824  u, type_of_shift);
825  non_local_scalar_loop_fuzz scalar_two_loop_fuzz(t_length,Nsamp, u,
826  type_of_shift);
827 
828  threelink_pseudoscalar_loop_fuzz eta3_loop_fuzz(t_length,Nsamp,
829  u,type_of_shift);
830  fourlink_pseudoscalar_loop_fuzz eta4_loop_fuzz(t_length,Nsamp,
831  u,type_of_shift);
832 
833  fourlink_pseudoscalar_kilcup_loop_fuzz eta4_kilcup_loop_fuzz(t_length,
834  Nsamp, u,
835  type_of_shift);
836 
837 
838 
839  // set up the loop code
840 
841  // Seed the RNG with the cfg number for now
842  QDP::Seed seed;
843  seed = CFGNO;
844  RNG::setrn(seed);
845 
846 
847  for(int i = 0; i < Nsamp; ++i){
848  psi = zero; // note this is ``zero'' and not 0
849  RNG::savern(seed);
850  QDPIO::cout << "SEED = " << seed << std::endl;
851 
852  QDPIO::cout << "Noise sample: " << i << std::endl;
853 
854  // Fill the volume with random noise
855  coverage_fraction = fill_volume_source(q_source, volume_source,
856  t_length, &src_tslice,
857  &src_color_ind, &src_parity_ind,
858  &src_corner_ind, src_seperation,
859  j_decay);
861  {
862  StopWatch swatch;
863  swatch.start();
864 
865  // Compute the solution std::vector for the particular source
866  res = (*qprop)(psi, q_source);
867 
868  swatch.stop();
869  double time_in_sec = swatch.getTimeInSeconds();
870  QDPIO::cout << "ks_fuzz_loops::PROF INVERTER [" << i << "] " << time_in_sec << " sec" << std::endl;
871 
872 
873 
874  }
875 
876 
877  push(xml_out,"Qprop_noise");
878  write(xml_out, "Noise_number" , i);
879  write(xml_out, "RsdCG" , RsdCG);
880  write(xml_out, "n_count", res.n_count);
881  write(xml_out, "Seed" , seed);
882  pop(xml_out);
883 
884  {
885  StopWatch swatch;
886  swatch.start();
887 
888  fuzz_smear(u_smr, psi,psi_fuzz, fuzz_width, j_decay) ;
889 
890  swatch.stop();
891  double time_in_sec = swatch.getTimeInSeconds();
892  QDPIO::cout << "ks_fuzz_loops::PROF fuzz_smear [" << i << "] " << time_in_sec << " sec" << std::endl;
893 
894 
895  }
896 
897  {
898  StopWatch swatch;
899  swatch.start();
900 
901 
902  // compute the un-fuzzed operators
903  scalar_one_loop.compute(q_source,psi,i);
904  scalar_one_kilcup_loop.compute(psi, i, Mass);
905  scalar_two_loop.compute(q_source,psi,i);
906  scalar_four_loop.compute(q_source,psi,i);
907  scalar_four_kilcup_loop.compute(psi,i,Mass);
908  eta3_loop.compute(q_source,psi,i);
909  eta4_loop.compute(q_source,psi,i);
910  eta4_kilcup_loop.compute(psi,i, Mass);
911  eta0_loop.compute(q_source,psi,i);
912 
913  // now compute the fuzzed operators
914  scalar_one_loop_fuzz.compute(q_source,psi_fuzz,i);
915  scalar_one_kilcup_loop_fuzz.compute(psi_fuzz, psi, i, Mass);
916  scalar_two_loop_fuzz.compute(q_source,psi_fuzz,i) ;
917  eta3_loop_fuzz.compute(q_source,psi_fuzz,i) ;
918  eta4_loop_fuzz.compute(q_source,psi_fuzz,i) ;
919  eta4_kilcup_loop_fuzz.compute(psi_fuzz,psi,i, Mass);
920 
921 
922  swatch.stop();
923  double time_in_sec = swatch.getTimeInSeconds();
924  QDPIO::cout << "ks_fuzz_loops::PROF compute [" << i << "] " << time_in_sec << " sec" << std::endl;
925 
926 
927  }
928 
929 
930  if(loop_checkpoint){
931  StopWatch swatch;
932  swatch.start();
933 
934  //write each measurement to the XML file
935 
936  scalar_one_loop.dump(xml_out,i) ;
937  scalar_one_kilcup_loop.dump(xml_out,i);
938  scalar_two_loop.dump(xml_out,i) ;
939  scalar_four_loop.dump(xml_out,i) ;
940  scalar_four_kilcup_loop.dump(xml_out,i) ;
941  eta3_loop.dump(xml_out,i) ;
942  eta4_loop.dump(xml_out,i) ;
943  eta4_kilcup_loop.dump(xml_out,i) ;
944  eta0_loop.dump(xml_out,i) ;
945 
946  scalar_one_loop_fuzz.dump(xml_out,i) ;
947  scalar_one_kilcup_loop_fuzz.dump(xml_out,i);
948  scalar_two_loop_fuzz.dump(xml_out,i) ;
949  eta3_loop_fuzz.dump(xml_out,i) ;
950  eta4_loop_fuzz.dump(xml_out,i) ;
951  eta4_kilcup_loop_fuzz.dump(xml_out,i) ;
952 
953  swatch.stop();
954  double time_in_sec = swatch.getTimeInSeconds();
955  QDPIO::cout << "ks_fuzz_loops::PROF CHECKPOINT [" << i << "] " << time_in_sec << " sec" << std::endl;
956 
957 
958  }
959 
960 
961  } // Nsamples
962 
963  // write output from the loop calc
964  {
965  StopWatch swatch;
966  swatch.start();
967 
968  scalar_one_loop.dump(xml_out) ;
969  scalar_one_kilcup_loop.dump(xml_out);
970  scalar_two_loop.dump(xml_out) ;
971  scalar_four_loop.dump(xml_out) ;
972  scalar_four_kilcup_loop.dump(xml_out) ;
973  eta3_loop.dump(xml_out) ;
974  eta4_loop.dump(xml_out) ;
975  eta4_kilcup_loop.dump(xml_out) ;
976 
977  eta0_loop.dump(xml_out) ;
978 
979  scalar_one_loop_fuzz.dump(xml_out);
980  scalar_one_kilcup_loop_fuzz.dump(xml_out);
981  scalar_two_loop_fuzz.dump(xml_out) ;
982  eta3_loop_fuzz.dump(xml_out) ;
983  eta4_loop_fuzz.dump(xml_out) ;
984  eta4_kilcup_loop_fuzz.dump(xml_out) ;
985 
986 
987  swatch.stop();
988  double time_in_sec = swatch.getTimeInSeconds();
989  QDPIO::cout << "ks_fuzz_loops::FINAL IO " << time_in_sec << " sec" << std::endl;
990 
991  }
992 
993  if(binary_loop_checkpoint )
994  {
995  // BINARY DUMP
996  StopWatch swatch;
997  swatch.start();
998 
999  scalar_one_loop.binary_dump(binary_name) ;
1000  scalar_one_kilcup_loop.binary_dump(binary_name) ;
1001  scalar_two_loop.binary_dump(binary_name) ;
1002  scalar_four_loop.binary_dump(binary_name);
1003  scalar_four_kilcup_loop.binary_dump(binary_name);
1004  eta3_loop.binary_dump(binary_name) ;
1005  eta4_loop.binary_dump(binary_name) ;
1006  eta4_kilcup_loop.binary_dump(binary_name) ;
1007  eta0_loop.binary_dump(binary_name) ;
1008 
1009  scalar_one_loop_fuzz.binary_dump(binary_name) ;
1010  scalar_one_kilcup_loop_fuzz.binary_dump(binary_name) ;
1011  scalar_two_loop_fuzz.binary_dump(binary_name) ;
1012  eta3_loop_fuzz.binary_dump(binary_name) ;
1013  eta4_loop_fuzz.binary_dump(binary_name) ;
1014  eta4_kilcup_loop_fuzz.binary_dump(binary_name) ;
1015 
1016  swatch.stop();
1017  double time_in_sec = swatch.getTimeInSeconds();
1018  QDPIO::cout << "ks_fuzz_loops::BINARY IO " << time_in_sec << " sec" << std::endl;
1019 
1020  }
1021 
1022 
1023  // end of this section
1024  pop(xml_out);
1025 
1026  }
1027 
1028 
1029  /**********************************************************************/
1030 
1031 
1034  LatticeStaggeredFermion & q_source1,
1035  LatticeStaggeredFermion & psi1 ,
1036  const multi1d<LatticeColorMatrix> & u,
1037  XMLWriter & xml_out,
1038  bool gauge_shift,
1039  bool sym_shift,
1040  bool loop_checkpoint,
1041  int t_length,
1042  Real Mass,
1043  int Nsamp,
1044  Real RsdCG,
1045  int CFGNO,
1046  VolSrc_type volume_source,
1047  int src_seperation,
1048  int j_decay){
1049 
1050 
1051  LatticeStaggeredFermion psi2 ;
1052  LatticeStaggeredFermion q_source2 ;
1053 
1054 
1055  push(xml_out,"local_loops_s");
1056 
1057  // write common parameters
1058  write(xml_out, "Mass" , Mass);
1059 
1060  write_out_source_type(xml_out, volume_source);
1061 
1062  write(xml_out, "Number_of_samples" , Nsamp);
1063 
1064  Stag_shift_option type_of_shift ;
1065  if( gauge_shift ){
1066  if(sym_shift){
1067  type_of_shift = SYM_GAUGE_INVAR ;
1068  }else{
1069  type_of_shift = GAUGE_INVAR ;
1070  }
1071  }else{
1072  if(sym_shift){
1073  type_of_shift = SYM_NON_GAUGE_INVAR ;
1074  }else{
1075  type_of_shift = NON_GAUGE_INVAR ;
1076  }
1077  }
1078 
1079  // set up the loop code
1080  local_scalar_loop scalar_one_loop(t_length,Nsamp,
1081  u,type_of_shift) ;
1082  local_scalar_kilcup_loop scalar_one_kilcup_loop(t_length, Nsamp,
1083  u, type_of_shift);
1084  non_local_scalar_loop scalar_two_loop(t_length,Nsamp,
1085  u,type_of_shift) ;
1086  fourlink_scalar_loop scalar_four_loop(t_length, Nsamp,
1087  u, type_of_shift);
1088  fourlink_scalar_kilcup_loop scalar_four_kilcup_loop(t_length, Nsamp,
1089  u, type_of_shift);
1090  threelink_pseudoscalar_loop eta3_loop(t_length,Nsamp,
1091  u,type_of_shift) ;
1092  fourlink_pseudoscalar_loop eta4_loop(t_length,Nsamp,
1093  u,type_of_shift) ;
1094 
1095  fourlink_pseudoscalar_kilcup_loop eta4_kilcup_loop(t_length,Nsamp,
1096  u,type_of_shift) ;
1097 
1098  zerolink_pseudoscalar_loop eta0_loop(t_length, Nsamp,
1099  u,type_of_shift) ;
1100 
1101  fourlink_pseudoscalar_stoch_conn eta4_conn(t_length,Nsamp,
1102  u,type_of_shift) ;
1103 
1104  // Seed the RNG with the cfg number for now
1105  QDP::Seed seed;
1106  seed = CFGNO;
1107  RNG::setrn(seed);
1108 
1109  int src_tslice=0;
1110  int src_color_ind = 0;
1111  int src_parity_ind = 0;
1112  int src_corner_ind =0;
1113 
1114  int src_tslice2=0;
1115  int src_color_ind2 = 0;
1116  int src_parity_ind2 = 0;
1117  int src_corner_ind2 =0;
1118 
1119  Real coverage_fraction;
1120 
1121  for(int i = 0; i < Nsamp; ++i){
1122  psi1 = zero; // note this is ``zero'' and not 0
1123  psi2 = zero; // note this is ``zero'' and not 0
1124  RNG::savern(seed);
1125  QDPIO::cout << "SEED = " << seed << std::endl;
1126 
1127  QDPIO::cout << "Noise sample: " << i << std::endl;
1128 
1129  // Fill the volume with random noise
1130  coverage_fraction = fill_volume_source(q_source1, volume_source,
1131  t_length, &src_tslice,
1132  &src_color_ind, &src_parity_ind,
1133  &src_corner_ind, src_seperation,
1134  j_decay);
1135 
1136  //fill 2nd source for stochastic connected correlators
1137  coverage_fraction = fill_volume_source(q_source2, volume_source,
1138  t_length, &src_tslice2,
1139  &src_color_ind2, &src_parity_ind2,
1140  &src_corner_ind2, src_seperation,
1141  j_decay);
1142 
1143  // Compute the solution std::vector for the particular source
1144  SystemSolverResults_t res1 = (*qprop)(psi1, q_source1);
1145  SystemSolverResults_t res2 = (*qprop)(psi2, q_source2);
1146 
1147  push(xml_out,"Qprop_noise");
1148  write(xml_out, "Noise_number" , i);
1149  write(xml_out, "RsdCG" , RsdCG);
1150  write(xml_out, "n_count1", res1.n_count);
1151  write(xml_out, "n_count2", res2.n_count);
1152  write(xml_out, "Seed" , seed);
1153  pop(xml_out);
1154 
1155 
1156  scalar_one_loop.compute(q_source1,psi1,i) ;
1157  scalar_one_kilcup_loop.compute(psi1,i, Mass);
1158  scalar_two_loop.compute(q_source1,psi1,i) ;
1159  scalar_four_loop.compute(q_source1,psi1,i);
1160  scalar_four_kilcup_loop.compute(psi1,i,Mass);
1161  eta3_loop.compute(q_source1,psi1,i) ;
1162  eta4_loop.compute(q_source1,psi1,i) ;
1163  eta4_kilcup_loop.compute(psi1,i, Mass);
1164  eta0_loop.compute(q_source1,psi1,i) ;
1165 
1166  eta4_conn.compute(q_source1, q_source2, psi1, psi2, i) ;
1167 
1168  if(loop_checkpoint){
1169  //write each measurement to the XML file
1170 
1171  scalar_one_loop.dump(xml_out,i) ;
1172  scalar_one_kilcup_loop.dump(xml_out,i);
1173  scalar_two_loop.dump(xml_out,i) ;
1174  scalar_four_loop.dump(xml_out,i) ;
1175  scalar_four_kilcup_loop.dump(xml_out,i) ;
1176  eta3_loop.dump(xml_out,i) ;
1177  eta4_loop.dump(xml_out,i) ;
1178  eta4_kilcup_loop.dump(xml_out,i) ;
1179  eta0_loop.dump(xml_out,i) ;
1180  }
1181 
1182  // MUST checkpoint stochastic connected correlator measurements!
1183  eta4_conn.dump(xml_out,i) ;
1184  QDPIO::cout << __func__ << ": OUTNOW!!!!" << std::endl;
1185  } // Nsamples
1186 
1187 
1188  // write output from the loop calc
1189  scalar_one_loop.dump(xml_out) ;
1190  scalar_one_kilcup_loop.dump(xml_out);
1191  scalar_two_loop.dump(xml_out) ;
1192  scalar_four_loop.dump(xml_out) ;
1193  scalar_four_kilcup_loop.dump(xml_out) ;
1194  eta3_loop.dump(xml_out) ;
1195  eta4_loop.dump(xml_out) ;
1196  eta4_kilcup_loop.dump(xml_out) ;
1197  eta0_loop.dump(xml_out) ;
1198 
1199  // no point in dumping connected correlator measurements,
1200  // since corrs must be formed with each noise source, but what the hell,
1201  // maybe there is a diagnostic use for it.
1202  QDPIO::cout << __func__ << ": OUTNOW!!!!" << std::endl;
1203  eta4_conn.dump(xml_out) ;
1204 
1205  QDPIO::cout << __func__ << ": OUTNOW!!!!" << std::endl;
1206  // end of this section
1207  pop(xml_out);
1208 
1209  }
1210 
1211  /**********************************************************************/
1212 
1213  /**********************************************************************/
1214 
1215  // fuzz the loops
1216 
1219  LatticeStaggeredFermion & q_source,
1220  LatticeStaggeredFermion & psi ,
1221  LatticeStaggeredFermion & psi_fuzz ,
1222  const multi1d<LatticeColorMatrix> & u,
1223  const multi1d<LatticeColorMatrix> & u_smr,
1224  XMLWriter & xml_out,
1225  bool gauge_shift,
1226  bool sym_shift,
1227  bool loop_checkpoint,
1228  int t_length,
1229  Real Mass,
1230  int Nsamp,
1231  Real RsdCG,
1232  int CFGNO,
1233  VolSrc_type volume_source,
1234  int fuzz_width,
1235  int src_seperation,
1236  int j_decay){
1237 
1238  int src_tslice=0;
1239  int src_color_ind = 0;
1240  int src_parity_ind = 0;
1241  int src_corner_ind =0;
1242 
1243  int src_tslice2=0;
1244  int src_color_ind2 = 0;
1245  int src_parity_ind2 = 0;
1246  int src_corner_ind2 =0;
1247 
1248  Real coverage_fraction;
1249  int j;
1250  int fuzz_index;
1251 
1252  LatticeStaggeredFermion psi2 ;
1253  LatticeStaggeredFermion q_source2 ;
1254 
1255  bool fuzz_sink=false;
1256  bool fuzz_src=false;
1257 
1258  push(xml_out,"fuzz_loops_s");
1259 
1260  // write common parameters
1261  write(xml_out, "Mass" , Mass);
1262 
1263 
1264  write_out_source_type(xml_out, volume_source);
1265 
1266 
1267  write(xml_out, "Number_of_samples" , Nsamp);
1268  write(xml_out, "fuzz_width" , fuzz_width);
1269 
1270 
1271  Stag_shift_option type_of_shift ;
1272  if( gauge_shift ){
1273  if(sym_shift){
1274  type_of_shift = SYM_GAUGE_INVAR ;
1275  }else{
1276  type_of_shift = GAUGE_INVAR ;
1277  }
1278  }else{
1279  if(sym_shift){
1280  type_of_shift = SYM_NON_GAUGE_INVAR ;
1281  }else{
1282  type_of_shift = NON_GAUGE_INVAR ;
1283  }
1284  }
1285 
1286  // set up the loop code
1287  local_scalar_loop scalar_one_loop(t_length,Nsamp,
1288  u,type_of_shift) ;
1289  local_scalar_kilcup_loop scalar_one_kilcup_loop(t_length, Nsamp,
1290  u,type_of_shift);
1291  non_local_scalar_loop scalar_two_loop(t_length,Nsamp,
1292  u,type_of_shift) ;
1293  fourlink_scalar_loop scalar_four_loop(t_length, Nsamp,
1294  u, type_of_shift);
1295  fourlink_scalar_kilcup_loop scalar_four_kilcup_loop(t_length, Nsamp,
1296  u, type_of_shift);
1297  threelink_pseudoscalar_loop eta3_loop(t_length,Nsamp,
1298  u,type_of_shift) ;
1299  fourlink_pseudoscalar_loop eta4_loop(t_length,Nsamp,
1300  u,type_of_shift) ;
1301 
1302  fourlink_pseudoscalar_kilcup_loop eta4_kilcup_loop(t_length, Nsamp,
1303  u,type_of_shift) ;
1304 
1305  // for test purposes
1306  zerolink_pseudoscalar_loop eta0_loop(t_length, Nsamp,
1307  u,type_of_shift);
1308 
1309 
1310 
1311  // fuzzed loops here
1312  local_scalar_loop_fuzz scalar_one_loop_fuzz(t_length, Nsamp, u,
1313  type_of_shift);
1314  local_scalar_kilcup_loop_fuzz scalar_one_kilcup_loop_fuzz(t_length, Nsamp,
1315  u, type_of_shift);
1316  non_local_scalar_loop_fuzz scalar_two_loop_fuzz(t_length,Nsamp, u,
1317  type_of_shift);
1318  threelink_pseudoscalar_loop_fuzz eta3_loop_fuzz(t_length,Nsamp,
1319  u,type_of_shift);
1320  fourlink_pseudoscalar_loop_fuzz eta4_loop_fuzz(t_length,Nsamp,
1321  u,type_of_shift);
1322 
1323  fourlink_pseudoscalar_kilcup_loop_fuzz eta4_kilcup_loop_fuzz(t_length,
1324  Nsamp, u,
1325  type_of_shift);
1326 
1327 
1328  // fourlink_pseudoscalar_stoch_conn eta4_conn(t_length,Nsamp,
1329  // u,type_of_shift) ;
1330 
1331 
1332  fourlink_pseudoscalar_stoch_conn eta4_conn(t_length,Nsamp,
1333  u,type_of_shift,
1334  false,false) ;
1335 
1336  fourlink_pseudoscalar_stoch_conn eta4_conn_FsrcLsink(t_length,Nsamp,
1337  u,type_of_shift,
1338  true,false) ;
1339 
1340  fourlink_pseudoscalar_stoch_conn eta4_conn_LsrcFsink(t_length,Nsamp,
1341  u,type_of_shift,
1342  false,true) ;
1343 
1344  // set up the loop code
1345 
1346  // Seed the RNG with the cfg number for now
1347  QDP::Seed seed;
1348  seed = CFGNO;
1349  RNG::setrn(seed);
1350 
1351 
1352  for(int i = 0; i < Nsamp; ++i){
1353  psi = zero; // note this is ``zero'' and not 0
1354  psi2 = zero; // note this is ``zero'' and not 0
1355 
1356  RNG::savern(seed);
1357  QDPIO::cout << "SEED = " << seed << std::endl;
1358 
1359  QDPIO::cout << "Noise sample: " << i << std::endl;
1360 
1361  // Fill the volume with random noise
1362  coverage_fraction = fill_volume_source(q_source, volume_source,
1363  t_length, &src_tslice,
1364  &src_color_ind, &src_parity_ind,
1365  &src_corner_ind, src_seperation,
1366  j_decay);
1367 
1368  //fill 2nd source for stochastic connected correlators
1369  coverage_fraction = fill_volume_source(q_source2, volume_source,
1370  t_length, &src_tslice2,
1371  &src_color_ind2, &src_parity_ind2,
1372  &src_corner_ind2, src_seperation,
1373  j_decay);
1374 
1375 
1376  // Compute the solution std::vector for the particular source
1377  SystemSolverResults_t res = (*qprop)(psi, q_source);
1378  SystemSolverResults_t res2 = (*qprop)(psi2, q_source2);
1379 
1380  push(xml_out,"Qprop_noise");
1381  write(xml_out, "Noise_number" , i);
1382  write(xml_out, "RsdCG" , RsdCG);
1383  write(xml_out, "n_count", res.n_count);
1384  write(xml_out, "Seed" , seed);
1385  pop(xml_out);
1386 
1387  fuzz_smear(u_smr, psi,psi_fuzz, fuzz_width, j_decay) ;
1388 
1389 
1390  // compute the un-fuzzed operators
1391  scalar_one_loop.compute(q_source,psi,i);
1392  scalar_one_kilcup_loop.compute(psi, i, Mass);
1393  scalar_two_loop.compute(q_source,psi,i);
1394  scalar_four_loop.compute(q_source,psi,i);
1395  scalar_four_kilcup_loop.compute(psi,i,Mass);
1396  eta3_loop.compute(q_source,psi,i);
1397  eta4_loop.compute(q_source,psi,i);
1398  eta4_kilcup_loop.compute(psi,i, Mass);
1399  eta0_loop.compute(q_source,psi,i);
1400 
1401  // eta4_conn.compute(q_source, q_source2, psi, psi2, i) ;
1402 
1403  eta4_conn.compute(q_source, q_source2,
1404  psi, psi2,
1405  u_smr, fuzz_width, j_decay, i);
1406 
1407  eta4_conn_FsrcLsink.compute(q_source, q_source2,
1408  psi, psi2,
1409  u_smr, fuzz_width, j_decay, i);
1410  eta4_conn_LsrcFsink.compute(q_source, q_source2,
1411  psi, psi2,
1412  u_smr, fuzz_width, j_decay, i);
1413 
1414  // now compute the fuzzed operators
1415  scalar_one_loop_fuzz.compute(q_source,psi_fuzz,i);
1416  scalar_one_kilcup_loop_fuzz.compute(psi_fuzz, psi, i, Mass);
1417  scalar_two_loop_fuzz.compute(q_source,psi_fuzz,i) ;
1418  eta3_loop_fuzz.compute(q_source,psi_fuzz,i) ;
1419  eta4_loop_fuzz.compute(q_source,psi_fuzz,i) ;
1420  eta4_kilcup_loop_fuzz.compute(psi_fuzz,psi,i, Mass);
1421 
1422 
1423 
1424 
1425  if(loop_checkpoint){
1426  //write each measurement to the XML file
1427 
1428  scalar_one_loop.dump(xml_out,i) ;
1429  scalar_one_kilcup_loop.dump(xml_out,i);
1430  scalar_two_loop.dump(xml_out,i) ;
1431  scalar_four_loop.dump(xml_out,i) ;
1432  scalar_four_kilcup_loop.dump(xml_out,i);
1433  eta3_loop.dump(xml_out,i) ;
1434  eta4_loop.dump(xml_out,i) ;
1435  eta4_kilcup_loop.dump(xml_out,i) ;
1436  eta0_loop.dump(xml_out,i) ;
1437 
1438  scalar_one_loop_fuzz.dump(xml_out,i) ;
1439  scalar_one_kilcup_loop_fuzz.dump(xml_out,i);
1440  scalar_two_loop_fuzz.dump(xml_out,i) ;
1441  eta3_loop_fuzz.dump(xml_out,i) ;
1442  eta4_loop_fuzz.dump(xml_out,i) ;
1443  eta4_kilcup_loop_fuzz.dump(xml_out,i) ;
1444  }
1445 
1446  // MUST checkpoint stochastic connected correlator measurements!
1447  eta4_conn.dump(xml_out,i) ;
1448  eta4_conn_FsrcLsink.dump(xml_out,i) ;
1449  eta4_conn_LsrcFsink.dump(xml_out,i) ;
1450 
1451  } // Nsamples
1452 
1453  // write output from the loop calc
1454  scalar_one_loop.dump(xml_out) ;
1455  scalar_one_kilcup_loop.dump(xml_out);
1456  scalar_two_loop.dump(xml_out) ;
1457  scalar_four_loop.dump(xml_out) ;
1458  scalar_four_kilcup_loop.dump(xml_out) ;
1459  eta3_loop.dump(xml_out) ;
1460  eta4_loop.dump(xml_out) ;
1461  eta4_kilcup_loop.dump(xml_out) ;
1462 
1463  eta0_loop.dump(xml_out) ;
1464 
1465  scalar_one_loop_fuzz.dump(xml_out) ;
1466  scalar_one_kilcup_loop_fuzz.dump(xml_out);
1467  scalar_two_loop_fuzz.dump(xml_out) ;
1468  eta3_loop_fuzz.dump(xml_out) ;
1469  eta4_loop_fuzz.dump(xml_out) ;
1470  eta4_kilcup_loop_fuzz.dump(xml_out) ;
1471 
1472  // no point in dumping connected correlator measurements,
1473  // since corrs must be formed with each noise source, but what the hell,
1474  // maybe there is a diagnostic use for it.
1475  QDPIO::cout << __func__ << ": OUTNOW!!!!" << std::endl;
1476  eta4_conn.dump(xml_out) ;
1477  eta4_conn_FsrcLsink.dump(xml_out) ;
1478  eta4_conn_LsrcFsink.dump(xml_out) ;
1479 
1480  // end of this section
1481  pop(xml_out);
1482 
1483  }
1484 
1485 
1486  /**********************************************************************/
1487 
1488 
1489 } // end namespace Chroma
Class for counted reference semantics.
Definition: handle.h:33
Class for local (zero-link) scalar loop, with VKVR and fuzzing.
void compute(LatticeStaggeredFermion &psi_fuzz, LatticeStaggeredFermion &psi, int isample, Real mass)
Do the measurement (with VKVR)
Class for local (zero-link) scalar loop using VKVR trick.
void compute(LatticeStaggeredFermion &psi, int isample, Real mass)
Do the measurement (needs psi, not q_source!)
void compute(LatticeStaggeredFermion &q_source, LatticeStaggeredFermion &psi, int isample)
void compute(LatticeStaggeredFermion &q_source, LatticeStaggeredFermion &psi, int isample)
void compute(LatticeStaggeredFermion &q_source, LatticeStaggeredFermion &psi, int isample)
void compute(LatticeStaggeredFermion &q_source, LatticeStaggeredFermion &psi, int isample)
void dump(XMLWriter &xml_out)
Definition: loops_s.h:29
void binary_dump(std::string start_name)
Definition: loops_s.h:155
void dump(XMLWriter &xml_out)
Definition: stoch_conn_s.h:32
Dilute Gaussian sources.
Class structure for fermion actions.
Fuzzed sources.
void read(XMLReader &xml, const std::string &path, AsqtadFermActParams &param)
Read parameters.
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
void gaussian_on_timeslice(LatticeStaggeredFermion &a, int slice, int mu)
Volume source of complex Z2 noise.
void fuzz_smear(const multi1d< LatticeColorMatrix > &ufuzz, const T &psi, T &psifuzz, int length, int j_decay)
apply a fuzz_smear operator to a lattice field
Definition: fuzz_smear.cc:75
void gaussian_on_corner(LatticeStaggeredFermion &a, int corner_index)
Diluted Gauusian-source.
void z2_src(LatticeFermion &a)
Z2-source.
Definition: z2_src.cc:53
void gaussian_corner_on_mod_dbl_slice(LatticeStaggeredFermion &a, int corner_index, int slice, int mu, int seperation)
Diluted Gauusian-source.
void gaussian_on_parity(LatticeStaggeredFermion &a, int parity)
Diluted Gauusian-source.
void gaussian_color_src_on_slice(LatticeStaggeredFermion &a, int color_index, int slice, int mu)
Diluted Gauusian-source.
void gaussian_color_src_on_parity(LatticeStaggeredFermion &a, int color_index, int parity)
Diluted Gauusian-source.
void gaussian_color_src(LatticeStaggeredFermion &a, int color_index)
Diluted Gauusian-source.
void gaussian_parity_src_on_slice(LatticeStaggeredFermion &a, int parity, int slice, int mu)
Diluted Gauusian-source.
void gaussian_on_mod_timeslice(LatticeStaggeredFermion &a, int slice, int mu, int seperation)
Diluted Gauusian-source.
void gaussian_corner_on_dbl_slice(LatticeStaggeredFermion &a, int corner_index, int slice, int mu)
Diluted Gauusian-source.
void gaussian_color_src_on_mod_slice(LatticeStaggeredFermion &a, int color_index, int slice, int mu, int seperation)
Diluted Gauusian-source.
Hadronic observables.
Wrapper routines for computing loops with staggeref fermions.
unsigned j
Definition: ldumul_w.cc:35
void setrn(int iseed[4])
Definition: make_seeds.cc:38
void savern(int iseed[4])
Definition: make_seeds.cc:46
int j_decay
Definition: meslate.cc:22
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
gaussian(aux)
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > enum InvType invType const multi1d< Real > & RsdCG
Definition: pbg5p_w.cc:30
QDP_error_exit("too many BiCG iterations", n_count, rsd_sq, cp, c, re_rvr, im_rvr, re_a, im_a, re_b, im_b)
static multi1d< LatticeColorMatrix > u
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > & Mass
Definition: pbg5p_w.cc:29
push(xml_out,"Condensates")
void ks_local_loops_and_stoch_conn(Handle< SystemSolver< LatticeStaggeredFermion > > &qprop, LatticeStaggeredFermion &q_source1, LatticeStaggeredFermion &psi1, const multi1d< LatticeColorMatrix > &u, XMLWriter &xml_out, bool gauge_shift, bool sym_shift, bool loop_checkpoint, int t_length, Real Mass, int Nsamp, Real RsdCG, int CFGNO, VolSrc_type volume_source, int src_seperation, int j_decay)
@ SYM_NON_GAUGE_INVAR
int i
Definition: pbg5p_w.cc:55
void write_out_source_type(XMLWriter &xml_out, VolSrc_type volume_source)
void ks_fuzz_loops_X(Handle< SystemSolver< LatticeStaggeredFermion > > &qprop, LatticeStaggeredFermion &q_source, LatticeStaggeredFermion &psi, LatticeStaggeredFermion &psi_fuzz, const multi1d< LatticeColorMatrix > &u, const multi1d< LatticeColorMatrix > &u_smr, XMLWriter &xml_out, bool gauge_shift, bool sym_shift, bool loop_checkpoint, int t_length, Real Mass, int Nsamp, Real RsdCG, int CFGNO, VolSrc_type volume_source, int fuzz_width, int src_seperation, int j_decay)
void ks_fuzz_loops_stoch_conn(Handle< SystemSolver< LatticeStaggeredFermion > > &qprop, LatticeStaggeredFermion &q_source, LatticeStaggeredFermion &psi, LatticeStaggeredFermion &psi_fuzz, const multi1d< LatticeColorMatrix > &u, const multi1d< LatticeColorMatrix > &u_smr, XMLWriter &xml_out, bool gauge_shift, bool sym_shift, bool loop_checkpoint, int t_length, Real Mass, int Nsamp, Real RsdCG, int CFGNO, VolSrc_type volume_source, int fuzz_width, int src_seperation, int j_decay)
LatticeFermion psi
Definition: mespbg5p_w.cc:35
Real fill_volume_source(LatticeStaggeredFermion &q_source, VolSrc_type volume_source, int t_length, int *p_src_tslice, int *p_src_color_ind, int *p_src_parity_ind, int *p_src_corner_ind, int src_seperation, int j_decay)
void ks_fuzz_loops(Handle< SystemSolver< LatticeStaggeredFermion > > &qprop, LatticeStaggeredFermion &q_source, LatticeStaggeredFermion &psi, LatticeStaggeredFermion &psi_fuzz, const multi1d< LatticeColorMatrix > &u, const multi1d< LatticeColorMatrix > &u_smr, XMLWriter &xml_out, bool gauge_shift, bool sym_shift, bool loop_checkpoint, int t_length, Real Mass, int Nsamp, Real RsdCG, int CFGNO, VolSrc_type volume_source, int fuzz_width, int src_seperation, int j_decay, bool binary_loop_checkpoint, std::string binary_name)
pop(xml_out)
void ks_local_loops(Handle< SystemSolver< LatticeStaggeredFermion > > &qprop, LatticeStaggeredFermion &q_source, LatticeStaggeredFermion &psi, const multi1d< LatticeColorMatrix > &u, XMLWriter &xml_out, bool gauge_shift, bool sym_shift, bool loop_checkpoint, int t_length, Real Mass, int Nsamp, Real RsdCG, int CFGNO, VolSrc_type volume_source, int src_seperation, int j_decay)
@ T_DILUTE_GAUSS
Definition: enum_loops_s.h:19
@ C_DILUTE_GAUSS
Definition: enum_loops_s.h:20
@ C_MOD_T_DILUTE_GAUSS
Definition: enum_loops_s.h:29
@ COR_MOD_DBL_T_DILUTE_GAUSS
Definition: enum_loops_s.h:28
@ CORNER_DILUTE_GAUSS
Definition: enum_loops_s.h:26
@ GAUSSIAN
Definition: enum_loops_s.h:18
@ MOD_T_DILUTE_GAUSS
Definition: enum_loops_s.h:25
@ PT_DILUTE_GAUSS
Definition: enum_loops_s.h:24
@ COR_DBL_T_DILUTE_GAUSS
Definition: enum_loops_s.h:27
@ P_DILUTE_GAUSS
Definition: enum_loops_s.h:21
@ CP_DILUTE_GAUSS
Definition: enum_loops_s.h:23
@ CT_DILUTE_GAUSS
Definition: enum_loops_s.h:22
Double zero
Definition: invbicg.cc:106
::std::string string
Definition: gtest.h:1979
Point source construction.
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Volume source of Z2 noise.