HILA
Loading...
Searching...
No Matches
stout_smear_nchm.h
1#ifndef STOUT_SMEAR_NCHM_H_
2#define STOUT_SMEAR_NCHM_H_
3
4#include "hila.h"
5
7
8/////////////////////////////////////////////////////////////////////////////
9/// Do stout smearing for the gauge fields
10/// U' = exp( -coeff * Proj_to_Alg(U * Staples) ) * U
11
12/**
13 * @brief Compute sum of staples around all links
14 * @details Compute staple sums around all positively oriented (i.e. corresponding
15 * plaquette sums are obtained by multiplying the staple sums by the corresponding
16 * positively oriented link variables)
17
18 * @tparam T Matrix element type
19 * @param U input gauge field of type T
20 * @param staples output vector/gauge field of type T
21 * @return void
22 */
23template <typename T>
24void nchm_staplesums(const GaugeField<T> &U, out_only GaugeField<T> &staples) {
25 Field<T> tF1, lw1;
26 foralldir(d1) {
27 foralldir(d2) if (d2 != d1) {
28 U[d2].start_gather(d1, ALL);
29 }
30 onsites(ALL) staples[d1][X] = 0;
31 }
32
33 foralldir(d1) foralldir(d2) if (d2 != d1) {
34 onsites(ALL) {
35 // tF1[X] = (U[d1][X] * U[d2][X + d1]).dagger();
36 mult_aa(U[d1][X], U[d2][X + d1], tF1[X]);
37 // lw1[X] = tF1[X] * U[d2][X];
38 mult(tF1[X], U[d2][X], lw1[X]);
39 }
40
41 lw1.start_gather(-d2, ALL);
42
43 onsites(ALL) {
44 // staple[d2][X] += U[d1][X + d2] * tF1[X];
45 mult_add(U[d1][X + d2], tF1[X], staples[d2][X]);
46 }
47 onsites(ALL) {
48 staples[d1][X] += lw1[X - d2];
49 }
50 }
51}
52
53template <typename T, typename atype = hila::arithmetic_type<T>>
54void nchm_stout_smear1(const GaugeField<T> &U, out_only GaugeField<T> &stout,
55 out_only std::array<Field<int>, NDIM> &niter, atype coeff) {
56 nchm_staplesums(U, stout);
57 foralldir(d) {
58 onsites(ALL) {
59 stout[d][X] = altexp((U[d][X] * stout[d][X]).project_to_algebra_scaled(-coeff).expand(),
60 niter[d][X]) * U[d][X];
61 }
62 }
63}
64
65template <typename T, typename atype = hila::arithmetic_type<T>>
66void nchm_stout_smear1(const GaugeField<T> &U, out_only GaugeField<T> &stout,
67 out_only VectorField<T> &stap,
68 out_only std::array<Field<int>, NDIM> &niter, atype coeff) {
69 nchm_staplesums(U, stap);
70 foralldir(d) {
71 onsites(ALL) {
72 stout[d][X] = altexp((U[d][X] * stap[d][X]).project_to_algebra_scaled(-coeff).expand(),
73 niter[d][X]) *
74 U[d][X];
75 }
76 }
77}
78
79template <typename T, typename atype = hila::arithmetic_type<T>>
80void nchm_stout_smear(const GaugeField<T> &U, out_only std::vector<GaugeField<T>> &stoutlist,
81 out_only std::vector<GaugeField<T>> &staplist,
82 out_only std::vector<std::array<Field<int>, NDIM>> &niterlist, atype coeff) {
83 // performs nst stout smearing steps on the gauge field U and stores the gague field obtained
84 // after the i-th smearing step in stoutlist[i]
85 stoutlist[0] = U; // original field
86 for (int i = 1; i < stoutlist.size(); ++i) {
87 nchm_stout_smear1(stoutlist[i - 1], stoutlist[i], staplist[i - 1], niterlist[i - 1], coeff);
88 }
89}
90
91template <typename T, typename atype = hila::arithmetic_type<T>>
92void nchm_stout_smear1(const GaugeField<T> &U, out_only GaugeField<T> &stout, atype coeff) {
93 nchm_staplesums(U, stout);
94 foralldir(d) {
95 onsites(ALL) {
96 stout[d][X] =
97 altexp((U[d][X] * stout[d][X]).project_to_algebra_scaled(-coeff).expand()) *
98 U[d][X];
99 }
100 }
101}
102
103template <typename T, typename atype = hila::arithmetic_type<T>>
104void nchm_stout_smear1(const GaugeField<T> &U, out_only GaugeField<T> &stout,
105 out_only GaugeField<T> &stap, atype coeff) {
106 nchm_staplesums(U, stap);
107 foralldir(d) {
108 onsites(ALL) {
109 stout[d][X] =
110 altexp((U[d][X] * stap[d][X]).project_to_algebra_scaled(-coeff).expand()) * U[d][X];
111 }
112 }
113}
114
115template <typename T, typename atype = hila::arithmetic_type<T>>
116void nchm_stout_smear(const GaugeField<T> &U, GaugeField<T> &stout, atype coeff, int iter) {
117 GaugeField<T> tmp;
118 if (iter % 2 == 0) {
119 stout = U;
120 for (int i = 0; i < iter / 2; ++i) {
121 nchm_stout_smear1(stout, tmp, coeff);
122 nchm_stout_smear1(tmp, stout, coeff);
123 }
124 } else {
125 tmp = U;
126 for (int i = 0; i < iter / 2; ++i) {
127 nchm_stout_smear1(tmp, stout, coeff);
128 nchm_stout_smear1(stout, tmp, coeff);
129 }
130 nchm_stout_smear1(tmp, stout, coeff);
131 }
132}
133
134template <typename T, typename atype = hila::arithmetic_type<T>>
135void nchm_stout_smear(const GaugeField<T> &U, GaugeField<T> &stout,
136 std::vector<std::array<Field<int>, NDIM>> &niterlist, atype coeff) {
137 GaugeField<T> tmp;
138 int iter = niterlist.size();
139 if (iter % 2 == 0) {
140 stout = U;
141 for (int i = 0; i < iter / 2; ++i) {
142 nchm_stout_smear1(stout, tmp, niterlist[2 * i], coeff);
143 nchm_stout_smear1(tmp, stout, niterlist[2 * i + 1], coeff);
144 }
145 } else {
146 tmp = U;
147 for (int i = 0; i < iter / 2; ++i) {
148 nchm_stout_smear1(tmp, stout, niterlist[2 * i], coeff);
149 nchm_stout_smear1(stout, tmp, niterlist[2 * i + 1], coeff);
150 }
151 nchm_stout_smear1(tmp, stout, niterlist[iter - 1], coeff);
152 }
153}
154
155
156template <typename T, int N = T::rows(), int NSTP = 2 * (NDIM - 1),
157 typename atype = hila::arithmetic_type<T>>
158void nchm_stout_smear_force(const std::vector<GaugeField<T>> &stoutlist,
159 const std::vector<VectorField<T>> &staplist,
160 const VectorField<Algebra<T>> &K, out_only VectorField<Algebra<T>> &KS,
161 std::vector<std::array<Field<int>, NDIM>> &niterlist, atype coeff) {
162 // uses the list stoutlist[] of smeared gauge fields to compute the pullback of the
163 // algebra-valued force field K under the smearing and returns the force acting on the unsmeared
164 // link variables as algebra-valued field KS
165 // Note: our definition of the force field is different from the one used in
166 // [arXiv:hep-lat/0311018v1], in order to match the force field representation used by our HMC
167 // implementation.
169 Field<T> K21, K22, K23, K24;
170
171 foralldir(d1) onsites(ALL) KS[d1][X] = K[d1][X];
172
173 for (int i = stoutlist.size() - 2; i >= 0; --i) {
174 const GaugeField<T> &U0 = stoutlist[i + 1];
175 const GaugeField<T> &U = stoutlist[i];
176 const VectorField<T> &staps = staplist[i];
177 const std::array<Field<int>, NDIM> &niter = niterlist[i];
178
179 foralldir(d1) {
180 //get_stout_staples(U, d1, stapl, staps);
181 onsites(ALL) {
182 // compute stout smearing operator and its derivatives:
183
184 // temp. variables:
185 T mtexp;
186 T mdtexp;
187 // turn staple sum into plaquette sum by multiplying with link variable:
188 //T tplaqs = U[d1][X] * staps[X];
189 T tplaqs = U[d1][X] * staps[d1][X];
190 // the following function computes first for X = -coeff * tplaqs the smearing
191 // operator Q = exp(X) and its derivatives dQ/dX[][], and uses these to
192 // compute the two matrices:
193 // mtexp = Q.dagger() * KS[d1][X].expand() * Q
194 // and
195 // mdtexp[i][j] = trace(Q.dagger * KS[d1][X].expand() * dQ/dX[j][i]) :
196 mult_exp(tplaqs.project_to_algebra_scaled(-coeff).expand(),
197 U[d1][X] * U0[d1][X].dagger() * KS[d1][X].expand(), mtexp, mdtexp, niter[d1][X]);
198
199 // set K1[d1][X] to be the equivalent of the \Lambda matrix from eq.(73) in
200 // [arXiv:hep-lat/0311018v1]:
201 K1[d1][X] = mdtexp.project_to_algebra_scaled(coeff).expand();
202
203 // equivalent of first line and first term on second line of eq.(75) in
204 // [arXiv:hep-lat/0311018v1]:
205 KS[d1][X] = (mtexp - tplaqs * K1[d1][X]).project_to_algebra();
206
207 // multiply K1[d1] by U[d1]:
208 K1[d1][X] *= U[d1][X];
209 }
210 }
211
212 // equivalent of remaining terms of eq.(75) in [arXiv:hep-lat/0311018v1]:
213 foralldir(d1) foralldir(d2) if (d1 != d2) {
214 onsites(ALL) {
215 T U2, U4, tM1;
216
217 U2 = U[d2][X + d1];
218 U4 = U[d2][X].dagger();
219
220 tM1 = U2 * U[d1][X + d2].dagger();
221 K21[X] = U4 * K1[d1][X] * tM1;
222 tM1 *= U4;
223 K22[X] = tM1 * K1[d1][X];
224 K24[X] = K1[d1][X] * tM1;
225
226 tM1 = U2 * K1[d1][X + d2].dagger() * U4;
227 K23[X] = U[d1][X] * tM1;
228 K22[X] += tM1 * U[d1][X];
229
230 K24[X] += K23[X];
231 }
232
233 onsites(ALL) {
234 KS[d2][X] -= (K22[X - d1] - K24[X]).project_to_algebra();
235 KS[d1][X] -= (K23[X] - K21[X - d2]).project_to_algebra();
236 }
237 }
238 }
239}
240
241#endif
Definition su2.h:7
The field class implements the standard methods for accessing Fields. Hilapp replaces the parity acce...
Definition field.h:61
Field< A > dagger() const
Returns dagger or Hermitian conjugate of Field depending on how it is defined for Field type T.
Definition field.h:1123
dir_mask_t start_gather(Direction d, Parity p=ALL) const
Definition field_comm.h:262
Gauge field class.
Definition gaugefield.h:22
#define foralldir(d)
Macro to loop over (all) Direction(s)
Definition coordinates.h:78
constexpr Parity ALL
bit pattern: 011
void mult_exp(const Matrix_t< n, m, T, MT > &mat, const Matrix_t< n, m, T, MT > &mmat, Matrix_t< n, m, T, MT > &r, Matrix_t< n, m, T, MT > &dr, const int order=20)
Calculate mmat*exp(mat) and trace(mmat*dexp(mat))
Definition matrix.h:2888
void mult_add(const Mt1 &a, const Mt2 &b, Mt3 &res)
compute product of two matrices and add result to existing matrix
Definition matrix.h:2437
void mult_aa(const Mt1 &a, const Mt2 &b, Mt3 &res)
compute hermitian conjugate of product of two matrices and write result to existing matrix
Definition matrix.h:2356
void mult(const Mt1 &a, const Mt2 &b, Mt3 &res)
compute product of two matrices and write result to existing matrix
Definition matrix.h:2327
Matrix_t< n, m, T, MT > altexp(const Matrix_t< n, m, T, MT > &mat, int &niter)
Calculate exp of a square matrix.
Definition matrix.h:2820