Bridge++  Ver. 2.0.4
afield-inc.h
Go to the documentation of this file.
1 
10 #ifndef QXS_AFIELD_INC_INCLUDED
11 #define QXS_AFIELD_INC_INCLUDED
12 
13 #include <cstdlib>
14 
15 #include "complexTraits.h"
16 #include "lib/ResourceManager/threadManager.h"
17 #include "lib_alt_QXS/inline/define_vlen.h"
18 #include "lib_alt_QXS/inline/afield_th-inc.h"
19 
20 
21 //====================================================================
22 template<typename REALTYPE>
23 void copy(AField<REALTYPE, QXS>& v, const AField<REALTYPE, QXS>& w)
24 {
25  v.copy(w);
26 }
27 
28 
29 //====================================================================
30 template<typename REALTYPE>
31 void copy(AField<REALTYPE, QXS>& v, const int ex,
32  const AField<REALTYPE, QXS>& w, const int ex_w)
33 {
34  v.copy(ex, w, ex_w);
35 }
36 
37 
38 //====================================================================
39 template<typename REALTYPE>
40 void axpy(AField<REALTYPE, QXS>& v, const int exv,
41  const typename AField<REALTYPE, QXS>::real_t a,
42  const AField<REALTYPE, QXS>& w, const int exw)
43 {
44  v.axpy(exv, a, w, exw);
45 }
46 
47 
48 //====================================================================
49 template<typename REALTYPE>
50 void axpy(AField<REALTYPE, QXS>& v,
51  const typename AField<REALTYPE, QXS>::real_t a,
52  const AField<REALTYPE, QXS>& w)
53 {
54  v.axpy(a, w);
55 }
56 
57 
58 //====================================================================
59 template<typename REALTYPE>
60 void axpy(AField<REALTYPE, QXS>& v, const typename AField<REALTYPE, QXS>::complex_t a, const AField<REALTYPE, QXS>& w)
61 {
62  v.axpy(real(a), imag(a), w);
63 }
64 
65 
66 //====================================================================
67 template<typename REALTYPE>
68 void axpy(AField<REALTYPE, QXS>& v, const int ex,
69  const typename AField<REALTYPE, QXS>::complex_t a,
70  const AField<REALTYPE, QXS>& w, const int ex_w)
71 {
72  v.axpy(ex, real(a), imag(a), w, ex_w);
73 }
74 
75 
76 //====================================================================
77 template<typename REALTYPE>
78 void aypx(const typename AField<REALTYPE, QXS>::real_t a,
79  AField<REALTYPE, QXS>& v, const AField<REALTYPE, QXS>& w)
80 {
81  v.aypx(a, w);
82 }
83 
84 
85 //====================================================================
86 template<typename REALTYPE>
87 void aypx(const typename AField<REALTYPE, QXS>::complex_t a, AField<REALTYPE, QXS>& v, const AField<REALTYPE, QXS>& w)
88 {
89  v.aypx(real(a), imag(a), w);
90 }
91 
92 
93 //====================================================================
94 template<typename REALTYPE>
95 void scal(AField<REALTYPE, QXS>& v, REALTYPE a)
96 {
97  v.scal(a);
98 }
99 
100 
101 //====================================================================
102 template<typename REALTYPE>
103 void scal(AField<REALTYPE, QXS>& v, const typename AField<REALTYPE, QXS>::complex_t a)
104 {
105  v.scal(real(a), imag(a));
106 }
107 
108 
109 //====================================================================
110 template<typename REALTYPE>
111 REALTYPE dot(AField<REALTYPE, QXS>& v, AField<REALTYPE, QXS>& w)
112 {
113  return v.dot(w);
114 }
115 
116 
117 //====================================================================
118 template<typename REALTYPE>
119 typename AField<REALTYPE, QXS>::complex_t dotc(const AField<REALTYPE, QXS>& v, const AField<REALTYPE, QXS>& w)
120 {
121  REALTYPE vw_r, vw_i;
122  v.dotc(vw_r, vw_i, w);
123  return typename AField<REALTYPE, QXS>::complex_t(vw_r, vw_i);
124 }
125 
126 
127 //====================================================================
128 template<class INDEX, class AFIELD>
129 void convert(INDEX& index, AFIELD& v, const Field& w)
130 {
131  int Nin = w.nin();
132  int Nvol = w.nvol();
133  int Nex = w.nex();
134  assert(v.check_size(Nin, Nvol, Nex));
135 
136  typename AFIELD::real_t *v2 = const_cast<AFIELD *>(&v)->ptr(0);
137 
138  int ith, nth, is, ns;
139  set_threadtask(ith, nth, is, ns, Nvol);
140 
141 #pragma omp barrier
142 
143  for (int ex = 0; ex < Nex; ++ex) {
144  for (int site = is; site < ns; ++site) {
145  for (int in = 0; in < Nin; ++in) {
146  int iw = in + Nin * (site + Nvol * ex);
147  int iv = index.idx(in, Nin, site, ex);
148  v2[iv] = w.cmp(iw);
149  }
150  }
151  }
152 
153 #pragma omp barrier
154 }
155 
156 
157 //====================================================================
158 template<class INDEX, class AFIELD>
159 void reverse(INDEX& index, Field& v, const AFIELD& w)
160 {
161  int Nin = w.nin();
162  int Nvol = w.nvol();
163  int Nex = w.nex();
164  assert(v.check_size(Nin, Nvol, Nex));
165 
166  int ith, nth, is, ns;
167  set_threadtask(ith, nth, is, ns, Nvol);
168 
169 #pragma omp barrier
170 
171  for (int ex = 0; ex < Nex; ++ex) {
172  for (int site = is; site < ns; ++site) {
173  for (int in = 0; in < Nin; ++in) {
174  int iw = in + Nin * (site + Nvol * ex);
175  int iv = index.idx(in, Nin, site, ex);
176  v.set(iv, double(w.cmp(iw)));
177  }
178  }
179  }
180 
181 #pragma omp barrier
182 }
183 
184 
185 //====================================================================
186 template<class INDEX, class FIELD>
187 void convert_spinor(INDEX& index, FIELD& v, const Field& w)
188 {
189  int Nin = w.nin();
190  int Nvol = w.nvol();
191  int Nex = w.nex();
192  assert(v.check_size(Nin, Nvol, Nex));
193 
194  int Nc = CommonParameters::Nc();
195  int Nd = CommonParameters::Nd();
196  assert(Nin == 2 * Nc * Nd);
197 
198  int ith, nth, is, ns;
199  set_threadtask(ith, nth, is, ns, Nvol);
200 
201 #pragma omp barrier
202 
203  for (int ex = 0; ex < Nex; ++ex) {
204  for (int site = is; site < ns; ++site) {
205  for (int id = 0; id < Nd; ++id) {
206  for (int ic = 0; ic < Nc; ++ic) {
207  int iwr = 2 * (ic + Nc * id) + Nin * (site + Nvol * ex);
208  int iwi = 1 + 2 * (ic + Nc * id) + Nin * (site + Nvol * ex);
209  int ivr = index.idx_SPr(ic, id, site, ex);
210  int ivi = index.idx_SPi(ic, id, site, ex);
211  v.e(ivr) = w.cmp(iwr);
212  v.e(ivi) = w.cmp(iwi);
213  }
214  }
215  }
216  }
217 
218 #pragma omp barrier
219 }
220 
221 
222 //====================================================================
223 template<class INDEX, class FIELD>
224 void convert_gauge(INDEX& index, FIELD& v, const Field& w)
225 {
226  int Nin = w.nin();
227  int Nvol = w.nvol();
228  int Nex = w.nex();
229  assert(v.check_size(Nin, Nvol, Nex));
230 
231  int Nc = CommonParameters::Nc();
232  assert(Nin == 2 * Nc * Nc);
233 
234  int ith, nth, is, ns;
235  set_threadtask(ith, nth, is, ns, Nvol);
236 
237 #pragma omp barrier
238 
239  for (int ex = 0; ex < Nex; ++ex) {
240  for (int site = is; site < ns; ++site) {
241  for (int ic2 = 0; ic2 < Nc; ++ic2) {
242  for (int ic1 = 0; ic1 < Nc; ++ic1) {
243  int iwr = 2 * (ic1 + Nc * ic2) + Nin * (site + Nvol * ex);
244  int iwi = 1 + 2 * (ic1 + Nc * ic2) + Nin * (site + Nvol * ex);
245  int ivr = index.idx_Gr(ic1, ic2, site, ex);
246  int ivi = index.idx_Gi(ic1, ic2, site, ex);
247  v.e(ivr) = w.cmp(iwr);
248  v.e(ivi) = w.cmp(iwi);
249  }
250  }
251  }
252  }
253 
254 #pragma omp barrier
255 }
256 
257 
258 //====================================================================
259 template<class INDEX2, class FIELD2, class INDEX1, class FIELD1>
260 void convert(INDEX2& index2, FIELD2& v2,
261  const INDEX1& index1, const FIELD1& v1)
262 {
263  int Nin = v1.nin();
264  int Nvol = v1.nvol();
265  int Nex = v1.nex();
266  assert(v2.check_size(Nin, Nvol, Nex));
267 
268  int ith, nth, is, ns;
269 
270  if ((sizeof(typename FIELD1::real_t) == 4) || (sizeof(typename FIELD2::real_t) == 4)) {
271  set_threadtask(ith, nth, is, ns, Nvol / VLENS);
272  for (int ex = 0; ex < Nex; ++ex) {
273  for (int vsite = is; vsite < ns; ++vsite) {
274  for (int in = 0; in < Nin; ++in) {
275  for (int vin = 0; vin < VLENS; ++vin) {
276  int site = VLENS * vsite + vin;
277  int iv1 = index1.idx(in, Nin, site, ex);
278  int iv2 = index2.idx(in, Nin, site, ex);
279  v2.e(iv2) = v1.cmp(iv1);
280  }
281  }
282  }
283  }
284  } else {
285  set_threadtask(ith, nth, is, ns, Nvol / VLEND);
286  for (int ex = 0; ex < Nex; ++ex) {
287  for (int vsite = is; vsite < ns; ++vsite) {
288  for (int in = 0; in < Nin; ++in) {
289  for (int vin = 0; vin < VLEND; ++vin) {
290  int site = VLEND * vsite + vin;
291  int iv1 = index1.idx(in, Nin, site, ex);
292  int iv2 = index2.idx(in, Nin, site, ex);
293  v2.e(iv2) = v1.cmp(iv1);
294  }
295  }
296  }
297  }
298  }
299 
300 
301 #pragma omp barrier
302 }
303 
304 
305 //====================================================================
306 template<class INDEX2, class FIELD2, class INDEX1, class FIELD1>
307 void convert_h(INDEX2& index2, FIELD2& v2,
308  const INDEX1& index1, const FIELD1& v1)
309 {
310  int Nin = v1.nin();
311  int Nvol = v1.nvol();
312  int Nex = v1.nex();
313  assert(v2.check_size(Nin, Nvol, Nex));
314 
315  int ith, nth, is, ns;
316  if ((sizeof(typename FIELD1::real_t) == 4) || (sizeof(typename FIELD2::real_t) == 4)) {
317  set_threadtask(ith, nth, is, ns, Nvol / VLENS);
318  for (int ex = 0; ex < Nex; ++ex) {
319  for (int vsite = is; vsite < ns; ++vsite) {
320  for (int in = 0; in < Nin; ++in) {
321  for (int vin = 0; vin < VLENS; ++vin) {
322  int site = VLENS * vsite + vin;
323  int iv1 = index1.idxh(in, Nin, site, ex);
324  int iv2 = index2.idxh(in, Nin, site, ex);
325  v2.e(iv2) = v1.cmp(iv1);
326  }
327  }
328  }
329  }
330  } else {
331  set_threadtask(ith, nth, is, ns, Nvol / VLEND);
332  for (int ex = 0; ex < Nex; ++ex) {
333  for (int vsite = is; vsite < ns; ++vsite) {
334  for (int in = 0; in < Nin; ++in) {
335  for (int vin = 0; vin < VLEND; ++vin) {
336  int site = VLEND * vsite + vin;
337  int iv1 = index1.idxh(in, Nin, site, ex);
338  int iv2 = index2.idxh(in, Nin, site, ex);
339  v2.e(iv2) = v1.cmp(iv1);
340  }
341  }
342  }
343  }
344  }
345 
346 #pragma omp barrier
347 }
348 
349 
350 //====================================================================
351 template<class INDEX, class FIELD>
352 void reverse(INDEX& index, Field& v, FIELD& w)
353 {
354  int Nin = v.nin();
355  int Nvol = v.nvol();
356  int Nex = v.nex();
357  w.check_size(Nin, Nvol, Nex);
358 
359  int ith, nth, is, ns;
360  set_threadtask(ith, nth, is, ns, Nvol);
361 
362 #pragma omp barrier
363 
364  for (int ex = 0; ex < Nex; ++ex) {
365  for (int site = is; site < ns; ++site) {
366  for (int in = 0; in < Nin; ++in) {
367  int iv = in + Nin * (site + Nvol * ex);
368  int iw = index.idx(in, Nin, site, ex);
369  v.set(iv, double(w.cmp(iw)));
370  }
371  }
372  }
373 
374 #pragma omp barrier
375 }
376 
377 
378 //====================================================================
379 template<class INDEX, class FIELD>
380 void reverse_spinor(INDEX& index, Field& v, FIELD& w)
381 {
382  int Nin = v.nin();
383  int Nvol = v.nvol();
384  int Nex = v.nex();
385  w.check_size(Nin, Nvol, Nex);
386 
387  int Nc = CommonParameters::Nc();
388  int Nd = CommonParameters::Nd();
389  assert(Nin == 2 * Nc * Nd);
390 
391  int ith, nth, is, ns;
392  set_threadtask(ith, nth, is, ns, Nvol);
393 
394 #pragma omp barrier
395 
396  for (int ex = 0; ex < Nex; ++ex) {
397  for (int site = is; site < ns; ++site) {
398  for (int id = 0; id < Nd; ++id) {
399  for (int ic = 0; ic < Nc; ++ic) {
400  int ivr = 2 * (ic + Nc * id) + Nin * (site + Nvol * ex);
401  int ivi = 1 + 2 * (ic + Nc * id) + Nin * (site + Nvol * ex);
402  int iwr = index.idx_SPr(ic, id, site, ex);
403  int iwi = index.idx_SPi(ic, id, site, ex);
404  v.set(ivr, double(w.cmp(iwr)));
405  v.set(ivi, double(w.cmp(iwi)));
406  }
407  }
408  }
409  }
410 
411 #pragma omp barrier
412 }
413 
414 
415 //====================================================================
416 template<class INDEX, class FIELD>
417 void reverse_gauge(INDEX& index, Field& v, FIELD& w)
418 {
419  int Nin = v.nin();
420  int Nvol = v.nvol();
421  int Nex = v.nex();
422  w.check_size(Nin, Nvol, Nex);
423 
424  int Nc = CommonParameters::Nc();
425  assert(Nin == 2 * Nc * Nc);
426 
427  int ith, nth, is, ns;
428  set_threadtask(ith, nth, is, ns, Nvol);
429 
430 #pragma omp barrier
431 
432  for (int ex = 0; ex < Nex; ++ex) {
433  for (int site = is; site < ns; ++site) {
434  for (int ic2 = 0; ic2 < Nc; ++ic2) {
435  for (int ic1 = 0; ic1 < Nc; ++ic1) {
436  int ivr = 2 * (ic1 + Nc * ic2) + Nin * (site + Nvol * ex);
437  int ivi = 1 + 2 * (ic1 + Nc * ic2) + Nin * (site + Nvol * ex);
438  int iwr = index.idx_Gr(ic1, ic2, site, ex);
439  int iwi = index.idx_Gi(ic1, ic2, site, ex);
440  v.set(ivr, double(w.cmp(iwr)));
441  v.set(ivi, double(w.cmp(iwi)));
442  }
443  }
444  }
445  }
446 
447 #pragma omp barrier
448 }
449 
450 
451 //====================================================================
452 template<typename REALTYPE>
453 REALTYPE norm2(const AField<REALTYPE, QXS>& v)
454 {
455  AField<REALTYPE, QXS> *vp = const_cast<AField<REALTYPE, QXS> *>(&v);
456  return vp->norm2();
457  // return v.norm2();
458 }
459 
460 
461 //====================================================================
462 template<typename REALTYPE>
463 void dotc_and_norm2(typename AField<REALTYPE, QXS>::complex_t& dotc,
464  REALTYPE& v_norm2, REALTYPE& w_norm2,
465  const AField<REALTYPE, QXS>& v,
466  const AField<REALTYPE, QXS>& w)
467 {
468  // returns <v|w>, <v|v>, <w|w>
469  assert(v.nex() == w.nex());
470  assert(v.nvol() == w.nvol());
471  assert(v.nin() == w.nin());
472 
473  dotc = v.dotc_and_norm2(v_norm2, w_norm2, w);
474 }
475 
476 
477 #define AFIELD_HAS_DOTC_AND_NORM2
478 
479 //============================================================END=====
480 #endif
scal
void scal(AField< REALTYPE, QXS > &v, REALTYPE a)
Definition: afield-inc.h:95
AField< REALTYPE, QXS >::nex
int nex() const
returning size of extra d.o.f.
Definition: afield.h:115
AField< REALTYPE, QXS >::axpy
void axpy(const real_t, const AField< real_t, QXS > &)
this is to be discarded.
Definition: afield-tmpl.h:154
VLEND
#define VLEND
Definition: define_vlen.h:42
AField< REALTYPE, QXS >::dot
real_t dot(const AField< real_t, QXS > &) const
Definition: afield-tmpl.h:523
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
AField< REALTYPE, QXS >
Definition: afield.h:35
convert
void convert(INDEX &index, AFIELD &v, const Field &w)
Definition: afield-inc.h:129
AField
Definition: afield_base.h:16
AField< REALTYPE, QXS >::nvol
int nvol() const
returning size of site d.o.f.
Definition: afield.h:112
Field::nex
int nex() const
Definition: field.h:128
Field::check_size
bool check_size(const int nin, const int nvol, const int nex) const
checking size parameters. [23 May 2016 H.Matsufuru]
Definition: field.h:135
convert_gauge
void convert_gauge(INDEX &index, FIELD &v, const Field &w)
Definition: afield-inc.h:224
reverse_gauge
void reverse_gauge(INDEX &index, Field &v, FIELD &w)
Definition: afield-inc.h:417
AField< REALTYPE, QXS >::norm2
real_t norm2(void) const
Definition: afield-tmpl.h:656
VLENS
#define VLENS
Definition: define_vlen.h:41
AField< REALTYPE, QXS >::copy
void copy(const Field &w)
Definition: afield-tmpl.h:71
real_t
double real_t
Definition: bridgeQXS_Clover_coarse_double.cpp:16
norm2
REALTYPE norm2(const AField< REALTYPE, QXS > &v)
Definition: afield-inc.h:453
Field::nin
int nin() const
Definition: field.h:126
Field::real_t
double real_t
Definition: field.h:51
dot
REALTYPE dot(AField< REALTYPE, QXS > &v, AField< REALTYPE, QXS > &w)
Definition: afield-inc.h:111
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
axpy
void axpy(AField< REALTYPE, QXS > &v, const int exv, const typename AField< REALTYPE, QXS >::real_t a, const AField< REALTYPE, QXS > &w, const int exw)
Definition: afield-inc.h:40
reverse_spinor
void reverse_spinor(INDEX &index, Field &v, FIELD &w)
Definition: afield-inc.h:380
threadManager.h
dotc
AField< REALTYPE, QXS >::complex_t dotc(const AField< REALTYPE, QXS > &v, const AField< REALTYPE, QXS > &w)
Definition: afield-inc.h:119
Field::nvol
int nvol() const
Definition: field.h:127
copy
void copy(AField< REALTYPE, QXS > &v, const AField< REALTYPE, QXS > &w)
Definition: afield-inc.h:23
AField< REALTYPE, QXS >::scal
void scal(const real_t)
Definition: afield-tmpl.h:453
AField< REALTYPE, QXS >::dotc
void dotc(real_t &, real_t &, const AField< real_t, QXS > &) const
Definition: afield-tmpl.h:579
Field::cmp
double cmp(const int jin, const int site, const int jex) const
Definition: field.h:143
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
reverse
void reverse(INDEX &index, Field &v, const AFIELD &w)
Definition: afield-inc.h:159
aypx
void aypx(const typename AField< REALTYPE, QXS >::real_t a, AField< REALTYPE, QXS > &v, const AField< REALTYPE, QXS > &w)
Definition: afield-inc.h:78
complex_t
ComplexTraits< double >::complex_t complex_t
Definition: afopr_Clover_coarse_double.cpp:23
complexTraits.h
Field
Container of Field-type object.
Definition: field.h:46
dotc_and_norm2
void dotc_and_norm2(typename AField< REALTYPE, QXS >::complex_t &dotc, REALTYPE &v_norm2, REALTYPE &w_norm2, const AField< REALTYPE, QXS > &v, const AField< REALTYPE, QXS > &w)
Definition: afield-inc.h:463
convert_h
void convert_h(INDEX2 &index2, FIELD2 &v2, const INDEX1 &index1, const FIELD1 &v1)
Definition: afield-inc.h:307
AField< REALTYPE, QXS >::aypx
void aypx(const int ex, const real_t ar, const real_t ai, const AField< real_t, QXS > &w, const int ex_w)
Definition: afield-tmpl.h:401
AField< REALTYPE, QXS >::nin
int nin() const
returning size of inner (on site) d.o.f.
Definition: afield.h:109
AField< REALTYPE, QXS >::dotc_and_norm2
complex_t dotc_and_norm2(real_t &, real_t &, const AField< real_t, QXS > &) const
Definition: afield-tmpl.h:707
convert_spinor
void convert_spinor(INDEX &index, FIELD &v, const Field &w)
Definition: afield-inc.h:187