Skip to content

Commit fc13eb0

Browse files
jing2lijing2li
committed
defect added
1 parent fd7f785 commit fc13eb0

5 files changed

Lines changed: 96 additions & 17 deletions

File tree

‎src/lib/get_rectangle_staples.c‎

Lines changed: 25 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@
2525
#include "get_rectangle_staples.h"
2626
#include "global.h"
2727
#include "su3.h"
28+
#include "ptbc.h"
2829

2930
void get_rectangle_staples(su3 *const v, const int x, const int mu) {
3031
get_rectangle_staples_general(v, x, mu, (const su3 *const *const)g_gauge_field);
@@ -51,14 +52,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
5152
_su3_times_su3(tmp1, *a, *b);
5253
z = g_iup[y][nu];
5354
c = &gf[z][mu];
55+
double const ptbc_fac0 = get_ptbc_coeff(x, nu) * get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, mu);
5456
_su3_times_su3(tmp2, tmp1, *c);
5557

5658
y = g_iup[x][mu];
5759
d = &gf[y][nu];
5860
z = g_iup[y][nu];
5961
e = &gf[z][nu];
62+
double const ptbc_fac1 = get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, nu);
6063
_su3_times_su3(tmp1, *d, *e);
61-
_su3_times_su3d_acc((*v), tmp2, tmp1);
64+
//_su3_times_su3d_acc((*v), tmp2, tmp1);
65+
_real_times_su3_times_su3d_acc((*v), tmp2, tmp1, ptbc_fac0 * ptbc_fac1);
6266

6367
/* 1 contr. starting idn[idn[x][nu]][nu]
6468
* e^+ d^+ a b c
@@ -73,14 +77,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
7377
a = &gf[z][mu];
7478
_su3d_times_su3(tmp1, *d, *a);
7579
e = &gf[y][nu];
80+
double const ptbc_fac2 = get_ptbc_coeff(z, nu) * get_ptbc_coeff(z, mu) * get_ptbc_coeff(y, nu);
7681
_su3d_times_su3(tmp2, *e, tmp1);
7782

7883
y = g_iup[z][mu];
7984
b = &gf[y][nu];
8085
z = g_iup[y][nu];
8186
c = &gf[z][nu];
87+
double const ptbc_fac3 = get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, nu);
8288
_su3_times_su3(tmp1, *b, *c);
83-
_su3_times_su3_acc((*v), tmp2, tmp1);
89+
//_su3_times_su3_acc((*v), tmp2, tmp1);
90+
_real_times_su3_times_su3_acc((*v), tmp2, tmp1, ptbc_fac2 * ptbc_fac3);
8491

8592
/* second contr. starting from x
8693
* a b c e^+ d^+
@@ -96,14 +103,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
96103
_su3_times_su3(tmp1, *a, *b);
97104
z = g_iup[y][mu];
98105
c = &gf[z][mu];
106+
double const ptbc_fac4 = get_ptbc_coeff(x, nu) * get_ptbc_coeff(y, mu) * get_ptbc_coeff(z, mu);
99107
_su3_times_su3(tmp2, tmp1, *c);
100108

101109
y = g_iup[x][mu];
102110
d = &gf[y][mu];
103111
z = g_iup[y][mu];
104112
e = &gf[z][nu];
113+
double const ptbc_fac5 = get_ptbc_coeff(y, mu) * get_ptbc_coeff(z, nu);
105114
_su3_times_su3(tmp1, *d, *e);
106-
_su3_times_su3d_acc((*v), tmp2, tmp1);
115+
//_su3_times_su3d_acc((*v), tmp2, tmp1);
116+
_real_times_su3_times_su3d_acc((*v), tmp2, tmp1, ptbc_fac4 * ptbc_fac5);
107117

108118
/* 1 contr. starting idn[x][nu]
109119
* d^+ a b c e^+
@@ -119,14 +129,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
119129
_su3d_times_su3(tmp1, *d, *a);
120130
z = g_iup[y][mu];
121131
b = &gf[z][mu];
132+
double const ptbc_fac6 = get_ptbc_coeff(y, nu) * get_ptbc_coeff(y, mu) * get_ptbc_coeff(z, mu);
122133
_su3_times_su3(tmp2, tmp1, *b);
123134

124135
y = g_iup[z][mu];
125136
c = &gf[y][nu];
126137
z = g_iup[x][mu];
127138
e = &gf[z][mu];
139+
double const ptbc_fac7 = get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, mu);
128140
_su3_times_su3d(tmp1, *c, *e);
129-
_su3_times_su3_acc((*v), tmp2, tmp1);
141+
//_su3_times_su3_acc((*v), tmp2, tmp1);
142+
_real_times_su3_times_su3_acc((*v), tmp2, tmp1, ptbc_fac6 * ptbc_fac7);
130143

131144
/* 1 contr. starting idn[idn[x][mu]][nu]
132145
* e^+ d^+ a b c
@@ -142,14 +155,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
142155
a = &gf[z][mu];
143156
_su3d_times_su3(tmp1, *d, *a);
144157
e = &gf[y][mu];
158+
double const ptbc_fac8 = get_ptbc_coeff(z, nu) * get_ptbc_coeff(z, mu) * get_ptbc_coeff(y, mu);
145159
_su3d_times_su3(tmp2, *e, tmp1);
146160

147161
y = g_idn[x][nu];
148162
b = &gf[y][mu];
149163
z = g_iup[y][mu];
150164
c = &gf[z][nu];
165+
double const ptbc_fac9 = get_ptbc_coeff(y, mu) * get_ptbc_coeff(z, nu);
151166
_su3_times_su3(tmp1, *b, *c);
152-
_su3_times_su3_acc((*v), tmp2, tmp1);
167+
//_su3_times_su3_acc((*v), tmp2, tmp1);
168+
_real_times_su3_times_su3_acc((*v), tmp2, tmp1, ptbc_fac8 * ptbc_fac9);
153169

154170
/* 1 contr. starting idn[x][mu]
155171
* d^+ a b c e^+
@@ -165,14 +181,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
165181
a = &gf[y][nu];
166182
_su3d_times_su3(tmp1, *d, *a);
167183
b = &gf[z][mu];
184+
double const ptbc_fac10 = get_ptbc_coeff(y, mu) * get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, mu);
168185
_su3_times_su3(tmp2, tmp1, *b);
169186

170187
y = g_iup[x][mu];
171188
e = &gf[y][nu];
172189
z = g_iup[x][nu];
173190
c = &gf[z][mu];
191+
double const ptbc_fac11 = get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, mu);
174192
_su3_times_su3d(tmp1, *c, *e);
175-
_su3_times_su3_acc((*v), tmp2, tmp1);
193+
//_su3_times_su3_acc((*v), tmp2, tmp1);
194+
_real_times_su3_times_su3_acc((*v), tmp2, tmp1, ptbc_fac10 * ptbc_fac11);
176195
}
177196
}
178197
}

‎src/lib/get_staples.c‎

Lines changed: 19 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@
2929
#include "start.h"
3030
#include "su3.h"
3131
#include "su3adj.h"
32+
#include "ptbc.h"
3233

3334
void get_staples(su3* const staple, const int x, const int mu, const su3** in_gauge_field) {
3435
int iy;
@@ -41,20 +42,24 @@ void get_staples(su3* const staple, const int x, const int mu, const su3** in_ga
4142
w1 = &in_gauge_field[x][k];
4243
w2 = &in_gauge_field[g_iup[x][k]][mu];
4344
w3 = &in_gauge_field[g_iup[x][mu]][k];
45+
double const ptbc_fac0 = get_ptbc_coeff(x, k) * get_ptbc_coeff(g_iup[x][k], mu) * get_ptbc_coeff(g_iup[x][mu], k);
4446

4547
/* st = w2 * w3^d */
4648
_su3_times_su3d(st, *w2, *w3);
4749
/* v = v + w1 * st */
48-
_su3_times_su3_acc(*staple, *w1, st);
50+
//_su3_times_su3_acc(*staple, *w1, st);
51+
_real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac0);
4952

5053
iy = g_idn[x][k];
5154
w1 = &in_gauge_field[iy][k];
5255
w2 = &in_gauge_field[iy][mu];
5356
w3 = &in_gauge_field[g_iup[iy][mu]][k];
57+
double const ptbc_fac1 = get_ptbc_coeff(iy, k) * get_ptbc_coeff(iy, mu) * get_ptbc_coeff(g_iup[iy][mu], k);
5458
/* st = w2 * w3 */
5559
_su3_times_su3(st, *w2, *w3);
5660
/* v = v + w1^d * st */
57-
_su3d_times_su3_acc(*staple, *w1, st);
61+
//_su3d_times_su3_acc(*staple, *w1, st);
62+
_real_times_su3d_times_su3_acc(*staple, *w1, st, ptbc_fac1);
5863
}
5964
}
6065
}
@@ -71,20 +76,24 @@ void get_spacelike_staples(su3* const staple, const int x, const int mu,
7176
w1 = &in_gauge_field[x][k];
7277
w2 = &in_gauge_field[g_iup[x][k]][mu];
7378
w3 = &in_gauge_field[g_iup[x][mu]][k];
79+
double const ptbc_fac0 = get_ptbc_coeff(x, k) * get_ptbc_coeff(g_iup[x][k], mu) * get_ptbc_coeff(g_iup[x][mu], k);
7480

7581
/* st = w2 * w3^d */
7682
_su3_times_su3d(st, *w2, *w3);
7783
/* v = v + w1 * st */
78-
_su3_times_su3_acc(*staple, *w1, st);
84+
//_su3_times_su3_acc(*staple, *w1, st);
85+
_real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac0);
7986

8087
iy = g_idn[x][k];
8188
w1 = &in_gauge_field[iy][k];
8289
w2 = &in_gauge_field[iy][mu];
8390
w3 = &in_gauge_field[g_iup[iy][mu]][k];
91+
double const ptbc_fac1 = get_ptbc_coeff(iy, k) * get_ptbc_coeff(iy, mu) * get_ptbc_coeff(g_iup[iy][mu], k);
8492
/* st = w2 * w3 */
8593
_su3_times_su3(st, *w2, *w3);
8694
/* v = v + w1^d * st */
87-
_su3d_times_su3_acc(*staple, *w1, st);
95+
//_su3d_times_su3_acc(*staple, *w1, st);
96+
_real_times_su3d_times_su3_acc(*staple, *w1, st, ptbc_fac1);
8897
}
8998
}
9099
}
@@ -101,19 +110,23 @@ void get_timelike_staples(su3* const staple, const int x, const int mu,
101110
w1 = &in_gauge_field[x][k];
102111
w2 = &in_gauge_field[g_iup[x][k]][mu];
103112
w3 = &in_gauge_field[g_iup[x][mu]][k];
113+
double const ptbc_fac0 = get_ptbc_coeff(x, k) * get_ptbc_coeff(g_iup[x][k], mu) * get_ptbc_coeff(g_iup[x][mu], k);
104114

105115
/* st = w2 * w3^d */
106116
_su3_times_su3d(st, *w2, *w3);
107117
/* v = v + w1 * st */
108-
_su3_times_su3_acc(*staple, *w1, st);
118+
//_su3_times_su3_acc(*staple, *w1, st);
119+
_real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac0);
109120

110121
iy = g_idn[x][k];
111122
w1 = &in_gauge_field[iy][k];
112123
w2 = &in_gauge_field[iy][mu];
113124
w3 = &in_gauge_field[g_iup[iy][mu]][k];
125+
double const ptbc_fac1 = get_ptbc_coeff(iy, k) * get_ptbc_coeff(iy, mu) * get_ptbc_coeff(g_iup[iy][mu], k);
114126
/* st = w2 * w3 */
115127
_su3_times_su3(st, *w2, *w3);
116128
/* v = v + w1^d * st */
117-
_su3d_times_su3_acc(*staple, *w1, st);
129+
//_su3d_times_su3_acc(*staple, *w1, st);
130+
_real_times_su3d_times_su3_acc(*staple, *w1, st, ptbc_fac1);
118131
}
119132
}

‎src/lib/measure_gauge_action.c‎

Lines changed: 14 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,7 @@
3939
#include "measure_gauge_action.h"
4040
#include "su3.h"
4141
#include "su3adj.h"
42+
#include "ptbc.h"
4243

4344
double measure_plaquette(const su3 *const *const gf) {
4445
static double res;
@@ -73,8 +74,11 @@ double measure_plaquette(const su3 *const *const gf) {
7374
v = &gf[ix][mu2];
7475
w = &gf[ix2][mu1];
7576
_su3_times_su3(pr2, *v, *w);
77+
// compute local parallel tempering factor: fac_1 * fac_2 * fac_3 * fac_4
78+
double const ptbc_fac = get_ptbc_coeff(ix, mu1) * get_ptbc_coeff(ix1, mu2) * get_ptbc_coeff(ix, mu2) * get_ptbc_coeff(ix2, mu1);
79+
//double const ptbc_fac = 1.;
7680
_trace_su3_times_su3d(ac, pr1, pr2);
77-
tr = ac + kc;
81+
tr = ac*ptbc_fac + kc;
7882
ts = tr + ks;
7983
tt = ts - ks;
8084
ks = ts;
@@ -135,9 +139,12 @@ double measure_gauge_action(const su3 *const *const gf, const double lambda) {
135139
v = &gf[ix][mu2];
136140
w = &gf[ix2][0];
137141
_su3_times_su3(pr2, *v, *w);
142+
// parallel tempering factor = fac_1 * fac_2 * fac_3 * fac_4 computed in 2 parts
143+
double const ptbc_fac = get_ptbc_coeff(ix, 0) * get_ptbc_coeff(ix1, mu2) * get_ptbc_coeff(ix, mu2) * get_ptbc_coeff(ix2, 0);
138144
_trace_su3_times_su3d(ac, pr1, pr2);
139145
ac *= (1 + lambda);
140-
tr = ac + kc;
146+
tr = ac*ptbc_fac + kc;
147+
//tr = ac + kc;
141148
ts = tr + ks;
142149
tt = ts - ks;
143150
ks = ts;
@@ -146,6 +153,7 @@ double measure_gauge_action(const su3 *const *const gf, const double lambda) {
146153
// magnetic part
147154
for (int mu1 = 1; mu1 < 3; mu1++) {
148155
ix1 = g_iup[ix][mu1];
156+
double const ptbc_fac1 = get_ptbc_coeff(ix, mu1);
149157
for (int mu2 = mu1 + 1; mu2 < 4; mu2++) {
150158
ix2 = g_iup[ix][mu2];
151159
v = &gf[ix][mu1];
@@ -154,9 +162,11 @@ double measure_gauge_action(const su3 *const *const gf, const double lambda) {
154162
v = &gf[ix][mu2];
155163
w = &gf[ix2][mu1];
156164
_su3_times_su3(pr2, *v, *w);
165+
double const ptbc_fac2 = get_ptbc_coeff(ix1, mu2) * get_ptbc_coeff(ix, mu2) * get_ptbc_coeff(ix2, mu1);
157166
_trace_su3_times_su3d(ac, pr1, pr2);
158-
ac *= (1 - lambda);
159-
tr = ac + kc;
167+
//ac *= (1 - lambda);
168+
tr = ac*(1 - lambda)*(ptbc_fac2*ptbc_fac1) + kc;
169+
//tr = ac + kc;
160170
ts = tr + ks;
161171
tt = ts - ks;
162172
ks = ts;

‎src/lib/measure_rectangles.c‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,7 @@
4545
#include "measure_rectangles.h"
4646
#include "su3.h"
4747
#include "su3adj.h"
48+
#include "ptbc.h"
4849

4950
double measure_rectangles(const su3 **const gf) {
5051
static double res;
@@ -86,6 +87,7 @@ double measure_rectangles(const su3 **const gf) {
8687
_su3_times_su3(tmp, *v, *w);
8788
v = &gf[k][nu];
8889
_su3_times_su3(pr1, tmp, *v);
90+
double const ptbc_fac0 = get_ptbc_coeff(i, mu) * get_ptbc_coeff(j, nu) * get_ptbc_coeff(k, nu);
8991
/*
9092
->
9193
^
@@ -100,12 +102,14 @@ double measure_rectangles(const su3 **const gf) {
100102
_su3_times_su3(tmp, *v, *w);
101103
v = &gf[k][mu];
102104
_su3_times_su3(pr2, tmp, *v);
105+
double const ptbc_fac1 = get_ptbc_coeff(i, nu) * get_ptbc_coeff(j, nu) * get_ptbc_coeff(k, mu);
103106

104107
/* Trace it */
105108
_trace_su3_times_su3d(ac, pr1, pr2);
106109
/* printf("i mu nu: %d %d %d, ac = %e\n", i, mu, nu, ac); */
107110
/* Kahan summation */
108-
tr = ac + kc;
111+
tr = ac*(ptbc_fac0 * ptbc_fac1) + kc;
112+
//tr = ac + kc;
109113
ts = tr + ks;
110114
tt = ts - ks;
111115
ks = ts;

‎src/lib/su3.h‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -490,6 +490,17 @@ typedef struct {
490490
(u).c21 += (v).c20 * (w).c01 + (v).c21 * (w).c11 + (v).c22 * (w).c21; \
491491
(u).c22 += (v).c20 * (w).c02 + (v).c21 * (w).c12 + (v).c22 * (w).c22;
492492

493+
#define _real_times_su3_times_su3_acc(u, v, w, r) \
494+
(u).c00 += (r) * ((v).c00 * (w).c00 + (v).c01 * (w).c10 + (v).c02 * (w).c20); \
495+
(u).c01 += (r) * ((v).c00 * (w).c01 + (v).c01 * (w).c11 + (v).c02 * (w).c21); \
496+
(u).c02 += (r) * ((v).c00 * (w).c02 + (v).c01 * (w).c12 + (v).c02 * (w).c22); \
497+
(u).c10 += (r) * ((v).c10 * (w).c00 + (v).c11 * (w).c10 + (v).c12 * (w).c20); \
498+
(u).c11 += (r) * ((v).c10 * (w).c01 + (v).c11 * (w).c11 + (v).c12 * (w).c21); \
499+
(u).c12 += (r) * ((v).c10 * (w).c02 + (v).c11 * (w).c12 + (v).c12 * (w).c22); \
500+
(u).c20 += (r) * ((v).c20 * (w).c00 + (v).c21 * (w).c10 + (v).c22 * (w).c20); \
501+
(u).c21 += (r) * ((v).c20 * (w).c01 + (v).c21 * (w).c11 + (v).c22 * (w).c21); \
502+
(u).c22 += (r) * ((v).c20 * (w).c02 + (v).c21 * (w).c12 + (v).c22 * (w).c22);
503+
493504
#define _su3_times_su3d(u, v, w) \
494505
(u).c00 = (v).c00 * conj((w).c00) + (v).c01 * conj((w).c01) + (v).c02 * conj((w).c02); \
495506
(u).c01 = (v).c00 * conj((w).c10) + (v).c01 * conj((w).c11) + (v).c02 * conj((w).c12); \
@@ -512,6 +523,17 @@ typedef struct {
512523
(u).c21 += (v).c20 * conj((w).c10) + (v).c21 * conj((w).c11) + (v).c22 * conj((w).c12); \
513524
(u).c22 += (v).c20 * conj((w).c20) + (v).c21 * conj((w).c21) + (v).c22 * conj((w).c22);
514525

526+
#define _real_times_su3_times_su3d_acc(u, v, w, r) \
527+
(u).c00 += (r) * ((v).c00 * conj((w).c00) + (v).c01 * conj((w).c01) + (v).c02 * conj((w).c02)); \
528+
(u).c01 += (r) * ((v).c00 * conj((w).c10) + (v).c01 * conj((w).c11) + (v).c02 * conj((w).c12)); \
529+
(u).c02 += (r) * ((v).c00 * conj((w).c20) + (v).c01 * conj((w).c21) + (v).c02 * conj((w).c22)); \
530+
(u).c10 += (r) * ((v).c10 * conj((w).c00) + (v).c11 * conj((w).c01) + (v).c12 * conj((w).c02)); \
531+
(u).c11 += (r) * ((v).c10 * conj((w).c10) + (v).c11 * conj((w).c11) + (v).c12 * conj((w).c12)); \
532+
(u).c12 += (r) * ((v).c10 * conj((w).c20) + (v).c11 * conj((w).c21) + (v).c12 * conj((w).c22)); \
533+
(u).c20 += (r) * ((v).c20 * conj((w).c00) + (v).c21 * conj((w).c01) + (v).c22 * conj((w).c02)); \
534+
(u).c21 += (r) * ((v).c20 * conj((w).c10) + (v).c21 * conj((w).c11) + (v).c22 * conj((w).c12)); \
535+
(u).c22 += (r) * ((v).c20 * conj((w).c20) + (v).c21 * conj((w).c21) + (v).c22 * conj((w).c22));
536+
515537
#define _su3d_times_su3(u, v, w) \
516538
(u).c00 = conj((v).c00) * (w).c00 + conj((v).c10) * (w).c10 + conj((v).c20) * (w).c20; \
517539
(u).c01 = conj((v).c00) * (w).c01 + conj((v).c10) * (w).c11 + conj((v).c20) * (w).c21; \
@@ -534,6 +556,17 @@ typedef struct {
534556
(u).c21 += conj((v).c02) * (w).c01 + conj((v).c12) * (w).c11 + conj((v).c22) * (w).c21; \
535557
(u).c22 += conj((v).c02) * (w).c02 + conj((v).c12) * (w).c12 + conj((v).c22) * (w).c22;
536558

559+
#define _real_times_su3d_times_su3_acc(u, v, w, r) \
560+
(u).c00 += (r) * (conj((v).c00) * (w).c00 + conj((v).c10) * (w).c10 + conj((v).c20) * (w).c20); \
561+
(u).c01 += (r) * (conj((v).c00) * (w).c01 + conj((v).c10) * (w).c11 + conj((v).c20) * (w).c21); \
562+
(u).c02 += (r) * (conj((v).c00) * (w).c02 + conj((v).c10) * (w).c12 + conj((v).c20) * (w).c22); \
563+
(u).c10 += (r) * (conj((v).c01) * (w).c00 + conj((v).c11) * (w).c10 + conj((v).c21) * (w).c20); \
564+
(u).c11 += (r) * (conj((v).c01) * (w).c01 + conj((v).c11) * (w).c11 + conj((v).c21) * (w).c21); \
565+
(u).c12 += (r) * (conj((v).c01) * (w).c02 + conj((v).c11) * (w).c12 + conj((v).c21) * (w).c22); \
566+
(u).c20 += (r) * (conj((v).c02) * (w).c00 + conj((v).c12) * (w).c10 + conj((v).c22) * (w).c20); \
567+
(u).c21 += (r) * (conj((v).c02) * (w).c01 + conj((v).c12) * (w).c11 + conj((v).c22) * (w).c21); \
568+
(u).c22 += (r) * (conj((v).c02) * (w).c02 + conj((v).c12) * (w).c12 + conj((v).c22) * (w).c22);
569+
537570
#define _su3_minus_const_times_im_trace_su3(w, c, v) \
538571
(w).c00 -= I * c * (cimag((v).c00) + cimag((v).c11) + cimag((v).c22)); \
539572
(w).c11 -= I * c * (cimag((v).c00) + cimag((v).c11) + cimag((v).c22)); \

0 commit comments

Comments
 (0)