CHROMA
pbg5p_w.cc
Go to the documentation of this file.
1 namespace Chroma {
2 
3 
4 /* $Id: pbg5p_w.cc,v 3.0 2006-04-03 04:59:04 edwards Exp $ */
5 
6 /* PBG5P - Calculates noise estimator for the trace used in the
7 /* following constructions */
8 
9 /* <psibar_psi(x)> = Tr D^{-1} */
10 /* <psibar psi(x) psibar psi(y)> = Tr D^{-2} - Tr D{-1} */
11 /* <psibar gamma_5 psi psibar gamma_5 psi> = = Tr H^{-2} */
12 
13 /* What is computed are the traces so iso(flavor) triplets and singlets can be
14 /* reconstructed */
15 
16 /* This routine is specific to Wilson and really Overlap fermions! */
17 
18 /* u -- gauge field ( Read ) */
19 /* nhit -- number of Gaussian hits ( Read ) */
20 /* Mass -- array of quark mass values ( Read ) */
21 /* numMass -- number of quark masses ( Read ) */
22 /* RsdCG -- the CG accuracy used for multishift CG ( Read ) */
23 
24 include(types.mh)
25 
26 void Pbg5p(XMLWriter& xml_out,
27  const WilsonTypeFermAct< multi1d<LatticeFermion> >& S_f,
29  const multi1d<Real>& Mass,
30  enum InvType invType,
31  const multi1d<Real>& RsdCG,
32  int MaxCG,
33  int nhit)
34 {
35  START_CODE();
36 
37  int n_congrd;
39  multi1d<Double> TrDinv(numMass);
40  multi1d<Double> TrHinv_sq(numMass);
41  multi1d<Double> TrDinv_avg(numMass);
42  multi1d<Double> TrHinv_sq_avg(numMass);
43  LatticeFermion eta;
44  LatticeFermion eta_sq;
45  LatticeFermion tmp;
46  LatticeFermion b;
47  LatticeFermion b_orig;
48  LatticeReal lrtrace_aux;
49 
52  Double one;
55  int i;
56 
57  int G5 = Ns*Ns - 1;
59 
60 
63  n_congrd = 0;
64  n_congrd_tot = 0;
65 
66 
67  for(int ihit = 1; ihit <= nhit; ++ihit)
68  {
69  /* Fill with random gaussian noise such that < b_dagger * b > = 1 per d.o.f. */
70  gaussian(b);
71 
72 #if 1
73  /* If using Schroedinger functional, zero out the boundaries */
74  S_f.getFermBC().modifyF(b);
75 #endif
76 
77  /* Make the source chiral */
78  LatticeFermion tmp = chiralProjectPlus(b);
79  b = tmp;
80 
81  /* Prop. solution */
82  /* eta = D^(-1) * b */
83  /* eta_sq = (D^dag*D)^(-1) * b */
84  /* The inverse of the un-preconditioned matrix is needed here. This */
85  /* is given by Qprop! */
86  /* Construct */
87  /* Tr Dinv = < b_dag * eta > */
88  /* Tr Hinv_sq = < b_dag * eta_sq > */
89 
90  TrDinv = 0;
91  TrHinv_sq = 0;
92 
93  n_congrd = 0;
94  eta = 0;
95  eta_sq = 0;
96  b_orig = b;
97 
98  for(i=0; i < numMass; ++i)
99  {
100  gaussian(tmp);
101  CHIRAL_PROJECT(tmp, isign, eta, REPLACE);
102  CHIRAL_PROJECT(tmp, isign, eta_sq, REPLACE);
103 
104  Qprop (u, b, Mass[i], RsdCG[i], eta, n_congrd);
105  b = b_orig;
106 
107  tmp = Gamma(G5)*eta;
108  Qprop (u, tmp, Mass[i], RsdCG[i], eta_sq, n_congrd);
109  tmp = Gamma(G5)*eta_sq;
110  eta_sq = tmp;
111 
112  TrDinv[i] += real(innerProduct(b, eta));
113  TrHinv_sq[i] += real(innerProduct(b, eta_sq));
114  }
115 
116 
117  /* Normalize and print */
118  ddummy1 = Double(1) / Double(vol);
119  for(i=0; i < numMass; ++i)
120  {
121  TrDinv[i] *= ddummy1;
122  TrHinv_sq[i] *= ddummy1;
123  }
124 
125  if (nhit > 1)
126  {
127  push(xml_out,"Pbpg5_trace_hit");
128  write(xml_out, "ihit", ihit);
129  write(xml_out, "nhit", nhit);
130  write(xml_out, "TrDinv", TrDinv);
131  write(xml_out, "TrHinv_sq", TrHinv_sq);
132  write(xml_out, "n_congrd", n_congrd);
133  pop(xml_out);
134  }
135 
136  /* Add on to running average */
137  TrDinv_avg += TrDinv;
140  }
141 
142 
143  /* Normalize */
144  ddummy1 = Double(1) / Double(nhit);
145  for(i=0; i < numMass; ++i)
146  {
147  TrDinv_avg[i] *= ddummy1;
148  TrHinv_sq_avg[i] *= ddummy1;
149  }
150 
151  push(xml_out,"Pbp5g_trace_avg");
152  write(xml_out, "nhit", nhit);
153  write(xml_out, "Mass", Mass);
154  write(xml_out, "TrDinv_avg", TrDinv_avg);
155  write(xml_out, "TrHinv_sq_avg", TrHinv_sq_avg);
156  write(xml_out, "n_congrd_tot", n_congrd_tot);
157  pop(xml_out);
158 
159 
160  END_CODE();
161 }
162 
163 } // end namespace Chroma
Support class for fermion actions and linear operators.
Definition: state.h:28
Class for counted reference semantics.
Definition: handle.h:33
Wilson-like fermion actions.
Definition: fermact.orig.h:344
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
BinaryReturn< C1, C2, FnInnerProduct >::Type_t innerProduct(const QDPSubType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
gaussian(aux)
Double two
Definition: pbg5p_w.cc:53
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > enum InvType invType const multi1d< Real > & RsdCG
Definition: pbg5p_w.cc:30
int G5
Definition: pbg5p_w.cc:57
static multi1d< LatticeColorMatrix > u
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > & Mass
Definition: pbg5p_w.cc:29
TrHinv_sq_avg
Definition: pbg5p_w.cc:62
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
int n_congrd
Definition: mespbg5p_w.cc:24
Double one
Definition: invbicg.cc:105
LatticeFermion b_orig
Definition: pbg5p_w.cc:47
push(xml_out,"Condensates")
Double ddummy1
Definition: mespbg5p_w.cc:44
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > enum InvType invType const multi1d< Real > int MaxCG
Definition: pbg5p_w.cc:32
int i
Definition: pbg5p_w.cc:55
Double mass
Definition: pbg5p_w.cc:54
@ PLUS
Definition: chromabase.h:45
int ihit
Definition: mespbg5p_w.cc:46
LatticeReal lrtrace_aux
Definition: mespbg5p_w.cc:40
pop(xml_out)
LatticeFermion eta
Definition: mespbg5p_w.cc:37
int n_congrd_tot
Definition: pbg5p_w.cc:38
START_CODE()
LatticeFermion eta_sq
Definition: pbg5p_w.cc:44
TrDinv_avg
Definition: pbg5p_w.cc:61
multi1d< Double > TrHinv_sq(numMass)
Complex b
Definition: invbicg.cc:96
int nhit
Definition: mespbg5p_w.cc:26
include(types.mh) LINEAR_OPERATOR(A)
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
multi1d< Double > TrDinv(numMass)
const WilsonTypeFermAct< multi1d< LatticeFermion > > & S_f
Definition: pbg5p_w.cc:27
Double ddummy2
Definition: mespbg5p_w.cc:45
FloatingPoint< double > Double
Definition: gtest.h:7351