CHROMA
lg5eps_double_pass_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Overlap-pole operator
3  */
4 
5 
6 #include "chromabase.h"
8 #include "meas/eig/gramschm.h"
9 
10 
11 namespace Chroma
12 {
13 //! Apply the GW operator onto a source std::vector
14 /*! \ingroup linop
15  *
16  * This routine applies the 4D GW operator onto a source
17  * std::vector. The coeffiecients for the approximation get
18  * wired into the class by the constructor and should
19  * come fromt fermion action.
20  *
21  * The operator applied is:
22  * D = (1/2)[ (1+m) + (1-m)gamma_5 sgn(H_w) ] psi
23  * or D^{dag} = (1/2)[ (1+m) + (1-m) sgn(H_w) gamma_5 psi
24  *
25  *
26  * \param chi result std::vector (Write)
27  * \param psi source std::vector (Read)
28  * \param isign Hermitian Conjugation Flag
29  * ( PLUS = no dagger| MINUS = dagger ) (Read)
30  */
31 void lg5eps_double_pass::operator() (LatticeFermion& chi,
32  const LatticeFermion& psi,
33  enum PlusMinus isign) const
34 {
36 }
37 
38 
39 //! Apply the GW operator onto a source std::vector
40 /*! \ingroup linop
41  *
42  * This routine applies the 4D GW operator onto a source
43  * std::vector. The coeffiecients for the approximation get
44  * wired into the class by the constructor and should
45  * come fromt fermion action.
46  *
47  * The operator applied is:
48  * D = (1/2)[ (1+m) + (1-m)gamma_5 sgn(H_w) ] psi
49  * or D^{dag} = (1/2)[ (1+m) + (1-m) sgn(H_w) gamma_5 psi
50  *
51  *
52  * \param chi result std::vector (Write)
53  * \param psi source std::vector (Read)
54  * \param isign Hermitian Conjugation Flag
55  * ( PLUS = no dagger| MINUS = dagger ) (Read)
56  */
57 void lg5eps_double_pass::operator() (LatticeFermion& chi,
58  const LatticeFermion& psi,
59  enum PlusMinus isign, Real epsilon) const
60 {
61  START_CODE();
62 
63  LatticeFermion tmp1, tmp2;
64 
65  // Gamma_5
66  int G5 = Ns*Ns - 1;
67 
68  chi = zero;
69 
70  switch (isign)
71  {
72  case PLUS:
73  // Non-Dagger: psi is source and tmp1
74  // chi := gamma_5 * (gamma_5 * mass + sgn(H)) * Psi
75  tmp1 = psi;
76  break;
77 
78  case MINUS:
79  // Dagger: apply gamma_5 to source psi to make tmp1
80  // chi := (mass + sgn(H) * gamma_5) * Psi
81  tmp1 = Gamma(G5) * psi;
82  break;
83 
84  default:
85  QDP_error_exit("unknown isign value", isign);
86  }
87 
88 
89 
90  // Project out eigenvectors of source if desired
91  // chi += func(lambda) * EigVec * <EigVec, psi>
92  // Usually "func(.)" is sgn(.); it is precomputed in EigValFunc.
93  // for all the eigenvalues
94  if (NEig > 0)
95  {
96  Complex cconsts;
97 
98  for(int i = 0; i < NEig; ++i)
99  {
100  cconsts = innerProduct(EigVec[i], tmp1);
101  tmp1 -= EigVec[i] * cconsts;
102 
103  cconsts *= EigValFunc[i];
104  chi += EigVec[i] * cconsts;
105  }
106  }
107 
108  // tmp1 <- H * Projected psi
109  // <- gamma_5 * M * psi
110  (*M)(tmp2, tmp1, PLUS);
111  tmp1 = Gamma(G5) * tmp2;
112 
113  Double c = norm2(tmp1);
114 
115  /* If exactly 0 norm, then solution must be 0 (for pos. def. operator) */
116  if (toBool(c == 0))
117  {
118  chi = zero;
119  END_CODE();
120  return;
121  }
122 
123  // Multiply in P(0) -- this may well be 0 for type 0 rational approximations
124  chi += tmp1 * constP;
125 
126 
127  // Effectively the multi mass solve is on tmp1
128  // *******************************************************************
129  // Solve (MdagM + rootQ_n) chi_n = H * tmp1
130  LatticeFermion Ap;
131  LatticeFermion r;
132  LatticeFermion p;
133 
134 
135  Double rsd_sq = norm2(psi) * epsilon*epsilon; // || psi ||^2 * epsilon^2
136 
137  Double cp;
138  Double d; // InnerProduct
139 
140 
141 
142 
143  multi1d<Double> a(MaxCG+1); // Alpha for unshifted (isz) system
144  multi1d<Double> b(MaxCG+1); // Beta for unshifted system
145 
146  bool convP;
147  int iters_taken;
148 
149  // By default, rootQ(isz) is considered the smallest shift
150  int isz = numroot-1; // isz identifies system with smalles shift
151 
152 
153 
154 
155  // r[0] := p[0] := tmp1
156  r = tmp1;
157  p = tmp1;
158 
159 
160  cp = 0;
161  convP = false;
162 
163 
164  Double alpha_minus_one=1;
165  b[0] = 0; // b[0] -- do we ever use ?
166 
167 
168 
169  int k;
170  multi1d<Double> c_iter(MaxCG+1);
171 
172  c_iter[0] = c;
173  // Do the iterations proper
174  // first pass
175  for(k = 0; k < MaxCG && ! convP ; ++k) {
176 
177  // Keep hold of the residuum
178 
179  (*MdagM)(Ap, p, PLUS);
180  Ap += p * rootQ[isz];
181 
182  // Project out eigenvectors
183  if (k % ReorthFreq == 0) {
184  GramSchm(Ap, EigVec, NEig, all);
185  }
186 
187  // d = < p, A.p >
188  d = innerProductReal(p, Ap); // 2 Nc Ns flops
189 
190 
191  // k+1 because a[0] corresponds to a_{-1}
192  a[k] = c/d;
193  r -= Real(a[k])*Ap;
194 
195  // Project out eigenvectors
196  if (k % ReorthFreq == 0) {
197  GramSchm (r, EigVec, NEig, all);
198  }
199 
200  cp = c;
201  c = norm2(r);
202 
203  c_iter[k+1] = c;
204  b[k+1] = c/cp;
205 
206  p = r + Real(b[k+1])*p;
207 
208  convP = toBool( c < rsd_sq );
209  }
210 
211  int niters = k;
212 
213  // OK First pass done. I now need to compute gamma_j and c_j
214  multi2d<Double> gamma(niters+1, numroot);
215  for(int s = 0; s < numroot; s++) {
216  gamma[0][s] = 1; // Really gamma_0
217  }
218 
219 
220  multi1d<bool> convPs(numroot);
221  multi1d<int> convIter(numroot);
222  convPs = false;
223 
224  for(int j=0; j < niters; j++) {
225  for(int s =0; s < numroot; s++) {
226 
227 
228  if( s != isz && ! convPs[s] ) {
229  Double a_minus;
230 
231  if( j == 0 ) {
232  // alpha_minus_one -- no need to store. Only need alpha_0,...
233  // in second pass
234  a_minus = Double(1);
235  }
236  else {
237  a_minus = a[j-1];
238  }
239 
240  Double ga = gamma[j][s];
241  Double ga_minus;
242  if( j == 0 ) {
243  // gamma[-1][s] -- no need to store. Only need gamma[0][s]...
244  // in second pass
245  ga_minus =Double(1);
246  }
247  else{
248  ga_minus = gamma[j-1][s];
249  }
250 
251  Double tmp_num = ga*ga_minus*a_minus;
252  Double tmp_den = ga_minus*a_minus*(Double(1) + a[j]*(rootQ[s] - rootQ[isz]));
253  tmp_den += a[j]*b[j]*(ga_minus - ga);
254 
255  gamma[j+1][s] = tmp_num/tmp_den;
256 
257  // If this system has converged at iter j+1, then dont update
258  // gamma-s anymore and note the convergence point. Updating
259  // ad infinitum causes underflow.
260  if( toBool( gamma[j+1][s]*gamma[j+1][s]*c_iter[j+1] < rsd_sq ) ) {
261  convPs[s] = true;
262  convIter[s] = j+1;
263  }
264 
265  }
266 
267 
268  }
269 
270 
271  }
272 
273 
274  multi1d<Double> sumC(niters+1);
275 
276  for(int j=0; j<=niters; j++) {
277 
278  sumC[j] = 0;
279 
280  for(int m=0; m < niters-j; m++) {
281 
282  Double qsum=resP[isz];
283 
284  for(int s =0; s < numroot; s++) {
285 
286  if( s != isz ) {
287 
288  // Only add gammas which are unconverged.
289  // Converged gamma-s are not updated beyond convergence
290  // which would cause problems (with underflows)
291  // gamma[convIters[s]] is last valid gamma
292  if( toBool( (j+m+1) <= convIter[s] ) ) {
293  qsum += resP[s]*gamma[m+j+1][s]*gamma[m+j][s]/gamma[j][s];
294  }
295 
296  }
297  }
298 
299 
300  Double delta_m = Double(1);
301 
302  if( m > 0 ) {
303  for(int k=1; k <= m; k++) {
304  delta_m *= b[j+k];
305  }
306  }
307 
308 
309  sumC[j] += a[j+m]*delta_m*qsum;
310  }
311  }
312 
313 
314  // Second pass Lanczos
315 
316  // Initialise r_0, p_0
317  r=tmp1;
318  p=tmp1;
319 
320  convP = false;
321 
322  for(k=0; k < niters && !convP ; k++) {
323  (*MdagM)(Ap, p, PLUS);
324  Ap += p * rootQ[isz];
325 
326  // Project out eigenvectors
327  if (k % ReorthFreq == 0) {
328  GramSchm(Ap, EigVec, NEig, all);
329  }
330 
331 
332  // Update chi
333  chi += Real(sumC[k])*r;
334 
335  // Update r
336  r -= Real(a[k])*Ap;
337 
338  // Project out eigenvectors
339  if (k % ReorthFreq == 0) {
340  GramSchm (r, EigVec, NEig, all);
341  }
342 
343 #if 0
344 
345  // This early termination criterion has to be abandoned when
346  // using relaxed solevers
347  Double chi_norm_new = epsilon*epsilon*norm2(chi);
348 
349  // Convergence criterion for total signum. Might be good enough
350  // without running to full niters
351  convP = toBool( c_iter[k+1]*sumC[k+1]*sumC[k+1] < chi_norm_new );
352 #endif
353  p = r + Real(b[k+1])*p;
354  }
355 
356 
357  // End of MULTI SHIFTERY
358 
359  // Now fix up the thing. Multiply in gamma5 if needed
360  // and then rescale to correct normalisation.
361 
362  if (isign == PLUS)
363  {
364  // chi := gamma_5 * (gamma_5 * mass + eps(H)) * Psi
365  tmp1 = Gamma(G5) * chi;
366  chi = tmp1;
367  }
368 
369  QDPIO::cout << "Overlap Inner Solve (lg5eps_double_pass): " << k << " iterations " << std::endl;
370 
371  END_CODE();
372 }
373 
374 } // End Namespace Chroma
375 
376 
Primary include file for CHROMA library code.
const multi1d< Real > EigValFunc
const multi1d< Real > rootQ
const multi1d< LatticeFermion > EigVec
const multi1d< Real > resP
Gramm-Schmidt orthogonolization.
void GramSchm(multi1d< LatticeFermion > &psi, const int Npsi, const multi1d< LatticeFermion > &vec, const int Nvec, const Subset &sub)
Gramm-Schmidt orthogonolization.
Definition: gramschm.cc:127
void operator()(LatticeFermion &chi, const LatticeFermion &psi, enum PlusMinus isign) const
Apply the operator onto a source std::vector.
unsigned j
Definition: ldumul_w.cc:35
Internal Overlap-pole operator.
static int m[4]
Definition: make_seeds.cc:16
Double tmp2
Definition: mesq.cc:30
int epsilon(int i, int j, int k)
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
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)
int G5
Definition: pbg5p_w.cc:57
Double c
Definition: invbicg.cc:108
Real rsd_sq
Definition: invbicg.cc:121
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
LatticeFermion psi
Definition: mespbg5p_w.cc:35
DComplex d
Definition: invbicg.cc:99
START_CODE()
Complex b
Definition: invbicg.cc:96
Double zero
Definition: invbicg.cc:106
int k
Definition: invbicg.cc:119
multi1d< LatticeFermion > s(Ncb)
FloatingPoint< double > Double
Definition: gtest.h:7351
int isz
Definition: pade_trln_w.cc:151