CHROMA
expsu3.cc
Go to the documentation of this file.
1 /*# $Id: expsu3.cc,v 3.1 2006-08-25 23:46:37 edwards Exp $ ($Date: 2006-08-25 23:46:37 $) */
2 /*# EXPSU3: */
3 /*# In place a = 1 + a + (1/2)*a^2 + ...+ (1/n!)*(a)^n n = 12 */
4 /*# a must be a lattice_complex32_colour_colour primitive */
5 
6 /*# If it does not satisfy the unitary condition, we eesu3 () to invoke */
7 /*# the exact exponentiator. */
8 /*# cflag decides whether we do the checking */
9 #include "chromabase.h"
10 
11 #include "util/gauge/expsu3.h"
12 
13 
14 namespace Chroma {
15 
16 void expsu3(LatticeColorMatrix& a, Reunitarize cflag)
17 {
18  START_CODE();
19 
20  LatticeColorMatrix aux1;
21  LatticeColorMatrix aux2;
22  LatticeColorMatrix aux3;
23  LatticeReal aa;
24  LatticeReal bb;
25  LatticeReal te1;
26  LatticeReal te2;
27  LatticeReal te3;
28  LatticeReal tmp;
29  LatticeReal x0_real;
30  LatticeReal x0_imag;
31  LatticeReal x1_real;
32  LatticeReal x1_imag;
33  LatticeReal x2_real;
34  LatticeReal x2_imag;
35  LatticeComplex x0;
36  LatticeComplex x1;
37  LatticeComplex x2;
38  LatticeBoolean bad;
39  Real ftmp1;
40  Real ftmp2;
41 
42  int numbad;
43  int i;
44 
45  if ( Nc != 3 )
46  QDP_error_exit("can only handle SU[3]", Nc);
47 
48 
49  /*# aux1 = a*a */
50  aux1 = a * a;
51 
52  /*# aa = -tr(a*a)/2 */
53  aa = -Real(0.5)*real(trace(aux1));
54 
55  /*# aux2 = a*a*a */
56  aux2 = aux1 * a;
57 
58  /*# bb = -Im( tr(a*a*a)/3 ) */
59  bb = -imag(trace(aux2))/Real(3);
60 
61  ftmp1 = Real(4);
62  ftmp2 = Real(396);
63  te1 = aa * ftmp1; /* multiply by 4.0 */
64  tmp=ftmp2;
65  te1 -= tmp; /* subtract 396.0 */
66 
67  ftmp2 = Real(23760);
68  te2 = aa * te1; /* multiply by te1 */
69  tmp=ftmp2;
70  te2 += tmp; /* add 23760.0 */
71 
72  ftmp2 = Real(665280);
73  te1 = aa * te2; /* multiply by te1 */
74  tmp=ftmp2;
75  te1 -= tmp; /* subtract 665280.0 */
76 
77  te2 = bb * bb; /* multiply by bb */
78 
79  te3 = te1;
80  te3 += te2; /* te3 = te1 + te2 */
81 
82  ftmp2 = Real(479001600);
83  te1 = te2 * te3; /* multiply by te1 */
84  tmp=ftmp2;
85  te1 += tmp; /* add 479001600.0 */
86 
87  ftmp1 = Real(2.08767569878681E-09);
88  x0_real = te1 * ftmp1; /* multiply by 2.08767569878681e-09 */
89 
90 
91 
92  ftmp1 = -Real(3);
93  ftmp2 = Real(110);
94  te1 = aa * ftmp1; /* multiply by -3.0 */
95  tmp=ftmp2;
96  te1 += tmp; /* add 110.0.0 */
97 
98  x0_imag = te1 * te2; /* multiply by te2 */
99 
100  tmp=ftmp2;
101  te1 = tmp; /* te1 = 110.0 - aa */
102  te1 -= aa;
103 
104  ftmp2 = Real(7920);
105  te3 = aa * te1; /* multiply by te1 */
106  tmp=ftmp2;
107  te3 -= tmp; /* add 7920.0 */
108 
109  ftmp2 = Real(332640);
110  te1 = aa * te3; /* multiply by te3 */
111  tmp=ftmp2;
112  te1 += tmp; /* add 332640.0 */
113 
114  ftmp2 = Real(6652800);
115  te3 = aa * te1; /* multiply by te3 */
116  tmp=ftmp2;
117  te3 -= tmp; /* subtract 6652800.0 */
118 
119  te1 = te3;
120  te1 += x0_imag; /* te1 = te3 + x0_imag */
121 
122  te3 = bb * te1; /* multiply by te1 */
123 
124  ftmp1 = Real(2.505210838544172E-08);
125  x0_imag = te3 * ftmp1; /* multiply by 2.505210838544172e-08 */
126 
127 
128 
129  ftmp1 = -Real(6);
130  ftmp2 = Real(330);
131  te1 = aa * ftmp1; /* multiply by -6.0 */
132  tmp=ftmp2;
133  te1 += tmp; /* add 330.0 */
134 
135  ftmp2 = Real(7920);
136  te3 = aa * te1; /* multiply by te1 */
137  tmp=ftmp2;
138  te3 -= tmp; /* subtract 7920.0 */
139 
140  x1_real = te3 * te2; /* x1_real = te3 * te2 */
141 
142  ftmp2 = Real(110);
143  tmp=ftmp2;
144  te1 = tmp; /* te1 = 110.0 - aa */
145  te1 -= aa;
146 
147  ftmp2 = Real(7920);
148  te3 = aa * te1; /* multiply by te1 */
149  tmp=ftmp2;
150  te3 -= tmp; /* subtract 7920.0 */
151 
152  ftmp2 = Real(332640);
153  te1 = aa * te3; /* multiply by te1 */
154  tmp=ftmp2;
155  te1 += tmp; /* add 332640.0 */
156 
157  ftmp2 = Real(6652800);
158  te3 = aa * te1; /* multiply by te1 */
159  tmp=ftmp2;
160  te3 -= tmp; /* subtract 6652800.0 */
161 
162  ftmp2 = Real(39916800);
163  te1 = aa * te3; /* multiply by te3 */
164  tmp=ftmp2;
165  te1 += tmp; /* add 39916800.0 */
166 
167  x1_real += te1; /* x1_real += te1 */
168 
169  ftmp1 = Real(2.505210838544172E-08);
170  x1_real = x1_real * ftmp1; /* multiply by 2.505210838544172e-08 */
171 
172 
173 
174  ftmp1 = -Real(4);
175  ftmp2 = Real(132);
176  te1 = aa * ftmp1; /* multiply by -4.0 */
177  tmp=ftmp2;
178  te1 += tmp; /* add 132.0 */
179 
180  x1_imag = te1 * te2; /* x1_imag = te1 * te2 */
181 
182  ftmp1 = -Real(5);
183  ftmp2 = Real(528);
184  te1 = aa * ftmp1; /* multiply by -5.0 */
185  tmp=ftmp2;
186  te1 += tmp; /* add 528.0 */
187 
188  ftmp2 = Real(35640);
189  te3 = aa * te1; /* multiply by te1 */
190  tmp=ftmp2;
191  te3 -= tmp; /* subtract 35640.0 */
192 
193  ftmp2 = Real(1330560);
194  te1 = aa * te3; /* multiply by te1 */
195  tmp=ftmp2;
196  te1 += tmp; /* add 1330560.0 */
197 
198  ftmp2 = Real(19958400);
199  te3 = aa * te1; /* multiply by te1 */
200  tmp=ftmp2;
201  te3 -= tmp; /* subtract 19958400.0 */
202 
203  te1 = x1_imag; /* te1 = x1_imag + te3 */
204  te1 += te3;
205 
206  ftmp1 = Real(2.08767569878681E-09);
207  x1_imag = bb * te1; /* x1_imag = bb * te1 */
208  x1_imag = x1_imag * ftmp1; /* x1_imag *= 2.08767569878681e-09 */
209 
210 
211 
212  ftmp1 = -Real(6);
213  ftmp2 = Real(396);
214  te1 = aa * ftmp1; /* multiply by -6.0 */
215  tmp=ftmp2;
216  te1 += tmp; /* add 396.0 */
217 
218  ftmp2 = Real(11880);
219  te3 = aa * te1; /* multiply by te1 */
220  tmp=ftmp2;
221  te3 -= tmp; /* subtract 11880.0 */
222 
223  x2_real = te3 * te2; /* x2_real = te2 * te3 */
224 
225  ftmp2 = Real(132);
226  te1=ftmp2; /* te1 = 132.0 - aa */
227  te1 -= aa;
228 
229  ftmp2 = Real(11880);
230  te3 = aa * te1; /* multiply by te1 */
231  tmp=ftmp2;
232  te3 -= tmp; /* subtract 11880.0 */
233 
234  ftmp2 = Real(665280);
235  te1 = aa * te3; /* multiply by te3 */
236  tmp=ftmp2;
237  te1 += tmp; /* add 665280.0 */
238 
239  ftmp2 = Real(19958400);
240  te3 = aa * te1; /* multiply by te1 */
241  tmp=ftmp2;
242  te3 -= tmp; /* subtract 19958400.0 */
243 
244  ftmp2 = Real(239500800);
245  te1 = aa * te3; /* multiply by te3 */
246  tmp=ftmp2;
247  te1 += tmp; /* add 239500800.0 */
248 
249  ftmp1 = Real(2.08767569878681E-09);
250  x2_real += te1; /* x2_real += te1 */
251  x2_real = ftmp1 * x2_real; /* x2_real *= 2.08767569878681e-09 */
252 
253 
254 
255  ftmp1 = Real(4);
256  ftmp2 = Real(330);
257  te1 = ftmp1 * aa; /* multiply by 4.0 */
258  tmp=ftmp2;
259  te1 -= tmp; /* subtract 330.0 */
260 
261  ftmp2 = Real(15840);
262  te3 = aa * te1; /* multiply by te1 */
263  tmp=ftmp2;
264  te3 += tmp; /* add 15840.0 */
265 
266  ftmp2 = Real(332640);
267  te1 = aa * te3; /* multiply by te3 */
268  tmp=ftmp2;
269  te1 -= tmp; /* subtract 332640.0 */
270 
271  te3 = te1; /* te3 = te1 + te2 */
272  te3 += te2;
273 
274  ftmp1 = Real(2.505210838544172E-08);
275  x2_imag = bb * te3; /* x2_imag = bb * te3 */
276  x2_imag = ftmp1 * x2_imag; /* x2_imag *= 2.505210838544172e-08 */
277 
278 
279  /*# Build lattice complex variables out of two lattice real variables. */
280  x0 = cmplx(x0_real,x0_imag);
281  x1 = cmplx(x1_real,x1_imag);
282  x2 = cmplx(x2_real,x2_imag);
283 
284 
285  aux2 = a * x1;
286  aux2 += aux1 * x2;
287 
288  aux3 = 1;
289  aux2 += x0 * aux3;
290 
291  if ( cflag == REUNITARIZE_LABEL ) {
292  reunit (aux2, bad, numbad, REUNITARIZE_LABEL);
293  if ( numbad > 0 ) {
294 
295  QDP_error_exit("found some matrices violating unitarity", numbad);
296  /*# eesun (a, aux1, aa, bb, aux2); */
297  }
298  }
299  a = aux2;
300 
301  END_CODE();
302 }
303 
304 
305 } // End namespace Chroma
chromabase.h
Primary include file for CHROMA library code.
numbad
int numbad
Definition: cool.cc:28
bad
LatticeBoolean bad
Definition: cool.cc:22
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::QDP_error_exit
QDP_error_exit("too many BiCG iterations", n_count, rsd_sq, cp, c, re_rvr, im_rvr, re_a, im_a, re_b, im_b)
Chroma::tmp
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
Chroma::i
int i
Definition: pbg5p_w.cc:55
Chroma::expsu3
void expsu3(LatticeColorMatrix &a, Reunitarize cflag)
Definition: expsu3.cc:16
Chroma::reunit
void reunit(LatticeColorMatrixF3 &xa)
Definition: reunit.cc:467
Chroma::a
Complex a
Definition: invbicg.cc:95
Chroma::Reunitarize
Reunitarize
Definition: reunit.h:29
Chroma::REUNITARIZE_LABEL
@ REUNITARIZE_LABEL
Definition: reunit.h:29
Chroma::START_CODE
START_CODE()
ftmp1
Double ftmp1
Definition: topol.cc:29
ftmp2
Double ftmp2
Definition: topol.cc:30
Generated by doxygen 1.9.1