HILA
Loading...
Searching...
No Matches
test_CG.cpp
1#include "test.h"
2
3#include "dirac/staggered.h"
4#include "dirac/wilson.h"
5#include "dirac/Hasenbusch.h"
6#include "dirac/conjugate_gradient.h"
7
8#define N 3
9
10void test_gamma_matrices() {
11 Wilson_vector<N, double> w1, w2, w3;
12 SU<N> U;
13 U.random();
14 w1.gaussian_random();
15
16#if NDIM == 4
17 w2 = w1 - gamma5 * (gamma5 * w1);
18 assert(w2.squarenorm() < 0.0001 && "g5*g5 = 1");
19#endif
20
21 foralldir(d) {
22 half_Wilson_vector<N, double> h1;
23 w2 = w1 - gamma_matrix[d] * (gamma_matrix[d] * w1);
24 assert(w2.squarenorm() < 0.0001 && "gamma_d*gamma_d = 1");
25
26 w2 = w1 + gamma_matrix[d] * w1;
27 h1 = half_Wilson_vector(w1, d, 1);
28 double diff = w2.squarenorm() - 2 * h1.squarenorm();
29 assert(diff * diff < 0.0001 && "half_Wilson_vector projection +1 norm");
30
31 w3 = (U * h1).expand(d, 1) - U * w2;
32 assert(w3.squarenorm() < 0.0001 && "half_wilson_vector expand");
33
34 w2 = w1 - gamma_matrix[d] * w1;
35 h1 = half_Wilson_vector(w1, d, -1);
36 diff = w2.squarenorm() - 2 * h1.squarenorm();
37 assert(diff * diff < 0.0001 && "half_Wilson_vector projection -1 norm");
38
39 w3 = (U * h1).expand(d, -1) - U * w2;
40 assert(w3.squarenorm() < 0.0001 && "half_wilson_vector expand");
41 }
42}
43
44int main(int argc, char **argv) {
45
46#if NDIM == 1
47 const CoordinateVector nd = {64};
48#elif NDIM == 2
49 const CoordinateVector nd = {32, 8};
50#elif NDIM == 3
51 const CoordinateVector nd = {16, 8, 8};
52#elif NDIM == 4
53 const CoordinateVector nd = {16, 8, 8, 8};
54#endif
55 hila::initialize(argc, argv);
56 lattice.setup(nd);
57
58 hila::seed_random(2);
59
60 test_gamma_matrices();
61
62 Field<SU<N>> U[NDIM];
63 foralldir(d){onsites(ALL){U[d][X] = 1;
64}
65}
66
67// Check conjugate of the staggered Dirac operator
68{
69 hila::out0 << "Checking with dirac_staggered\n";
70 using dirac = dirac_staggered<SU<N>>;
71 dirac D(0.1, U);
72 Field<SU_vector<N, double>> a, b, Db, Ddaggera, DdaggerDb;
73 onsites(ALL) {
74 a[X].gaussian_random();
75 b[X].gaussian_random();
76 }
77
78 double diffre = 0, diffim = 0;
79 D.apply(b, Db);
80 D.dagger(a, Ddaggera);
81 onsites(ALL) {
82 diffre += a[X].dot(Db[X]).re - Ddaggera[X].dot(b[X]).re;
83 diffim += a[X].dot(Db[X]).im - Ddaggera[X].dot(b[X]).im;
84 }
85
86 assert(diffre * diffre < 1e-16 && "test dirac_staggered");
87 assert(diffim * diffim < 1e-16 && "test dirac_staggered");
88
89 // Now run CG on Ddaggera. b=1/D a -> Db = a
90 CG<dirac> inverse(D);
91 onsites(ALL) {
92 a[X].gaussian_random();
93 }
94 b[ALL] = 0;
95 D.dagger(a, Ddaggera);
96 inverse.apply(Ddaggera, b);
97 D.apply(b, Db);
98
99 diffre = 0;
100 onsites(ALL) { diffre += squarenorm(a[X] - Db[X]); }
101 assert(diffre * diffre < 1e-16 && "test D (DdgD)^-1 Ddg");
102
103 // Run CG on a, multiply by DdaggerD and check the same
104 inverse.apply(a, b);
105 D.apply(b, Db);
106 D.dagger(Db, DdaggerDb);
107
108 diffre = 0;
109 onsites(ALL) { diffre += squarenorm(a[X] - DdaggerDb[X]); }
110 assert(diffre * diffre < 1e-16 && "test DdgD (DdgD)^-1");
111
112 // The other way around, multiply first
113 onsites(ALL) {
114 b[X].gaussian_random();
115 }
116 D.apply(b, Db);
117 D.dagger(Db, DdaggerDb);
118 a[ALL] = 0;
119 inverse.apply(DdaggerDb, a);
120
121 diffre = 0;
122 onsites(ALL) { diffre += squarenorm(a[X] - b[X]); }
123 assert(diffre * diffre < 1e-16 && "test (DdgD)^-1 DdgD");
124}
125
126// Check conjugate of the wilson Dirac operator
127{
128 hila::out0 << "Checking with Dirac_Wilson\n";
129 using dirac = Dirac_Wilson<SU<N, double>>;
130 dirac D(0.05, U);
131 Field<Wilson_vector<N, double>> a, b, Db, Ddaggera, DdaggerDb;
132#if NDIM > 3
133 a.set_boundary_condition(e_t, hila::bc::ANTIPERIODIC);
134 b.copy_boundary_condition(a);
135 Db.copy_boundary_condition(a);
136 Ddaggera.copy_boundary_condition(a);
137 DdaggerDb.copy_boundary_condition(a);
138#endif
139
140 onsites(ALL) {
141 a[X].gaussian_random();
142 b[X].gaussian_random();
143 }
144
145 double diffre = 0, diffim = 0;
146 D.apply(b, Db);
147 D.dagger(a, Ddaggera);
148 onsites(ALL) {
149 diffre += a[X].dot(Db[X]).re - Ddaggera[X].dot(b[X]).re;
150 diffim += a[X].dot(Db[X]).im - Ddaggera[X].dot(b[X]).im;
151 }
152
153 assert(diffre * diffre < 1e-16 && "test Dirac_Wilson");
154 assert(diffim * diffim < 1e-16 && "test Dirac_Wilson");
155
156 // Now run CG on Ddaggera. b=1/D a -> Db = a
157 CG<dirac> inverse(D);
158 b[ALL] = 0;
159 onsites(ALL) {
160 a[X].gaussian_random();
161 }
162 inverse.apply(a, b);
163 D.dagger(b, Db);
164 D.apply(Db, DdaggerDb);
165
166 diffre = 0;
167 onsites(ALL) { diffre += squarenorm(a[X] - DdaggerDb[X]); }
168 assert(diffre * diffre < 1e-8 && "test DdgD (DdgD)^-1");
169
170 // Now run CG on Ddaggera. b=1/D a -> Db = a
171 b[ALL] = 0;
172 onsites(ALL) {
173 a[X].gaussian_random();
174 }
175 D.dagger(a, Ddaggera);
176 inverse.apply(Ddaggera, b);
177 D.apply(b, Db);
178
179 diffre = 0;
180 onsites(ALL) { diffre += squarenorm(a[X] - Db[X]); }
181 assert(diffre * diffre < 1e-8 && "test D (DdgD)^-1 Ddg");
182}
183
184// Check conjugate of the even-odd preconditioned staggered Dirac operator
185{
186 hila::out0 << "Checking with dirac_staggered_evenodd\n";
187 dirac_staggered_evenodd D(5.0, U);
188 Field<SU_vector<N, double>> a, b, Db, Ddaggera, DdaggerDb;
189 onsites(ALL) {
190 a[X].gaussian_random();
191 b[X].gaussian_random();
192 }
193
194 double diffre = 0, diffim = 0;
195 D.apply(b, Db);
196 D.dagger(a, Ddaggera);
197 onsites(EVEN) {
198 diffre += a[X].dot(Db[X]).re - Ddaggera[X].dot(b[X]).re;
199 diffim += a[X].dot(Db[X]).im - Ddaggera[X].dot(b[X]).im;
200 }
201
202 assert(diffre * diffre < 1e-16 && "test dirac_staggered_evenodd");
203 assert(diffim * diffim < 1e-16 && "test dirac_staggered_evenodd");
204
205 // Now run CG on Ddaggera. b=1/D a -> Db = a
206 CG inverse(D);
207 b[ALL] = 0;
208 inverse.apply(Ddaggera, b);
209 D.apply(b, Db);
210
211 onsites(EVEN) { diffre += squarenorm(a[X] - Db[X]); }
212 assert(diffre * diffre < 1e-8 && "test CG");
213}
214
215// Check conjugate of the even-odd preconditioned wilson Dirac operator
216{
217 hila::out0 << "Checking with Dirac_Wilson_evenodd\n";
218 using dirac = Dirac_Wilson_evenodd<SU<N>>;
219 dirac D(0.12, U);
220 Field<Wilson_vector<N, double>> a, b, Db, Ddaggera, DdaggerDb;
221#if NDIM > 3
222 a.set_boundary_condition(e_t, hila::bc::ANTIPERIODIC);
223 b.set_boundary_condition(e_t, hila::bc::ANTIPERIODIC);
224 Db.set_boundary_condition(e_t, hila::bc::ANTIPERIODIC);
225 Ddaggera.set_boundary_condition(e_t, hila::bc::ANTIPERIODIC);
226 DdaggerDb.set_boundary_condition(e_t, hila::bc::ANTIPERIODIC);
227#endif
228
229 a[ODD] = 0;
230 b[ODD] = 0;
231 onsites(EVEN) {
232 a[X].gaussian_random();
233 b[X].gaussian_random();
234 }
235
236 double diffre = 0, diffim = 0;
237 D.apply(b, Db);
238 D.dagger(a, Ddaggera);
239 onsites(EVEN) {
240 diffre += a[X].dot(Db[X]).re - Ddaggera[X].dot(b[X]).re;
241 diffim += a[X].dot(Db[X]).im - Ddaggera[X].dot(b[X]).im;
242 }
243
244 assert(diffre * diffre < 1e-16 && "test Dirac_Wilson_evenodd");
245 assert(diffim * diffim < 1e-16 && "test Dirac_Wilson_evenodd");
246
247 // Now check that D (DdgD)^-1 Ddg = 1
248 CG<dirac> inverse(D);
249 a[ALL] = 0;
250 b[ALL] = 0;
251 onsites(EVEN) {
252 a[X].gaussian_random();
253 }
254 D.dagger(a, Ddaggera);
255 inverse.apply(Ddaggera, b);
256 D.apply(b, Db);
257
258 diffre = 0;
259 onsites(EVEN) { diffre += squarenorm(a[X] - Db[X]); }
260 hila::out0 << "D inv Dg " << diffre << "\n";
261 assert(diffre * diffre < 1e-16 && "test D (DdgD)^-1 Ddg");
262
263 // Run CG on a, multiply by DdaggerD and check the same
264 onsites(EVEN) {
265 a[X].gaussian_random();
266 }
267 b[ALL] = 0;
268 inverse.apply(a, b);
269 D.apply(b, Db);
270 D.dagger(Db, DdaggerDb);
271
272 diffre = 0;
273 onsites(EVEN) { diffre += squarenorm(a[X] - DdaggerDb[X]); }
274 assert(diffre * diffre < 1e-16 && "test DdgD (DdgD)^-1");
275
276 // The other way around, multiply first
277 onsites(EVEN) {
278 b[X].gaussian_random();
279 }
280 D.apply(b, Db);
281 D.dagger(Db, DdaggerDb);
282 a[ALL] = 0;
283 inverse.apply(DdaggerDb, a);
284
285 diffre = 0;
286 onsites(EVEN) { diffre += squarenorm(a[X] - b[X]); }
287 assert(diffre * diffre < 1e-16 && "test (DdgD)^-1 DdgD");
288}
289
290// The the Hasenbusch operator
291{
292 hila::out0 << "Checking with Hasenbusch_operator\n";
293 using dirac_base = Dirac_Wilson_evenodd<SU<N>>;
294 dirac_base Dbase(0.1, U);
295 using dirac = Hasenbusch_operator<dirac_base>;
296 dirac D(Dbase, 0.1);
297 Field<Wilson_vector<N, double>> a, b, Db, Ddaggera, DdaggerDb;
298 Field<Wilson_vector<N, double>> sol;
299
300 a[ODD] = 0;
301 b[ODD] = 0;
302 onsites(EVEN) {
303 a[X].gaussian_random();
304 b[X].gaussian_random();
305 sol[X] = 0;
306 }
307
308 double diffre = 0, diffim = 0;
309 D.apply(b, Db);
310 D.dagger(a, Ddaggera);
311 onsites(EVEN) {
312 diffre += a[X].dot(Db[X]).re - Ddaggera[X].dot(b[X]).re;
313 diffim += a[X].dot(Db[X]).im - Ddaggera[X].dot(b[X]).im;
314 }
315
316 assert(diffre * diffre < 1e-16 && "test Hasenbusch_operator");
317 assert(diffim * diffim < 1e-16 && "test Hasenbusch_operator");
318
319 // Now run CG on Ddaggera. b=1/D a -> Db = a
320 CG<dirac> inverse(D);
321 b[ALL] = 0;
322 inverse.apply(Ddaggera, b);
323 D.apply(b, Db);
324
325 onsites(EVEN) { diffre += squarenorm(a[X] - Db[X]); }
326 assert(diffre * diffre < 1e-16 && "test CG");
327
328 // Run CG on a, multiply by DdaggerD and check the same
329 onsites(EVEN) {
330 a[X].gaussian_random();
331 }
332 b[ALL] = 0;
333 inverse.apply(a, b);
334 D.apply(b, Db);
335 D.dagger(Db, DdaggerDb);
336
337 diffre = 0;
338 onsites(EVEN) { diffre += squarenorm(a[X] - DdaggerDb[X]); }
339 assert(diffre * diffre < 1e-16 && "test DdgD (DdgD)^-1");
340
341 // The other way around, multiply first
342 onsites(EVEN) {
343 b[X].gaussian_random();
344 }
345 D.apply(b, Db);
346 D.dagger(Db, DdaggerDb);
347 a[ALL] = 0;
348 inverse.apply(DdaggerDb, a);
349
350 diffre = 0;
351 onsites(EVEN) { diffre += squarenorm(a[X] - b[X]); }
352 assert(diffre * diffre < 1e-16 && "test (DdgD)^-1 DdgD");
353}
354
355hila::finishrun();
356}
squarenorm
hila::arithmetic_type< T > squarenorm(const Array< n, m, T > &rhs)
Return square norm of Array.
Definition array.h:1019
CG
The conjugate gradient operator. Applies the inverse square of an operator on a vector.
Definition conjugate_gradient.h:19
CoordinateVector_t< int >
Dirac_Wilson_evenodd
Definition wilson.h:234
Dirac_Wilson
Definition wilson.h:146
Field
The field class implements the standard methods for accessing Fields. Hilapp replaces the parity acce...
Definition field.h:62
Field::set_boundary_condition
void set_boundary_condition(Direction dir, hila::bc bc)
Set the boundary condition in a given Direction (periodic or antiperiodic)
Definition field.h:579
Field::dagger
Field< A > dagger() const
Returns dagger or Hermitian conjugate of Field depending on how it is defined for Field type T.
Definition field.h:1188
Field::copy_boundary_condition
void copy_boundary_condition(const Field< A > &rhs)
Copy the boundary condition from another field.
Definition field.h:651
Hasenbusch_operator
Definition Hasenbusch.h:15
Lattice::setup
void setup(const CoordinateVector &siz)
lattice.setup(CoordinateVector size) - set up the base lattice. Called at the beginning of the progra...
Definition lattice.h:447
Matrix_t::gaussian_random
Mtype & gaussian_random(base_type width=1.0)
Fills Matrix with gaussian random elements from gaussian distribution.
Definition matrix.h:1415
SU
Class for SU(N) matrix.
Definition sun_matrix.h:110
SU::random
const SU & random(int nhits=16)
Generate random SU(N) matrix.
Definition sun_matrix.h:211
dirac_staggered_evenodd
Definition libraries/dirac/staggered.h:193
dirac_staggered
Definition libraries/dirac/staggered.h:92
EVEN
constexpr Parity EVEN
bit pattern: 001
Definition coordinates.h:190
foralldir
#define foralldir(d)
Macro to loop over (all) Direction(s)
Definition coordinates.h:80
ODD
constexpr Parity ODD
bit pattern: 010
Definition coordinates.h:192
ALL
constexpr Parity ALL
bit pattern: 011
Definition coordinates.h:194
hila::out0
std::ostream out0
This writes output only from main process (node 0)
hila::initialize
void initialize(int argc, char **argv)
Initial setup routines.
Definition initialize.cpp:60
hila::seed_random
void seed_random(uint64_t seed, bool device_rng=true)
Seed random generators with 64-bit unsigned value. On MPI shuffles the seed so that different MPI ran...
Definition random.cpp:86
hila::finishrun
void finishrun()
Normal, controlled exit - all nodes must call this. Prints timing information and information about c...
Definition initialize.cpp:302