HILA
Loading...
Searching...
No Matches
stout_smear_nch.h
1#ifndef STOUT_SMEAR_NCH_H_
2#define STOUT_SMEAR_NCH_H_
3
4#include "hila.h"
5
7
8/////////////////////////////////////////////////////////////////////////////
9/// Do stout smearing for the gauge fields using non-Cayley-Hamilton exponentiation
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 nch_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 nch_stout_smear1(const GaugeField<T> &U, out_only GaugeField<T> &stout, atype coeff) {
55 nch_staplesums(U, stout);
56 foralldir(d) {
57 onsites(ALL) {
58 stout[d][X] =
59 exp((U[d][X] * stout[d][X]).project_to_algebra_scaled(-coeff).expand()) * U[d][X];
60 stout[d][X].reunitarize();
61 }
62 }
63}
64
65template <typename T, typename atype = hila::arithmetic_type<T>>
66void nch_stout_smear1(const GaugeField<T> &U, out_only GaugeField<T> &stout,
67 out_only VectorField<T> &stap, atype coeff) {
68 nch_staplesums(U, stap);
69 foralldir(d) {
70 onsites(ALL) {
71 stout[d][X] =
72 exp((U[d][X] * stap[d][X]).project_to_algebra_scaled(-coeff).expand()) * U[d][X];
73 stout[d][X].reunitarize();
74 }
75 }
76}
77
78
79template <typename T, typename atype = hila::arithmetic_type<T>>
80void nch_stout_smear(const GaugeField<T> &U, GaugeField<T> &stout, atype coeff, int iter) {
81 GaugeField<T> tmp;
82 if (iter % 2 == 0) {
83 stout = U;
84 for (int i = 0; i < iter / 2; ++i) {
85 nch_stout_smear1(stout, tmp, coeff);
86 nch_stout_smear1(tmp, stout, coeff);
87 }
88 } else {
89 tmp = U;
90 for (int i = 0; i < iter / 2; ++i) {
91 nch_stout_smear1(tmp, stout, coeff);
92 nch_stout_smear1(stout, tmp, coeff);
93 }
94 nch_stout_smear1(tmp, stout, coeff);
95 }
96}
97
98
99template <typename T, typename atype = hila::arithmetic_type<T>>
100void nch_stout_smear(const GaugeField<T> &U, out_only std::vector<GaugeField<T>> &stoutlist,
101 out_only std::vector<VectorField<T>> &staplist, atype coeff) {
102 // performs nst stout smearing steps on the gauge field U and stores the gague field obtained
103 // after the i-th smearing step in stoutlist[i]
104 stoutlist[0] = U; // original field
105 for (int i = 1; i < stoutlist.size(); ++i) {
106 nch_stout_smear1(stoutlist[i - 1], stoutlist[i], staplist[i - 1], coeff);
107 }
108}
109
110template <typename T, typename atype = hila::arithmetic_type<T>>
111void nch_stout_smear_force(const std::vector<GaugeField<T>> &stoutlist,
112 const std::vector<VectorField<T>> &staplist,
113 const VectorField<Algebra<T>> &K, out_only VectorField<Algebra<T>> &KS,
114 atype coeff) {
115 // uses the list stoutlist[] of smeared gauge fields to compute the pullback of the
116 // algebra-valued force field K under the smearing and returns the force acting on the unsmeared
117 // link variables as algebra-valued field KS
118 // Note: our definition of the force field is different from the one used in
119 // [arXiv:hep-lat/0311018v1], in order to match the force field representation used by our HMC
120 // implementation.
122 Field<T> K21, K22, K23, K24;
123
124 foralldir(d1) onsites(ALL) KS[d1][X] = K[d1][X];
125
126 for (int i = stoutlist.size() - 2; i >= 0; --i) {
127 const GaugeField<T> &U0 = stoutlist[i + 1];
128 const GaugeField<T> &U = stoutlist[i];
129 const VectorField<T> &staps = staplist[i];
130
131 foralldir(d1) {
132 //get_stout_staples(U, d1, stapl, staps);
133 onsites(ALL) {
134 // compute stout smearing operator and its derivatives:
135
136 // temp. variables:
137 T mtexp;
138 T mdtexp;
139 // turn staple sum into plaquette sum by multiplying with link variable:
140 //T tplaqs = U[d1][X] * staps[X];
141 T tplaqs = U[d1][X] * staps[d1][X];
142 // the following function computes first for X = -coeff * tplaqs the smearing
143 // operator Q = exp(X) and its derivatives dQ/dX[][], and uses these to
144 // compute the two matrices:
145 // mtexp = Q.dagger() * KS[d1][X].expand() * Q
146 // and
147 // mdtexp[i][j] = trace(Q.dagger * KS[d1][X].expand() * dQ/dX[j][i]) :
148 mult_exp(tplaqs.project_to_algebra_scaled(-coeff).expand(),
149 U[d1][X] * U0[d1][X].dagger() * KS[d1][X].expand(), mtexp, mdtexp);
150
151 // set K1[d1][X] to be the equivalent of the \Lambda matrix from eq.(73) in
152 // [arXiv:hep-lat/0311018v1]:
153 K1[d1][X] = mdtexp.project_to_algebra_scaled(coeff).expand();
154
155 // equivalent of first line and first term on second line of eq.(75) in
156 // [arXiv:hep-lat/0311018v1]:
157 KS[d1][X] = (mtexp - tplaqs * K1[d1][X]).project_to_algebra();
158
159 // multiply K1[d1] by U[d1]:
160 K1[d1][X] *= U[d1][X];
161 }
162 }
163
164 // equivalent of remaining terms of eq.(75) in [arXiv:hep-lat/0311018v1]:
165 foralldir(d1) foralldir(d2) if (d1 != d2) {
166 onsites(ALL) {
167 T U2, U4, tM1;
168
169 U2 = U[d2][X + d1];
170 U4 = U[d2][X].dagger();
171
172 tM1 = U2 * U[d1][X + d2].dagger();
173 K21[X] = U4 * K1[d1][X] * tM1;
174 tM1 *= U4;
175 K22[X] = tM1 * K1[d1][X];
176 K24[X] = K1[d1][X] * tM1;
177
178 tM1 = U2 * K1[d1][X + d2].dagger() * U4;
179 K23[X] = U[d1][X] * tM1;
180 K22[X] += tM1 * U[d1][X];
181
182 K24[X] += K23[X];
183 }
184
185 onsites(ALL) {
186 KS[d2][X] -= (K22[X - d1] - K24[X]).project_to_algebra();
187 KS[d1][X] -= (K23[X] - K21[X - d2]).project_to_algebra();
188 }
189 }
190 }
191}
192
193#endif
Array< n, m, T > exp(Array< n, m, T > a)
Exponential.
Definition array.h:1059
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