CHROMA
pade_trln_w.cc
Go to the documentation of this file.
1 
2 #error "NOT FULLY CONVERTED"
3 
4 
5 include(types.mh)
6 
7 
8 /* Parameters:
9 
10  u -- The Gauge Field (Read)
11  chi -- Z2 noise std::vector (Read)
12  Ncb -- No of checkerboards (Read)
13  Nkappa -- The no of kappa values (Read)
14  Npade -- The order of the Pade approx (Read)
15  Kappas -- The array of Kappa values (Read)
16  Pade_c -- Array holding the denominator Pade coeffs (Read)
17  Pade_b0 -- b0 from the Pade Expansion (Read)
18  Pade_b -- Array holding numerator Pade coeffs (Read)
19  RsdMR -- Stopping Criterion for the MR algorithm (Read)
20  psibp -- An array to hold the psibar-psi values (Write)
21  trlns -- An array to hold the trace log values (Write)
22  n_count -- The number of MR iterations taken (n_count)
23 */
24 SUBROUTINE(PadeTrLn, u, chi, Ncb, Nkappa, Kappas, Npade,
27  n_count)
28 
29  /* Gauge field */
30  multi1d<LatticeColorMatrix> u(Nd);
31 
32  /* Z2 src */
33  multi1d<LatticeFermion> chi(Ncb);
34 
35  /* No of checkerboards */
36  int Ncb;
37 
38  /* No of Kappa values */
39  int Nkappa;
40 
41  /* No of Pade Coeffs */
42  int Npade;
43 
44  /* The array of Kappas */
45  multi1d<Real> Kappas(Nkappa);
46 
47  /* Index of largest kappa (smallest mass) */
48  /* The array of Pade_c */
49  multi1d<Real> Pade_c(Npade);
50  Real Pade_b0;
51  multi1d<Real> Pade_b(Npade);
52 
53  /* Stopping criterion for the Multi Mass MR solver */
54  Real RsdMR;
55 
56  /* Array for unimproved psibar psi estimates */
57  multi1d<Complex> psibp_unimp(Nkappa);
58 
59  /* Array of unimproved Trace Log estimates */
60  multi1d<Complex> trlns_unimp(Nkappa);
61 
62  /* Array for improved psibar psi estimates */
63  multi1d<Complex> psibp_imp(Nkappa);
64 
65  /* Array of improved Trace Log estimates */
66  multi1d<Complex> trlns_imp(Nkappa);
67 
68  /* No of solver iters taken */
69  int n_count;
70 
71  /* Order of unbiased subtraction */
72  int Nimpr;
73  multi1d<Real> Omega_impr(Nimpr);
74 
75 {
76  include(COMMON_DECLARATIONS)
77 
78  multi2d<LatticeFermion> psi(Ncb, Nshift); /* Array of solver solutions */
79 
80  multi1d<LatticeFermion> psi_test(Ncb); /* SOmething to check the residuals with */
81  multi1d<LatticeFermion> r(Ncb);
82 
83  multi1d<LatticeFermion> Dchi(Ncb);
84 
85  Double chi_norm; /* The norm of the source */
86  Double r_norm; /* The norm of the residuals */
87 
88  /* The shifted and unshifted operators */
89 
90  /* A dslash operator */
92 
93  LINEAR_OPERATOR(Mshifted);
94  PROTOTYPE(`Mshifted', `DESC', `DESC', `DESC', `VAL', `VAL')
95 
96  LINEAR_OPERATOR(Munshifted);
97  PROTOTYPE(`Munshifted', `DESC', `DESC', `DESC', `VAL', `VAL')
98 
99  /* The number of shifted solutions */
100  int Nshift; /* No of shifts for solver */
101  multi1d<Real> Shifts(Nshift); /* The array of shifts */
102 
103  multi1d<Real> RsdMR_array(Nshift); /* An array version of RsdMR needed by solver */
104 
105  Real Shiftmass; /* Useful for testing the solutions */
106 
107  Real MaxRsd;
108  /* Various Counters */
109  int i;
110  int j;
111  int l;
112  int kappa;
113  int cb;
114 
115  /* Two variables for tracing */
116  LatticeComplex lc;
117  LatticeComplex lct;
118  DComplex dcsum;
119  Complex csum;
120 
121  /* Temporaries */
122  multi1d<Real> first_terms(Nkappa);
123  multi1d<Complex> pade_terms(Npade);
124 
125 
126  /* Stuff for Trln */
127  Real TrI;
128  Real minusone;
129  Real tmp;
130  Real kappa_pow;
132 
133 
134  /* Dslash Bilinears */
135  multi1d<Complex> bilinear(Nimpr);
136 
137  /* Psibar Psi improvement terms */
138  multi2d<Complex> psib_offterms(Nkappa, Nimpr);
139  multi3d<Complex> trln_offterms(Nkappa, Npade, Nimpr);
140  multi2d<Complex> sum_trln_offterms(Nkappa, Nimpr);
141 
142  /* Glue improvement : 'Plaquette' */
147  Real f_plaq;
148  Complex c_plaq;
149 
150  /* index of the smallest shift */
151  int isz;
152  Real minshift;
153 
155 
156  minusone = - TO_REAL(1);
157 
158  /* First work out the shifts from the kappas and the sigmas */
159  Nshift = Nkappa * ( Npade + 1 );
160 
161  /* Allocate the space */
162 
163  /* Work out the shifts */
164  j=0;
165 
166  /* First just the 1/kappa's */
167  for(i = 0; i < Nkappa; i++) {
168  Shifts[j] = TO_REAL(1) / Kappas[i];
169  j++;
170  }
171 
172  /* Then the (1 + shift)/kappa's */
173  for(i = 0; i < Nkappa; i++) {
174  for(l = 0; l < Npade; l++) {
175  Shifts[j] = ( TO_REAL(1) + Pade_c[l] ) / Kappas[i];
176  j++;
177  }
178  }
179 
180  /* Set up the solver */
181 
183  /* push(xml_out,"pade_trln");
184 write(xml_out, "Nshift", Nshift);
185 write(xml_out, "Shifts", Shifts);
186 write(xml_out, "RsdMR", RsdMR);
187 write(xml_out, "MaxCG", MaxCG);
188 pop(xml_out); */
189 
190  /* Fill all Psi-s with 0 to start with */
191  psi = 0;
192 
193  /* Compute the norm of chi */
194  chi_norm=TO_REAL(0);
195 
196  for(cb=0; cb < Ncb; cb++) {
197  chi_norm += norm2(chi[cb]);
198  }
199  chi_norm = sqrt(chi_norm);
200 
201  /* Flick on fermionic Bc's */
203 
204 
205  Shiftmass = 0;
207 
208  CONSTRUCT_LINEAR_OPERATOR(Munshifted, ldunshifted, D);
209 
210  /* Work out what the smallest shift */
211  minshift = Shifts[0];
212  isz = 0;
213  for(i=1; i < Nshift; i++) {
214  if( Shifts[i] < minshift ) {
215  minshift = Shifts[i];
216  isz = i;
217  }
218  }
219 
220  /* Call the solver */
222 
223  /* push(xml_out,"Pade_MR_Solver");
224 write(xml_out, "n_count", n_count);
225 pop(xml_out); */
226 
227  FREE_LINEAR_OPERATOR(Munshifted);
228 
229 
230  for(i=0; i < Nshift; i++) {
231 
232  for(cb = 0; cb < Ncb; cb++) {
233  psi_test[cb] = psi[i][cb];
234  }
235 
236 
237  Shiftmass = Shifts[i];
238  CONSTRUCT_LINEAR_OPERATOR(Mshifted, ldshifted, D, Shiftmass);
239 
240  Mshifted (Mshifted, psi_test, r, Ncb, PLUS);
241 
242  FREE_LINEAR_OPERATOR(Mshifted);
243 
244  Shiftmass = - TO_REAL(1);
245  for(cb = 0; cb < Ncb; cb++) {
246  r[cb] = r[cb] * Shiftmass;
247  r[cb] = chi[cb] + r[cb];
248  }
249 
250 
251  r_norm = 0;
252  for(cb = 0; cb < Ncb; cb++) {
253  r_norm += norm2(r[cb]);
254  }
255 
256  r_norm = sqrt(r_norm);
257  r_norm /= chi_norm;
258  RsdMR_array[i] = FLOAT(r_norm);
259  }
260 
261  MaxRsd = RsdMR_array[0];
262  for( i=1; i < Nshift; i++) {
263  if( RsdMR_array[i] > MaxRsd ) {
264  MaxRsd = RsdMR_array[i];
265  }
266  }
267  PRINTF("TrLnSolve: n_count = %d, RsdMR = %16.8e, MaxRsd = %16.8e\n",
268  n_count, RsdMR, MaxRsd);
269 
270  /* Flip off Boundary Conditions */
272 
273 
274  /* Now make sure we have all converged */
275  /* push(xml_out,"pade_trln_rsd");
276 write(xml_out, "RsdMR_array", RsdMR_array);
277 pop(xml_out); */
278 
279 
280 
281  for(i=0; i < Nkappa; i++) {
282  tmp = TO_REAL(1) / Kappas[i];
283  for(cb = 0; cb < Ncb; cb++) {
284  psi[i][cb] = psi[i][cb] * tmp;
285  }
286  }
287 
288  j=Nkappa;
289  for(i=0; i < Nkappa; i++) {
290  tmp = TO_REAL(1) / Kappas[i];
291  for(l=0; l < Npade; l++) {
292  for(cb = 0; cb < Ncb; cb++) {
293  psi[j][cb] = psi[j][cb] * tmp;
294  }
295  j++;
296  }
297  }
298 
299  /* Work out dumb psibar psi estimates */
300 
301 
302  for(i = 0; i < Nkappa; i++) {
303 
304  lc = trace(adj[chi[0]] * psi[i][0]);
305 
306  for(cb = 1; cb < Ncb; cb++) {
307  lct = trace(adj[chi[cb]] * psi[i][cb]);
308  lc += lct;
309  }
310 
311  dcsum = sum(lc);
312 
313  psibp_unimp[i] = FLOAT(dcsum);
314 
315  }
316 
317 
318  /* Now try to compute TrLn = b0 TrI - sum Pade_c(l) chi (M + c(l))^{-1} chi */
319 
320  /* First find the first term: Pade_b0 * Tr I, where Tr I = Nc * Nd * V */
321  TrI = TO_REAL(Ns*Nc)*TO_REAL(2)*TO_REAL(vol_cb);
322  trlns_unimp = 1;
323  tmp = TrI * Pade_b0;
326 
327  /* Now get the remaining Npade terms */
328 
329  j = Nkappa;
330  for(i = 0 ; i < Nkappa; i++) {
331  for(l = 0; l < Npade; l++) {
332 
333  /* pade_diag_term(l) = < chi, M(kappa, pade_c) chi > */
334  lc = trace(adj[chi[0]] * psi[j][0]);
335  for(cb = 1; cb < Ncb; cb++) {
336  lct = trace(adj[chi[cb]] * psi[j][cb]);
337  lc += lct;
338  }
339 
340  /* Sum up the traces */
341  dcsum = sum(lc);
342 
343  /* Change to floating point */
344  csum = FLOAT(dcsum);
345  csum = Pade_b[l] * csum;
346  pade_terms[l] = csum * minusone;
347  trlns_unimp[i] += pade_terms[l];
348  j++;
349 
350  }
351  }
352 
353  /* Work out subtraction terms. First the estimators of the
354  traceless operators */
355 
356  /* 1 Work out plaquette for 4-th order term */
358 
359  f_plaq = FLOAT(w_plaq);
360 
361 
362 
363  bilinear = 0;
364 
365 
366 
367  /* Copy chi -> r */
368  for(cb=0; cb < Ncb; cb++) {
369  r[cb] = chi[cb];
370  }
371 
372 
373  /* Flip bc-s back on */
374  phfctr (u, FORWARD);
375 
376  for(i=0; i < Nimpr; i++) {
377 
378  /* Apply Dslash: Dchi = Dslash r */
379  for(cb=0; cb < Ncb; cb++) {
380  dslash (u, r[cb], Dchi[1-cb], PLUS, cb);
381  }
382 
383  /* Now trace Dchi with chi to get the bilinear(term) chi D chi */
384  lc = trace(adj[chi[0]] * Dchi[0]);
385  for(cb = 1; cb < Ncb; cb++) {
386  lct = trace(adj[chi[cb]] * Dchi[cb]);
387  lc += lct;
388  }
389 
390  /* Sum up the traces */
391  dcsum = sum(lc);
392 
393  /* Change to floating point */
394  bilinear[i] = FLOAT(dcsum);
395 
396  /* If i corresponds to an even powered term (recall i goes from 0) */
397  if ( (i + 1) % 2 == 0 ) {
398 
399  /* Deal with the even power appropriately */
400  switch( i + 1 ) {
401  case 2:
402  /* Second order term is traceless */
403  break;
404  case 4:
405 
406  /* Here we will put the plaquette */
407  c_plaq = 1;
408  tmp = TO_REAL(64)*f_plaq*TO_REAL(3)*TO_REAL(6)*TO_REAL(Ncb)*TO_REAL(vol_cb);
409  c_plaq = tmp * c_plaq;
410  bilinear[i] += c_plaq;
411  break;
412  default:
413 
414  /* Anything higher than 4th order we will ignore for now */
415  bilinear[i] = 0;
416  break;
417  }
418 
419  }
420 
421  /* Now copy Dchi back to R for the next iteration */
422  for(cb = 0; cb < Ncb; cb++) {
423  r[cb] = Dchi[cb];
424  }
425  }
426 
427  phfctr (u, BACKWARD);
428 
429  /* We are done with Dchi -> NOw dealing with the bilinears only */
430 
431  for(kappa=0; kappa < Nkappa; kappa++) {
432  for(i = 0; i < Nimpr; i++) {
433 
434  /* psib_offterm = kappa^{i+1} bilinear(i) */
435  kappa_pow = pow((double)Kappas[kappa], (double)TO_REAL(i+1));
436 
438 
439  for(l = 0; l < Npade; l++) {
440 
441  /* Pade_offterm(order) = kappa^{i+1}/(1 + c(order))^{i+2} */
442  tmp = TO_REAL(1) / (TO_REAL(1) + Pade_c[l]);
443  pade_c_pow = pow((double)tmp, (double)TO_REAL(i+2));
445  trln_offterms[i][l][kappa] = tmp * bilinear[i];
446  }
447  }
448  }
449 
450 
452 
453  for(kappa=0; kappa < Nkappa; kappa++) {
454  for(i = 0; i < Nimpr; i++) {
455  for(l = 0; l < Npade; l++) {
457  }
458 
459  /* Multiply each term by the subtraction coefficient omega */
461 
463  }
464  }
465 
466  /* push(xml_out,"PsiBarPsiOffterms");
467 write(xml_out, "psib_offterms", psib_offterms);
468 pop(xml_out);
469  push(xml_out,"TrlnOffterms");
470 write(xml_out, "sum_trln_offterms", sum_trln_offterms);
471 pop(xml_out); */
472 
473  /* Now construct improved estimates */
476 
477  for(kappa=0; kappa < Nkappa; kappa++) {
478 
479  for(i=0; i < Nimpr; i++) {
482  }
483  }
484 
485  /* Normalise psi bar psi */
486  TrI = TO_REAL(1)/ TrI;
487 
488  for(kappa=0; kappa < Nkappa; kappa++) {
491  }
492 
493 
494 
496 
497 
499 }
Nd
Definition: meslate.cc:74
@ PLUS
Definition: chromabase.h:45
dslash(u, eta, tmp, MINUS, 1)
FloatingPoint< double > Double
Definition: gtest.h:7351
Omega_impr
Definition: pade_trln_w.cc:25
START_CODE()
int cb
Definition: pade_trln_w.cc:113
Double chi_norm
Definition: pade_trln_w.cc:85
LatticeComplex lc
Definition: pade_trln_w.cc:116
int isz
Definition: pade_trln_w.cc:151
Nkappa
Definition: pade_trln_w.cc:24
Complex csum
Definition: pade_trln_w.cc:119
Double link
Definition: pade_trln_w.cc:146
LatticeComplex lct
Definition: pade_trln_w.cc:117
MInvMRm(Munshifted, chi, psi, chi_norm, Shifts, Nshift, isz, RsdMR_array, Ncb, n_count)
int Nshift
Definition: pade_trln_w.cc:100
Pade_b
Definition: pade_trln_w.cc:25
END_CODE()
multi3d< Complex > trln_offterms(Nkappa, Npade, Nimpr)
multi1d< LatticeFermion > psi_test(Ncb)
DComplex dcsum
Definition: pade_trln_w.cc:118
Complex c_plaq
Definition: pade_trln_w.cc:148
Real TrI
Definition: pade_trln_w.cc:127
RsdMR
Definition: pade_trln_w.cc:26
multi1d< LatticeFermion > r(Ncb)
chi
Definition: pade_trln_w.cc:24
Real MaxRsd
Definition: pade_trln_w.cc:107
bilinear
Definition: pade_trln_w.cc:363
Real minshift
Definition: pade_trln_w.cc:152
Real Shiftmass
Definition: pade_trln_w.cc:105
Double t_plaq
Definition: pade_trln_w.cc:145
LINEAR_OPERATOR(D)
u
Definition: pade_trln_w.cc:24
multi2d< Complex > psib_offterms(Nkappa, Nimpr)
Kappas
Definition: pade_trln_w.cc:24
psi
Definition: pade_trln_w.cc:191
int n_count
Definition: pade_trln_w.cc:69
PRINTF("TrLnSolve: n_count = %d, RsdMR = %16.8e, MaxRsd = %16.8e\n", n_count, RsdMR, MaxRsd)
multi1d< Real > RsdMR_array(Nshift)
trlns_imp
Definition: pade_trln_w.cc:26
FILL(RsdMR_array, RsdMR)
Real pade_c_pow
Definition: pade_trln_w.cc:131
psibp_imp
Definition: pade_trln_w.cc:26
multi1d< Real > Shifts(Nshift)
int l
Definition: pade_trln_w.cc:111
Nimpr
Definition: pade_trln_w.cc:25
Double w_plaq
Definition: pade_trln_w.cc:143
int kappa
Definition: pade_trln_w.cc:112
CONSTRUCT_LINEAR_OPERATOR(Munshifted, ldunshifted, D)
include(types.mh) SUBROUTINE(PadeTrLn
trlns_unimp
Definition: pade_trln_w.cc:26
Double s_plaq
Definition: pade_trln_w.cc:144
Real f_plaq
Definition: pade_trln_w.cc:147
Npade
Definition: pade_trln_w.cc:24
DestDslash(D, WILSON_DSLASH)
multi1d< Real > first_terms(Nkappa)
Real tmp
Definition: pade_trln_w.cc:129
psibp_unimp
Definition: pade_trln_w.cc:26
Double r_norm
Definition: pade_trln_w.cc:86
sum_trln_offterms
Definition: pade_trln_w.cc:451
int j
Definition: pade_trln_w.cc:110
FREE_LINEAR_OPERATOR(Munshifted)
MesPlq(u, w_plaq, s_plaq, t_plaq, link)
for(i=0;i< Nkappa;i++)
Definition: pade_trln_w.cc:167
int i
Definition: pade_trln_w.cc:109
Ncb
Definition: pade_trln_w.cc:24
Real minusone
Definition: pade_trln_w.cc:128
Real kappa_pow
Definition: pade_trln_w.cc:130
phfctr(u, FORWARD)
Pade_b0
Definition: pade_trln_w.cc:25
ConsDslash(D, u, Shiftmass, WILSON_DSLASH)
multi1d< Complex > pade_terms(Npade)
Pade_c
Definition: pade_trln_w.cc:25
multi1d< LatticeFermion > Dchi(Ncb)
#define FORWARD
Definition: primitives.h:82
#define WILSON_DSLASH
Definition: primitives.h:50
#define BACKWARD
Definition: primitives.h:83
Double sum
Definition: qtopcor.cc:37