CHROMA
mciter32.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief One heatbath interation of updating the gauge field configuration
3  */
4 
5 //#error "Not tested (or even compiled). However, reasonably well converted."
6 
7 //#error "THIS VERSION SHOULD BE USED IN PLACE OF MCITER.CC - IT IS MORE GENERAL. REMOVE THIS FILE"
8 
9 //#warning "THERE ARE INSTANCES OF shift(...,mu)*shift(...,mu) THAT SHOULD BE OPTIMIZED"
10 
11 #include "chromabase.h"
12 #include "util/gauge/reunit.h"
13 #include "update/heatbath/su3over.h"
14 #include "update/heatbath/su3hb.h"
15 #include "update/heatbath/su2_hb_update.h"
16 
17 #if 1
18 #warning "FIXME"
19 const int igluetmp = 1;
20 const Real GlueCoeffRT(0.0);
21 const Real GlueCoeffPG(0.0);
22 #endif
23 
24 
25 namespace Chroma {
26 
27 
28 //! One heatbath interation of updating the gauge field configuration
29 /*!
30  * \ingroup heatbath
31  *
32  * Make one interation of updating the gauge field configuration:
33  * for Wilson or Symmanzik improved pure gauge action:
34  * this consists of n_over overrelaxation sweeps followed
35  * by one heatbath sweep with nheat trials.
36  * In the case of SU(3), for each link we loop over the 3 SU(2) subgroups.
37  *
38  * Because of the option for the Symanzik improved action we need 2^(d+1)
39  * sublattices, deviding first into 2^d hypercubes and then
40  * checkerboarding those.
41  *
42  *
43  * Warning: this works only for Nc = 2 and 3 !
44  *
45  * \param u gauge field ( Modify )
46  * \param n_over number of overrelaxation sweeps ( Read )
47  * \param nheat number of heatbath trials ( Read )
48  * \param NmaxHB max number of heatbath hits ( Read )
49  * \param ntrials total number of individual heatbath trials ( Modify )
50  * \param nfails total number of individual heatbath failures ( Modify )
51  */
52 
53 void mciter32(multi1d<LatticeColorMatrix>& u,
54  int n_over, int nheat, int NmaxHB,
55  int& ntrials, int& nfails,
56  const Set& ss,
57  const multi3d<int>& neighsubl)
58 {
59  START_CODE();
60 
61  // Expect neighsubl to be declared multi3d<int> neighsubl(Nd,2,ss.numSubsets());
62 
63  LatticeColorMatrix tmp_1;
64  LatticeColorMatrix tmp_2;
65  LatticeColorMatrix u_staple;
66  LatticeColorMatrix u_staple2;
67 
68 
69 #if 0
70  igluetmp = fabs(GlueImp);
71  if (igluetmp > 2)
72  QDP_error_exit("Illegal value of GlueImp", GlueImp);
73 #endif
74 
75 
76  for(int iter = 0; iter <= n_over; ++iter)
77  {
78  for(int subl = 0; subl < ss.numSubsets(); ++subl)
79  for(int mu = 0; mu < Nd; ++mu)
80  {
81  u_staple[ss[subl]] = 0;
82 
83  int fsublm = neighsubl[mu][1][subl];
84  int bsublm = neighsubl[mu][0][subl];
85 
86  for(int nu = 0; nu < Nd; ++nu)
87  {
88  if(nu == mu) continue;
89 
90  int fsubln = neighsubl[nu][1][subl];
91  int bsubln = neighsubl[nu][0][subl];
92 
93  /* Forward staple */
94  /* tmp_1(x) = u(x+mu,nu)*u_dag(x+nu,mu) */
95  tmp_1[ss[subl]] = shift(u[nu], FORWARD, mu) * shift(adj(u[mu]), FORWARD, nu);
96 
97  /* u_staple(x) += tmp_1 * u_dag(x,nu)
98  += u(x+mu,nu)*u_dag(x+nu,mu)*u_dag(x,nu) */
99  u_staple[ss[subl]] += tmp_1 * adj(u[nu]);
100 
101  /* Backward staple */
102  int bfsubl = neighsubl[mu][1][bsubln];
103 
104  /* tmp_1(x) = u(x,mu)*u(x+mu,nu) */
105  tmp_1[ss[bsubln]] = u[mu] * shift(u[nu], FORWARD, mu);
106 
107  /* u_staple(x) += tmp_1_dag(x-nu) * u(x-nu,nu)
108  += u_dag(x+mu-nu,nu)*u_dag(x-nu,mu)*u(x-nu,nu) */
109  u_staple[ss[subl]] += shift(adj(tmp_1), BACKWARD, nu) * shift(u[nu], BACKWARD, nu);
110  }
111 
112  /* Symmanzik improvement terms */
113  if (igluetmp >= 1)
114  {
115  u_staple2 = 0;
116 
117  for(int nu = 0; nu < Nd; ++nu)
118  {
119  if(nu == mu) continue;
120 
121  int fsubln = neighsubl[nu][1][subl];
122  int bsubln = neighsubl[nu][0][subl];
123 
124  /* forward 1x2 rectangle in mu,nu directions, nu!=mu */
125  int fsublmn = neighsubl[mu][1][fsubln];
126  int fsubl2 = neighsubl[nu][1][fsubln];
127 
128  /* tmp_1(x) = u(x+mu,nu)*u_dag(x+nu,mu) */
129  tmp_1[ss[fsubln]] = shift(u[nu], FORWARD, mu) * shift(adj(u[mu]), FORWARD, nu);
130 
131  /* tmp_2(x) = tmp_1(x+nu)*u_dag(x,nu)
132  = u(x+mu,nu)*u_dag(x+nu,mu)*u_dag(x,nu) */
133  tmp_2[ss[fsubln]] = tmp_1 * adj(u[nu]);
134 
135  /* tmp_1(x) = tmp_2(x+nu)*u_dag(x,nu)
136  = u(x+nu+mu,nu)*u_dag(x+2nu,mu)*
137  u_dag(x+nu,nu)*u_dag(x,nu) */
138  tmp_1[ss[subl]] = shift(tmp_2, FORWARD, nu) * adj(u[nu]);
139 
140  /* u_staple2 += u(x+mu,nu)*tmp_1(x) */
141  u_staple2[ss[subl]] += shift(u[nu], FORWARD, mu) * tmp_1;
142 
143  /* backward 1x2 rectangle in mu,nu directions, nu!=mu */
144  int bsubl2 = neighsubl[nu][0][bsubln];
145  int bfsubl = neighsubl[mu][1][bsubl2];
146 
147  /* tmp_1(x) = u(x,mu)*u(x+mu,nu) */
148  tmp_1[ss[bsubl2]] = u[mu] * shift(u[nu], FORWARD, mu);
149 
150  /* tmp_2(x) = tmp_1_dag(x-nu) * u(x-nu,nu)
151  = u_dag(x+mu-nu,nu)*u_dag(x-nu,mu)*u(x-nu,nu) */
152  tmp_2[ss[bsubln]] = shift(adj(tmp_1), BACKWARD, nu) * shift(u[nu], BACKWARD, nu);
153 
154  /* tmp_1(x) = u_dag(x+mu,nu)*tmp_2(x)
155  = u_dag(x+mu,nu)*u_dag(x+mu-nu,nu)*
156  u_dag(x-nu,mu)*u(x-nu,nu) */
157  bfsubl = neighsubl[mu][1][bsubln];
158  tmp_1[ss[bsubln]] = shift(adj(u[nu]), FORWARD, mu) * tmp_2;
159 
160  /* u_staple2(x) += tmp_1(x-nu) * u(x-nu,nu) */
161  u_staple2[ss[subl]] += shift(tmp_1, BACKWARD, nu) * shift(u[nu], BACKWARD, nu);
162 
163  /* forward 2x1 rectangle in mu,nu directions open at bottom left */
164  fsubl2 = neighsubl[mu][1][fsublm];
165  fsublmn = neighsubl[nu][1][fsublm];
166 
167  /* tmp_1(x) = u(x,mu)*u(x+mu,nu) */
168  tmp_1[ss[fsublm]] = u[mu] * shift(u[nu], FORWARD, mu);
169 
170  /* tmp_2(x) = tmp_1(x)*u_dag(x+nu,mu)
171  = u(x,mu)*u(x+mu,nu)*u_dag(x+nu,mu) */
172  tmp_2[ss[fsublm]] = tmp_1 * shift(adj(u[mu]), FORWARD, nu);
173 
174  /* tmp_1(x) = tmp_2(x+mu)*u_dag(x+nu,mu)
175  = u(x+mu,mu)*u(x+2mu,nu)*
176  u_dag(x+mu+nu,mu)*u_dag(x+nu,mu) */
177  tmp_1[ss[subl]] = shift(tmp_2, FORWARD, mu) * shift(adj(u[mu]), FORWARD, nu);
178 
179  /* u_staple2(x) += tmp_1 * u_dag(x,nu) */
180  u_staple2[ss[subl]] += tmp_1 * adj(u[nu]);
181 
182  /* backward 2x1 rectangle in mu,nu directions open at top left */
183  bfsubl = neighsubl[mu][1][bsubln];
184  fsubl2 = neighsubl[mu][1][bfsubl];
185 
186  /* tmp_1(x) = u(x,mu)*u(x+mu,nu) */
187  tmp_1[ss[bfsubl]] = u[mu] * shift(u[nu], FORWARD, mu);
188 
189  /* tmp_2(x) = u(x,mu)*tmp_1(x+mu)
190  = u(x,mu)*u(x+mu,mu)*u(x+2mu,nu) */
191  tmp_2[ss[bsubln]] = u[mu] * shift(tmp_1, FORWARD, mu);
192 
193  /* tmp_1(x) = tmp_2_dag(x-nu)*u(x-nu,nu)
194  = u_dag(x-nu+2mu,nu)*u_dag(x-nu+mu,mu)*
195  u_dag(x-nu,mu)*u(x-nu,nu) */
196  tmp_1[ss[subl]] = shift(adj(tmp_2), BACKWARD, nu) * shift(u[nu], BACKWARD, nu);
197 
198  /* u_staple2(x) += u(x+mu,mu) * tmp_1 */
199  u_staple2[ss[subl]] += shift(u[mu], FORWARD, mu) * tmp_1;
200 
201  /* forward 2x1 rectangle in mu,nu directions open at bottom right */
202  /* tmp_1(x) = u(x,nu)*u(x+nu,mu) */
203  bfsubl = neighsubl[nu][1][bsublm];
204  tmp_1[ss[bsublm]] = u[nu] * shift(u[mu], FORWARD, nu);
205 
206  /* tmp_2(x) = tmp_1_dag(x-mu)*u(x-mu,mu)
207  = u_dag(x-mu+nu,mu)*u_dag(x-mu,nu)*u(x-mu,mu) */
208  tmp_2[ss[subl]] = shift(adj(tmp_1), BACKWARD, mu) * shift(u[mu], BACKWARD, mu);
209 
210  /* tmp_1(x) = u_dag(x+nu,mu)*tmp_2(x)
211  = u_dag(x+nu,mu)*u_dag(x-mu+nu,mu)*
212  u_dag(x-mu,nu)*u(x-mu,mu) */
213  tmp_1[ss[subl]] = shift(adj(u[mu]), FORWARD, nu) * tmp_2;
214 
215  /* u_staple2(x) += u(x+mu,nu) * tmp_1 */
216  u_staple2[ss[subl]] += shift(u[nu], FORWARD, mu) * tmp_1;
217 
218  /* backward 2x1 rectangle in mu,nu directions open at top right */
219  bfsubl = neighsubl[mu][1][bsubln];
220  int bsublmn = neighsubl[mu][0][bsubln];
221 
222  /* tmp_1(x) = u_dag(x,mu)*u(x,nu) */
223  tmp_1[ss[bsublmn]] = adj(u[mu]) * u[nu];
224 
225  /* tmp_2(x) = u_dag(x,mu)*tmp_1(x-mu)
226  = u_dag(x,mu)*u_dag(x-mu,mu)*u(x-mu,nu) */
227  tmp_2[ss[bsubln]] = adj(u[mu]) * shift(tmp_1, BACKWARD, mu);
228 
229  /* tmp_1(x) = u_dag(x+mu,nu)*tmp_2(x)
230  = u_dag(x+mu,nu)*u_dag(x,mu)*
231  u_dag(x-mu,mu)*u(x-mu,nu) */
232  tmp_1[ss[bsubln]] = shift(adj(u[nu]), FORWARD, mu) * tmp_2;
233 
234  /* u_staple2(x) += tmp_1(x-nu) * u(x-mu,mu) */
235  u_staple2[ss[subl]] += shift(tmp_1, BACKWARD, nu) * shift(u[mu], BACKWARD, mu);
236  }
237 
238  u_staple[ss[subl]] += u_staple2 * GlueCoeffRT;
239  }
240 
241  /* Symmanzik 6-link parallelograms */
242  if (igluetmp == 2)
243  {
244  u_staple2 = 0;
245 
246  /* In each cube, labeled by eta<nu<rho there are 4 ways
247  to get a 6-link parallelogram: (eta,nu,rho,-eta,-nu,-rho),
248  (nu,eta,rho,-nu,-eta,-rho), (eta,rho,nu,-eta,-rho,-nu)
249  and (eta,-nu,rho,-eta,nu,-rho). The first three can be
250  obtained by keeping only (eta,nu,rho,-eta,-nu,-rho)
251  while restricting only to eta<rho, nu!=eta, nu!=rho. */
252  for(int eta = 0; eta < Nd-1; ++eta)
253  for(int nu = 0; nu < Nd; ++nu)
254  {
255  if(nu == eta) continue;
256 
257  int fsuble = neighsubl[eta][1][subl];
258  int bsuble = neighsubl[eta][0][subl];
259  int fsubln = neighsubl[nu][1][subl];
260  int bsubln = neighsubl[nu][0][subl];
261 
262  for(int rho = eta+1; rho < Nd; ++rho)
263  {
264  if(rho == nu) continue;
265 
266  int fsublr = neighsubl[rho][1][subl];
267  int bsublr = neighsubl[rho][0][subl];
268 
269  if(mu == eta)
270  {
271  /* "update" link is first link */
272  /* tmp_1(x) = u(x,nu)*u(x+nu,mu) */
273  int fsublmn = neighsubl[nu][1][fsublr];
274  tmp_1[ss[fsublr]] = u[nu] * shift(u[mu], FORWARD, nu);
275 
276  /* tmp_2(x) = u(x,rho)*tmp_1(x+rho)
277  = u(x,rho)*u(x+rho,nu)*u(x+rho+nu,mu) */
278  tmp_2[ss[subl]] = u[rho] * shift(tmp_1, FORWARD, rho);
279 
280  /* tmp_1(x) = u(x,nu)*u(x+nu,rho) */
281  fsublmn = neighsubl[nu][1][fsublm];
282  tmp_1[ss[fsublm]] = u[nu] * shift(u[rho], FORWARD, nu);
283 
284  /* u_staple2(x) += tmp_1(x+mu) * tmp_2_dag(x) */
285  u_staple2[ss[subl]] += shift(tmp_1, FORWARD, mu) * adj(tmp_2);
286 
287  /* "update" link is fourth link */
288  int bsubl2 = neighsubl[rho][0][bsubln];
289  int bfsubl = neighsubl[mu][1][bsubl2];
290  int fsubl2 = neighsubl[nu][1][bfsubl];
291 
292  /* tmp_1(x) = u(x,nu)*u(x+nu,rho) */
293  tmp_1[ss[bfsubl]] = u[nu] * shift(u[rho], FORWARD, nu);
294 
295  /* tmp_2(x) = u(x,mu)*tmp_1(x+mu)
296  = u(x,mu)*u(x+mu,nu)*u(x+mu+nu,rho) */
297  tmp_2[ss[bsubl2]] = u[mu] * shift(tmp_1, FORWARD, mu);
298 
299  /* tmp_1(x) = tmp_2_dag(x-rho)*u(x-rho,rho)
300  = u_dag(x-rho+mu+nu,rho)*u_dag(x-rho+mu,nu)*
301  u_dag(x-rho,mu)*u(x-rho,rho) */
302  tmp_1[ss[bsubln]] = shift(adj(tmp_2), BACKWARD, rho) * shift(u[rho], BACKWARD, rho);
303 
304  /* u_staple2(x) += tmp_1(x-nu) * u(x-nu,nu) */
305  u_staple2[ss[subl]] += shift(tmp_1, BACKWARD, nu) * shift(u[nu], BACKWARD, nu);
306 
307  }
308  else if(mu == nu)
309  {
310  /* "update" link is second link */
311  int bfsubl = neighsubl[rho][1][bsuble];
312  int fsubl2 = neighsubl[mu][1][bfsubl];
313 
314  /* tmp_1(x) = u(x,mu)*u(x+mu,eta) */
315  tmp_1[ss[bfsubl]] = u[mu] * shift(u[eta], FORWARD, mu);
316 
317  /* tmp_2(x) = u(x,rho)*tmp_1(x+rho)
318  = u(x,rho)*u(x+rho,mu)*u(x+rho+mu,eta) */
319  tmp_2[ss[bsuble]] = u[rho] * shift(tmp_1, FORWARD, rho);
320 
321  /* tmp_1(x) = tmp_2_dag(x)*u(x,eta)
322  = u_dag(x+rho+mu,eta)*u_dag(x+rho,mu)*
323  u_dag(x,rho)*u(x,eta) */
324  tmp_1[ss[bsuble]] = adj(tmp_2) * u[eta];
325 
326  /* u_staple2(x) += u(x+mu,rho) * tmp_1(x-eta) */
327  u_staple2[ss[subl]] += shift(u[rho], FORWARD, mu) * shift(tmp_1, BACKWARD, eta);
328 
329  /* "update" link is fifth link */
330  bfsubl = neighsubl[eta][1][bsublr];
331  fsubl2 = neighsubl[mu][1][bfsubl];
332 
333  /* tmp_1(x) = u(x,mu)*u(x+mu,rho) */
334  tmp_1[ss[bfsubl]] = u[mu] * shift(u[rho], FORWARD, mu);
335 
336  /* tmp_2(x) = u(x,eta)*tmp_1(x+eta)
337  = u(x,eta)*u(x+eta,mu)*u(x+eta+mu,rho) */
338  tmp_2[ss[bsublr]] = u[eta] * shift(tmp_1, FORWARD, eta);
339 
340  /* tmp_1(x) = tmp_2_dag(x)*u(x,rho)
341  = u_dag(x+eta+mu,eta)*u_dag(x+eta,mu)*
342  u_dag(x,eta)*u(x,rho) */
343  tmp_1[ss[bsublr]] = adj(tmp_2) * u[rho];
344 
345  /* u_staple2(x) += u(x+mu,eta) * tmp_1(x-rho) */
346  u_staple2[ss[subl]] += shift(u[eta], FORWARD, mu) * shift(tmp_1, BACKWARD, rho);
347 
348  }
349  else if(mu == rho)
350  {
351  /* "update" link is third link */
352  int bsubl2 = neighsubl[eta][0][bsubln];
353  int bfsubl = neighsubl[mu][1][bsubl2];
354  int fsubl2 = neighsubl[nu][1][bfsubl];
355 
356  /* tmp_1(x) = u(x,nu)*u(x+nu,eta) */
357  tmp_1[ss[bfsubl]] = u[nu] * shift(u[eta], FORWARD, nu);
358 
359  /* tmp_2(x) = u(x,mu)*tmp_1(x+mu)
360  = u(x,mu)*u(x+mu,nu)*u(x+mu+nu,eta) */
361  tmp_2[ss[bsubl2]] = u[mu] * shift(tmp_1, FORWARD, mu);
362 
363  /* tmp_1(x) = tmp_2_dag(x-eta)*u(x-eta,eta)
364  = u_dag(x-eta+mu+nu,eta)*u_dag(x-eta+mu,nu)*
365  u_dag(x-eta,mu)*u(x-eta,eta) */
366  tmp_1[ss[bsubln]] = shift(adj(tmp_2), BACKWARD, eta) * shift(u[eta], BACKWARD, eta);
367 
368  /* u_staple2(x) += tmp_1(x-nu) * u(x-nu,nu) */
369  u_staple2[ss[subl]] += shift(tmp_1, BACKWARD, nu) * shift(u[nu], BACKWARD, nu);
370 
371  /* "update" link is sixth link */
372  /* tmp_1(x) = u(x,nu)*u(x+nu,mu) */
373  int fsublmn = neighsubl[nu][1][fsuble];
374  tmp_1[ss[fsuble]] = u[nu] * shift(u[mu], FORWARD, nu);
375 
376  /* tmp_2(x) = u(x,eta)*tmp_1(x+eta)
377  = u(x,eta)*u(x+eta,nu)*u(x+eta+nu,mu) */
378  tmp_2[ss[subl]] = u[eta] * shift(tmp_1, FORWARD, eta);
379 
380  /* tmp_1(x) = u(x,nu)*u(x+nu,eta) */
381  fsublmn = neighsubl[nu][1][fsublm];
382  tmp_1[ss[fsublm]] = u[nu] * shift(u[eta], FORWARD, nu);
383 
384  /* u_staple2(x) += tmp_1(x+mu) * tmp_2_dag(x) */
385  u_staple2[ss[subl]] += shift(tmp_1, FORWARD, mu) * adj(tmp_2);
386  }
387  }
388  }
389 
390  /* Finally the (eta,-nu,rho,-eta,nu,-rho) parallelogram */
391  for(int eta = 0; eta < Nd-2; ++eta)
392  for(int nu = eta+1; nu < Nd-1; ++nu)
393  {
394  int fsuble = neighsubl[eta][1][subl];
395  int bsuble = neighsubl[eta][0][subl];
396  int fsubln = neighsubl[nu][1][subl];
397  int bsubln = neighsubl[nu][0][subl];
398 
399  for(int rho = nu+1; rho < Nd; ++rho)
400  {
401  int fsublr = neighsubl[rho][1][subl];
402  int bsublr = neighsubl[rho][0][subl];
403 
404  if(mu == eta)
405  {
406  /* "update" link is first link */
407  /* tmp_1(x) = u_dag(x-nu,mu)*u(x-nu,nu) */
408  int bfsubl = neighsubl[nu][0][fsublr];
409  tmp_1[ss[fsublr]] = shift(adj(u[mu]), BACKWARD, nu) * shift(u[nu], BACKWARD, nu);
410 
411  /* tmp_2(x) = tmp_1(x+rho)*u_dag(x,rho)
412  = u_dag(x+rho-nu,mu)*u(x+rho-nu,nu)*u_dag(x,rho) */
413  tmp_2[ss[subl]] = shift(tmp_1, FORWARD, rho) * adj(u[rho]);
414 
415  /* tmp_1(x) = u_dag(x-nu,nu)*u(x-nu,rho) */
416  bfsubl = neighsubl[nu][0][fsublm];
417  tmp_1[ss[fsublm]] = shift(adj(u[nu]), BACKWARD, nu) * shift(u[rho], BACKWARD, nu);
418 
419  /* u_staple2(x) += tmp_1(x+mu) * tmp_2(x) */
420  u_staple2[ss[subl]] += shift(tmp_1, FORWARD, mu) * tmp_2;
421 
422  /* "update" link is forth link */
423  /* tmp_1(x) = u_dag(x-rho,mu)*u(x-rho,rho) */
424  bfsubl = neighsubl[rho][0][fsubln];
425  tmp_1[ss[fsubln]] = shift(adj(u[mu]), BACKWARD, rho) * shift(u[rho], BACKWARD, rho);
426 
427  /* tmp_2(x) = tmp_1(x+nu)*u_dag(x,nu)
428  = u_dag(x+nu-rho,mu)*u(x+nu-rho,rho)*u_dag(x,nu) */
429  tmp_2[ss[subl]] = shift(tmp_1, FORWARD, nu) * adj(u[nu]);
430 
431  /* tmp_1(x) = u_dag(x-rho,rho)*u(x-rho,nu) */
432  bfsubl = neighsubl[rho][0][fsublm];
433  tmp_1[ss[fsublm]] = shift(adj(u[rho]), BACKWARD, rho) * shift(u[nu], BACKWARD, rho);
434 
435  /* u_staple2(x) += tmp_1(x+mu) * tmp_2(x) */
436  u_staple2[ss[subl]] += shift(tmp_1, FORWARD, mu) * tmp_2;
437 
438  }
439  else if(mu == nu)
440  {
441  /* "update" link is second link */
442  /* tmp_1(x) = u_dag(x-eta,mu)*u(x-eta,eta) */
443  int bfsubl = neighsubl[eta][0][fsublr];
444  tmp_1[ss[fsublr]] = shift(adj(u[mu]), BACKWARD, eta) * shift(u[eta], BACKWARD, eta);
445 
446  /* tmp_2(x) = tmp_1(x+rho)*u_dag(x,rho)
447  = u_dag(x+rho-eta,mu)*u(x+rho-eta,eta)*
448  u_dag(x,rho) */
449  tmp_2[ss[subl]] = shift(tmp_1, FORWARD, rho) * adj(u[rho]);
450 
451  /* tmp_1(x) = u_dag(x-eta,eta)*u(x-eta,rho) */
452  bfsubl = neighsubl[eta][0][fsublm];
453  tmp_1[ss[fsublm]] = shift(adj(u[eta]), BACKWARD, eta) * shift(u[rho], BACKWARD, eta);
454 
455  /* u_staple2(x) += tmp_1(x+mu) * tmp_2(x) */
456  u_staple2[ss[subl]] += shift(tmp_1, FORWARD, mu) * tmp_2;
457 
458  /* "update" link is fifth link */
459  /* tmp_1(x) = u_dag(x-rho,mu)*u(x-rho,rho) */
460  bfsubl = neighsubl[rho][0][fsuble];
461  tmp_1[ss[fsuble]] = shift(adj(u[mu]), BACKWARD, rho) * shift(u[rho], BACKWARD, rho);
462 
463  /* tmp_2(x) = tmp_1(x+eta)*u_dag(x,eta)
464  = u_dag(x+eta-rho,mu)*u(x+eta-rho,rho)*
465  u_dag(x,eta) */
466  tmp_2[ss[subl]] = shift(tmp_1, FORWARD, eta) * adj(u[eta]);
467 
468  /* tmp_1(x) = u_dag(x-rho,rho)*u(x-rho,eta) */
469  bfsubl = neighsubl[rho][0][fsublm];
470  tmp_1[ss[fsublm]] = shift(adj(u[rho]), BACKWARD, rho) * shift(u[eta], BACKWARD, rho);
471 
472  /* u_staple2(x) += tmp_1(x+mu) * tmp_2(x) */
473  u_staple2[ss[subl]] += shift(tmp_1, FORWARD, mu) * tmp_2;
474 
475  }
476  else if(mu == rho)
477  {
478  /* "update" link is third link */
479  /* tmp_1(x) = u_dag(x-eta,mu)*u(x-eta,eta) */
480  int bfsubl = neighsubl[eta][0][fsubln];
481  tmp_1[ss[fsubln]] = shift(adj(u[mu]), BACKWARD, eta) * shift(u[eta], BACKWARD, eta);
482 
483  /* tmp_2(x) = tmp_1(x+nu)*u_dag(x,nu)
484  = u_dag(x+nu-eta,mu)*u(x+nu-eta,eta)*
485  u_dag(x,nu) */
486  tmp_2[ss[subl]] = shift(tmp_1, FORWARD, nu) * adj(u[nu]);
487 
488  /* tmp_1(x) = u_dag(x-eta,eta)*u(x-eta,nu) */
489  bfsubl = neighsubl[eta][0][fsublm];
490  tmp_1[ss[fsublm]] = shift(adj(u[eta]), BACKWARD, eta) * shift(u[nu], BACKWARD, eta);
491 
492  /* u_staple2(x) += tmp_1(x+mu) * tmp_2(x) */
493  u_staple2[ss[subl]] += shift(tmp_1, FORWARD, mu) * tmp_2;
494 
495  /* "update" link is sixth link */
496  /* tmp_1(x) = u_dag(x-nu,mu)*u(x-nu,nu) */
497  bfsubl = neighsubl[nu][0][fsuble];
498  tmp_1[ss[fsuble]] = shift(adj(u[mu]), BACKWARD, nu) * shift(u[nu], BACKWARD, nu);
499 
500  /* tmp_2(x) = tmp_1(x+eta)*u_dag(x,eta)
501  = u_dag(x+eta-nu,mu)*u(x+eta-nu,nu)*
502  u_dag(x,eta) */
503  tmp_2[ss[subl]] = shift(tmp_1, FORWARD, eta) * adj(u[eta]);
504 
505  /* tmp_1(x) = u_dag(x-nu,nu)*u(x-nu,eta) */
506  bfsubl = neighsubl[nu][0][fsublm];
507  tmp_1[ss[fsublm]] = shift(adj(u[nu]), BACKWARD, nu) * shift(u[eta], BACKWARD, nu);
508 
509  /* u_staple2(x) += tmp_1(x+mu) * tmp_2(x) */
510  u_staple2[ss[subl]] += shift(tmp_1, FORWARD, mu) * tmp_2;
511  }
512  }
513  }
514 
515  u_staple[ss[subl]] += u_staple2 * GlueCoeffPG;
516  }
517 
518 
519  if ( iter < n_over )
520  {
521  /* Do an overrelaxation step */
522  /*# Loop over SU(2) subgroup index */
523  for(int su2_index = 0; su2_index < Nc*(Nc-1)/2; ++su2_index)
524  su3over(u[mu], u_staple, su2_index, ss[subl]);
525  }
526  else
527  {
528  /* Do a heatbath step */
529  /*# Loop over SU(2) subgroup index */
530  for(int su2_index = 0; su2_index < Nc*(Nc-1)/2; ++su2_index)
531  {
532  int ntry;
533  int nfail;
534 
535  su2_hb_update(u[mu], u_staple,
536  Real(2.0/Nc),
537  su2_index, ss[subl],
538  NmaxHB);
539 
540 // su3hb(u[mu], u_staple, su2_index, nheat, ntry, nfail, ss[subl]);
541  ntrials = ntrials + ntry;
542  nfails = nfails + nfail;
543  }
544 
545 
546  /* Reunitarize */
547  // We seem to have different reunitarizes here. Would be better to have a subset version.
548 /* reunit (u[mu][subl], lbad, OPTION[REUNITARIZE_ERROR], numbad); */
549 // reunit (u[mu][subl]);
550  reunit (u[mu]);
551 
552  }
553 
554 #if 0
555  /* If using Schroedinger functional, reset the boundaries */
556  if ( SchrFun > 0 )
557  {
558  copymask(u[mu][subl], lSFmask[mu][subl], SFBndFld[mu][subl], REPLACE);
559  }
560 #endif
561  }
562  }
563 
564  END_CODE();
565 }
566 
567 } // end namespace Chroma
chromabase.h
Primary include file for CHROMA library code.
GlueImp
int GlueImp
mu
int mu
Definition: cool.cc:24
su2_index
int su2_index
Definition: cool.cc:27
nu
int nu
Definition: cool.cc:25
Chroma::su3over
void su3over(LatticeColorMatrix &u, const LatticeColorMatrix &w, int su2_index, const Subset &sub)
Do one SU(2) subgroup microcanonical overrelaxation update of SU(Nc) matrix.
Definition: su3over.cc:26
Chroma::mciter32
void mciter32(multi1d< LatticeColorMatrix > &u, int n_over, int nheat, int NmaxHB, int &ntrials, int &nfails, const Set &ss, const multi3d< int > &neighsubl)
One heatbath interation of updating the gauge field configuration.
Definition: mciter32.cc:53
igluetmp
const int igluetmp
Definition: mciter32.cc:19
GlueCoeffRT
const Real GlueCoeffRT(0.0)
GlueCoeffPG
const Real GlueCoeffPG(0.0)
tmp_2
LatticeColorMatrix tmp_2
Definition: meslate.cc:51
tmp_1
LatticeColorMatrix tmp_1
Definition: meslate.cc:50
Nd
Nd
Definition: meslate.cc:74
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::QDP_error_exit
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)
Chroma::u
static multi1d< LatticeColorMatrix > u
Definition: syssolver_linop_qop_mg_w.cc:39
Chroma::reunit
void reunit(LatticeColorMatrixF3 &xa)
Definition: reunit.cc:467
Chroma::eta
LatticeFermion eta
Definition: mespbg5p_w.cc:37
Chroma::START_CODE
START_CODE()
Chroma::su2_hb_update
void su2_hb_update(LatticeColorMatrix &u_mu, const LatticeColorMatrix &u_mu_staple, Double BetaMC, const int su2_index, const Subset &sub, const int NmaxHB)
Definition: su2_hb_update.cc:17
FORWARD
#define FORWARD
Definition: primitives.h:82
BACKWARD
#define BACKWARD
Definition: primitives.h:83
copymask
copymask(lcoord, lbit, ltmp_1, REPLACE)
reunit.h
Reunitarize in place a color matrix to SU(N)
su2_hb_update.h
su3hb.h
Do one SU(2) subgroup heatbath update of SU(Nc) matrix U with action.
su3over.h
Do one SU(2) subgroup microcanonical overrelaxation update of SU(Nc)
Generated by doxygen 1.9.1