HILA
Loading...
Searching...
No Matches
bulk_prevention_action.h
Go to the documentation of this file.
1/** @file bulk_prevention_action.h */
2
3#ifndef BULK_PREVENTION_ACTION_H_
4#define BULK_PREVENTION_ACTION_H_
5
6#include "hila.h"
7
8// functions for bulk-preventing action, cf. arXiv:2306.14319 (with n=2)
9template <typename T>
10T ch_inv(const T &U) {
11 // compute inverse of the square matrix U, using the
12 // Cayley-Hamilton theorem (Faddeev-LeVerrier algorithm)
13 T tB[2];
14 int ip = 0;
15 tB[ip] = 1.;
16 auto tc = trace(U);
17 tB[1 - ip] = U;
18 for (int k = 2; k <= T::size(); ++k) {
19 tB[1 - ip] -= tc;
20 mult(U, tB[1 - ip], tB[ip]);
21 tc = trace(tB[ip]) / k;
22 ip = 1 - ip;
23 }
24 return tB[ip] / tc;
25}
26
27template <typename T>
28T bp_UAmat(const T &U) {
29 // compute U*A(U) with the A-matrix from Eq. (B3) of arXiv:2306.14319 for n=2
30 T tA1 = U;
31 tA1 += 1.;
32 tA1 *= 0.5;
33 T tA2 = ch_inv(tA1);
34 mult(tA2, tA2.dagger(), tA1);
35 return U * tA1 * tA1 * tA2;
36}
37
38template <typename T>
39T bp_iOsqmat(const T &U) {
40 // compute matrix inside the trace on r.h.s. of Eq. (B1) of arXiv:2306.14319 for n=2
41 // (without identiy matrix subtraction)
42 T tA1 = U;
43 tA1 += 1.;
44 tA1 *= 0.5;
45 T tA2 = ch_inv(tA1);
46 mult(tA2, tA2.dagger(), tA1);
47 return tA1 * tA1;
48}
49
50template <typename group>
51double measure_s_bp(const GaugeField<group> &U) {
52 // measure the BP plaquette action
53 Reduction<double> plaq = 0;
54 plaq.allreduce(false).delayed(true);
55 foralldir(dir1) foralldir(dir2) if (dir1 < dir2) {
56 U[dir2].start_gather(dir1, ALL);
57 U[dir1].start_gather(dir2, ALL);
58 onsites(ALL) {
59 plaq += real(trace(bp_iOsqmat(U[dir1][X] * U[dir2][X + dir1] *
60 (U[dir2][X] * U[dir1][X + dir2]).dagger()))) /
61 group::size() -
62 1.0;
63 }
64 }
65 return plaq.value();
66}
67
68template <typename group, typename atype = hila::arithmetic_type<group>>
69void get_force_bp_add(const GaugeField<group> &U, VectorField<Algebra<group>> &K, atype eps = 1.0) {
70 // compute force for BP action for n=2 according to Eq. (B5) of arXiv:2306.14319 and add it to K
71 Field<group> fmatp;
72 Field<group> fmatmd1;
73 Field<group> fmatmd2;
74 atype teps = 2.0 * eps;
75 foralldir(d1) {
76 foralldir(d2) if (d1 < d2) {
77 U[d2].start_gather(d1, ALL);
78 U[d1].start_gather(d2, ALL);
79 onsites(ALL) {
80 fmatp[X] = bp_UAmat(U[d1][X] * U[d2][X + d1] * (U[d2][X] * U[d1][X + d2]).dagger());
81 fmatmd1[X] = (fmatp[X] * U[d2][X]).dagger() *
82 U[d2][X]; // parallel transport fmatp[X].dagger() to X+d2
83 fmatmd2[X] =
84 U[d1][X].dagger() * fmatp[X] * U[d1][X]; // parallel transport fmatp[X] to X+d1
85 }
86 fmatmd1.start_gather(-d2, ALL);
87 fmatmd2.start_gather(-d1, ALL);
88 onsites(ALL) {
89 K[d1][X] -= (fmatmd1[X - d2] + fmatp[X]).project_to_algebra_scaled(teps);
90 K[d2][X] -= (fmatmd2[X - d1] - fmatp[X]).project_to_algebra_scaled(teps);
91 }
92 }
93 }
94}
95
96template <typename group, typename atype = hila::arithmetic_type<group>>
97void get_force_bp(const GaugeField<group> &U, out_only VectorField<Algebra<group>> &K,
98 atype eps = 1.0) {
99 // compute force for BP action for n=2 according to Eq. (B5) of arXiv:2306.14319 and write it to
100 // K
101 foralldir(d1) {
102 K[d1][ALL] = 0;
103 }
104 get_force_bp_add(U, K, eps);
105}
106
107/*
108template <typename group, typename atype = hila::arithmetic_type<group>>
109void get_force_bp(const GaugeField<group> &U, out_only VectorField<group> &K) {
110 // compute force matrix for BP action for n=2 according to Eq. (B5) of arXiv:2306.14319 and add it to K
111 Field<group> fmatp;
112 Field<group> fmatmd1;
113 Field<group> fmatmd2;
114 foralldir(d1) {
115 foralldir(d2) if (d1 < d2) {
116 U[d2].start_gather(d1, ALL);
117 U[d1].start_gather(d2, ALL);
118 }
119 }
120 foralldir(d1) {
121 K[d1][ALL] = 0;
122 }
123 foralldir(d1) {
124 foralldir(d2) if (d1 < d2) {
125 onsites(ALL) {
126 fmatp[X] = bp_UAmat(U[d1][X] * U[d2][X + d1] * (U[d2][X] * U[d1][X + d2]).dagger());
127 fmatp[X] -= fmatp[X].dagger();
128 fmatmd1[X] =
129 U[d2][X].dagger() * fmatp[X] * U[d2][X]; // parallel transport fmatp[X] to X+d2
130 fmatmd2[X] =
131 U[d1][X].dagger() * fmatp[X] * U[d1][X]; // parallel transport fmatp[X] to X+d1
132 }
133 fmatmd1.start_gather(-d2, ALL);
134 fmatmd2.start_gather(-d1, ALL);
135 onsites(ALL) {
136 K[d1][X] -= (fmatp[X] - fmatmd1[X - d2]);
137 K[d2][X] -= (fmatmd2[X - d1] - fmatp[X]);
138 }
139 }
140 onsites(ALL) {
141 K[d1][X] -= K[d1][X].dagger();
142 K[d1][X] *= 0.5;
143 }
144 }
145}
146
147template <typename group, typename atype = hila::arithmetic_type<group>>
148void get_force_bp_add2(const GaugeField<group> &U, VectorField<Algebra<group>> &K, atype eps = 1.0) {
149 VectorField<group> tK;
150 get_force_bp(U, tK);
151 foralldir(d1) {
152 onsites(ALL) K[d1][X] += tK[d1][X].project_to_algebra_scaled(eps);
153 }
154}
155*/
156
157#endif
Array< n, m, hila::arithmetic_type< T > > real(const Array< n, m, T > &arg)
Return real part of Array.
Definition array.h:688
Definition su2.h:7
The field class implements the standard methods for accessing Fields. Hilapp replaces the parity acce...
Definition field.h:61
dir_mask_t start_gather(Direction d, Parity p=ALL) const
Definition field_comm.h:262
Gauge field class.
Definition gaugefield.h:22
Special reduction class: enables delayed and non-blocking reductions, which are not possible with the...
Definition reduction.h:14
const T value()
Return value of the reduction variable. Wait for the comms if needed.
Definition reduction.h:157
Reduction & allreduce(bool b=true)
allreduce(bool) turns allreduce on or off. By default on.
Definition reduction.h:130
Reduction & delayed(bool b=true)
deferred(bool) turns deferred on or off. By default turns on.
Definition reduction.h:148
Complex< T > dagger(const Complex< T > &val)
Return dagger of Complex number.
Definition cmplx.h:1223
#define foralldir(d)
Macro to loop over (all) Direction(s)
Definition coordinates.h:78
constexpr Parity ALL
bit pattern: 011
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