Bridge++  Ver.2.1.3
field.cpp
Go to the documentation of this file.
1 
14 #include "Field/field.h"
15 
16 #include <cstring>
17 using std::string;
18 
19 #include "bridge_defs.h"
20 #include "lib/ResourceManager/threadManager.h"
21 #include "lib/IO/bridgeIO.h"
22 using Bridge::vout;
23 
24 #include "Field/field_thread-inc.h"
25 
26 const std::string Field::class_name = "Field";
27 
28 //====================================================================
29 void Field::check()
30 {
31  // ThreadManager::assert_single_thread(class_name);
32  // vout.general("Field was constructed.\n");
33 }
34 
35 
36 //====================================================================
37 void Field::set(double a)
38 {
39  int ith, nth, is, ns;
40  set_threadtask(ith, nth, is, ns, m_Nvol);
41 
42  double *RESTRICT yp = this->ptr(0);
43 
44  for (int ex = 0; ex < m_Nex; ++ex) {
45  for (int site = is; site < ns; ++site) {
46  int kv = m_Nin * (site + m_Nvol * ex);
47  for (int in = 0; in < m_Nin; ++in) {
48  yp[in + kv] = a;
49  }
50  }
51  }
52 }
53 
54 
55 //====================================================================
56 void Field::setc(const double a)
57 {
58  if (m_element_type == Element_type::REAL) {
59  set(a);
60  } else if (m_element_type == Element_type::COMPLEX) {
61  dcomplex c = cmplx(a, 0.0);
62  setc(c);
63  } else {
64  vout.crucial("Error at %s: unsupported arg types.\n", __func__);
65  exit(EXIT_FAILURE);
66  }
67 }
68 
69 
70 //====================================================================
71 void Field::setc(const dcomplex a)
72 {
73  if (imag(a) == 0.0) {
74  return set(real(a));
75  } else if (m_element_type == Element_type::COMPLEX) {
76  double *RESTRICT yp = this->ptr(0);
77  double ar = real(a);
78  double ai = imag(a);
79  int Nin2 = m_Nin / 2;
80 
81  int ith, nth, is, ns;
82  set_threadtask(ith, nth, is, ns, m_Nvol);
83 
84  for (int ex = 0; ex < m_Nex; ++ex) {
85  for (int site = is; site < ns; ++site) {
86  int kv = m_Nin * (site + m_Nvol * ex);
87  for (int k = 0; k < Nin2; ++k) {
88  int kr = 2 * k;
89  int ki = 2 * k + 1;
90  yp[kr + kv] = ar;
91  yp[ki + kv] = ai;
92  }
93  }
94  }
95  } else if (m_element_type == Element_type::REAL) {
96  double ar = real(a);
97  double ai = imag(a);
98 
99  if (fabs(ai) < fabs(ar) * CommonParameters::epsilon_criterion()) {
100  return set(real(a));
101  } else {
102  vout.crucial("Error at %s: real vector and complex parameter.\n",
103  __func__);
104  exit(EXIT_FAILURE);
105  }
106  } else {
107  vout.crucial("Error at %s: unsupported arg types.\n", __func__);
108  exit(EXIT_FAILURE);
109  }
110 }
111 
112 
113 //====================================================================
114 double Field::norm2() const
115 {
116  const double *RESTRICT yp = this->ptr(0);
117 
118  int ith, nth, is, ns;
119  set_threadtask(ith, nth, is, ns, m_Nvol);
120 
121  double a = 0.0;
122 
123  for (int ex = 0; ex < m_Nex; ++ex) {
124  for (int site = is; site < ns; ++site) {
125  int kv = m_Nin * (site + m_Nvol * ex);
126  for (int in = 0; in < m_Nin; ++in) {
127  a += yp[in + kv] * yp[in + kv];
128  }
129  }
130  }
131 
132  ThreadManager::reduce_sum_global(a, ith, nth);
133 
134  return a;
135 }
136 
137 
138 //====================================================================
139 void Field::xI()
140 {
141  if (m_element_type == Element_type::REAL) {
142  vout.crucial("Error at %s: xI() is not available for real field.\n",
143  class_name.c_str());
144  exit(EXIT_FAILURE);
145  }
146 
147  double *RESTRICT yp = this->ptr(0);
148  int Nin2 = m_Nin / 2;
149 
150  int ith, nth, is, ns;
151  set_threadtask(ith, nth, is, ns, m_Nvol);
152 
153  for (int ex = 0; ex < m_Nex; ++ex) {
154  for (int site = is; site < ns; ++site) {
155  int kv = m_Nin * (site + m_Nvol * ex);
156  for (int k = 0; k < Nin2; ++k) {
157  int kr = 2 * k;
158  int ki = 2 * k + 1;
159  double yr = yp[kr + kv];
160  double yi = yp[ki + kv];
161  yp[kr + kv] = -yi;
162  yp[ki + kv] = yr;
163  }
164  }
165  }
166 }
167 
168 
169 //====================================================================
170 void Field::stat(double& Fave, double& Fmax, double& Fdev) const
171 {
172 #pragma omp barrier
173 
174  double sum = 0.0;
175  double sum2 = 0.0;
176  double vmax = 0.0;
177 
178  const double *RESTRICT yp = this->ptr(0);
179 
180  int ith, nth, is, ns;
181  set_threadtask(ith, nth, is, ns, m_Nvol);
182 
183  vmax = 0.0;
184  for (int ex = 0; ex < m_Nex; ++ex) {
185  for (int site = is; site < ns; ++site) {
186  double fst = 0.0;
187  for (int in = 0; in < m_Nin; ++in) {
188  double fv = yp[myindex(in, site, ex)];
189  fst += fv * fv;
190  }
191  sum2 += fst;
192  fst = sqrt(fst);
193  sum += fst;
194  if (fst > vmax) vmax = fst;
195  }
196  }
197 
198  ThreadManager::reduce_sum_global(sum, ith, nth);
199  ThreadManager::reduce_sum_global(sum2, ith, nth);
200  ThreadManager::reduce_max_global(vmax, ith, nth);
201 
202  double vfac = double(m_Nvol) * double(m_Nex)
203  * double(CommonParameters::NPE());
204  Fave = sum / vfac;
205  Fdev = sqrt(sum2 / vfac - Fave * Fave);
206  Fmax = vmax;
207 
208 #pragma omp barrier
209 }
210 
211 
212 //====================================================================
213 void copy(Field& y, const Field& x)
214 {
215  int Nin = y.nin();
216  int Nvol = y.nvol();
217  int Nex = y.nex();
218  assert(x.check_size(Nin, Nvol, Nex));
219 
220  double *yp = y.ptr(0);
221  const double *RESTRICT xp = x.ptr(0);
222 
223  int ith, nth, is, ns;
224  set_threadtask(ith, nth, is, ns, Nvol);
225 
226  for (int ex = 0; ex < Nex; ++ex) {
227  for (int site = is; site < ns; ++site) {
228  int kv = Nin * (site + Nvol * ex);
229  for (int in = 0; in < Nin; ++in) {
230  yp[in + kv] = xp[in + kv];
231  }
232  }
233  }
234 }
235 
236 
237 //====================================================================
238 void copy(Field& y, const int exy, const Field& x, const int exx)
239 {
240  int Nin = y.nin();
241  int Nvol = y.nvol();
242 
243  assert(x.nin() == Nin);
244  assert(x.nvol() == Nvol);
245 
246  double *yp = y.ptr(0, 0, exy);
247  const double *RESTRICT xp = x.ptr(0, 0, exx);
248 
249  int ith, nth, is, ns;
250  set_threadtask(ith, nth, is, ns, Nvol);
251 
252  for (int site = is; site < ns; ++site) {
253  int kv = Nin * site;
254  for (int in = 0; in < Nin; ++in) {
255  yp[in + kv] = xp[in + kv];
256  }
257  }
258 }
259 
260 
261 //====================================================================
262 void scal(Field& x, const double a)
263 {
264  int Nin = x.nin();
265  int Nvol = x.nvol();
266  int Nex = x.nex();
267 
268  double *RESTRICT xp = x.ptr(0);
269 
270  int ith, nth, is, ns;
271  set_threadtask(ith, nth, is, ns, Nvol);
272 
273  for (int ex = 0; ex < Nex; ++ex) {
274  for (int site = is; site < ns; ++site) {
275  int kv = Nin * (site + Nvol * ex);
276  for (int in = 0; in < Nin; ++in) {
277  xp[in + kv] *= a;
278  }
279  }
280  }
281 }
282 
283 
284 //====================================================================
285 void scal(Field& x, const int exx, const double a)
286 {
287  int Nin = x.nin();
288  int Nvol = x.nvol();
289 
290  double *RESTRICT xp = x.ptr(0, 0, exx);
291 
292  int ith, nth, is, ns;
293  set_threadtask(ith, nth, is, ns, Nvol);
294 
295  for (int site = is; site < ns; ++site) {
296  int kv = Nin * site;
297  for (int in = 0; in < Nin; ++in) {
298  xp[in + kv] *= a;
299  }
300  }
301 }
302 
303 
304 //====================================================================
305 void scal(Field& x, const dcomplex a)
306 {
307  if (imag(a) == 0.0) {
308  return scal(x, real(a));
309  } else if (x.field_element_type() == Element_type::COMPLEX) {
310  int Nin = x.nin();
311  int Nvol = x.nvol();
312  int Nex = x.nex();
313 
314  int ith, nth, is, ns;
315  set_threadtask(ith, nth, is, ns, Nvol);
316 
317  double *RESTRICT xp = x.ptr(0);
318  double ar = real(a);
319  double ai = imag(a);
320  int Nin2 = Nin / 2;
321 
322  for (int ex = 0; ex < Nex; ++ex) {
323  for (int site = is; site < ns; ++site) {
324  int kv = Nin * (site + Nvol * ex);
325  for (int k = 0; k < Nin2; ++k) {
326  int kr = 2 * k;
327  int ki = 2 * k + 1;
328  double xr = xp[kr + kv];
329  double xi = xp[ki + kv];
330  xp[kr + kv] = ar * xr - ai * xi;
331  xp[ki + kv] = ar * xi + ai * xr;
332  }
333  }
334  }
335  } else {
336  vout.crucial("Error at %s: real vector and complex parameter.\n",
337  __func__);
338  exit(EXIT_FAILURE);
339  }
340 }
341 
342 
343 //====================================================================
344 void scal(Field& x, const int exx, const dcomplex a)
345 {
346  if (imag(a) == 0.0) {
347  return scal(x, exx, real(a));
348  } else if (x.field_element_type() == Element_type::COMPLEX) {
349  int Nin = x.nin();
350  int Nvol = x.nvol();
351  int Nex = x.nex();
352 
353  int ith, nth, is, ns;
354  set_threadtask(ith, nth, is, ns, Nvol);
355 
356  double *RESTRICT xp = x.ptr(0, 0, exx);
357  double ar = real(a);
358  double ai = imag(a);
359  int Nin2 = Nin / 2;
360 
361  for (int site = is; site < ns; ++site) {
362  int kv = Nin * site;
363  for (int k = 0; k < Nin2; ++k) {
364  int kr = 2 * k;
365  int ki = 2 * k + 1;
366  double xr = xp[kr + kv];
367  double xi = xp[ki + kv];
368  xp[kr + kv] = ar * xr - ai * xi;
369  xp[ki + kv] = ar * xi + ai * xr;
370  }
371  }
372  } else {
373  vout.crucial("Error at %s: real vector and complex parameter.\n",
374  __func__);
375  exit(EXIT_FAILURE);
376  }
377 }
378 
379 
380 //====================================================================
381 void axpy(Field& y, const double a, const Field& x)
382 {
383  int Nin = y.nin();
384  int Nvol = y.nvol();
385  int Nex = y.nex();
386  assert(x.check_size(Nin, Nvol, Nex));
387 
388  double *RESTRICT yp = y.ptr(0);
389  const double *RESTRICT xp = x.ptr(0);
390 
391  int ith, nth, is, ns;
392  set_threadtask(ith, nth, is, ns, Nvol);
393 
394  for (int ex = 0; ex < Nex; ++ex) {
395  for (int site = is; site < ns; ++site) {
396  int kv = Nin * (site + Nvol * ex);
397  for (int in = 0; in < Nin; ++in) {
398  yp[in + kv] += a * xp[in + kv];
399  }
400  }
401  }
402 }
403 
404 
405 //====================================================================
406 void axpy(Field& y, const int exy,
407  const double a, const Field& x, const int exx)
408 {
409  assert(x.nin() == y.nin());
410  assert(x.nvol() == y.nvol());
411 
412  int Nin = y.nin();
413  int Nvol = y.nvol();
414 
415  double *RESTRICT yp = y.ptr(0, 0, exy);
416  const double *RESTRICT xp = x.ptr(0, 0, exx);
417 
418  int ith, nth, is, ns;
419  set_threadtask(ith, nth, is, ns, Nvol);
420 
421  for (int site = is; site < ns; ++site) {
422  int kv = Nin * site;
423  for (int in = 0; in < Nin; ++in) {
424  yp[in + kv] += a * xp[in + kv];
425  }
426  }
427 }
428 
429 
430 //====================================================================
431 void axpy(Field& y, const dcomplex a, const Field& x)
432 {
433  if (imag(a) == 0.0) {
434  return axpy(y, real(a), x);
435  } else if ((y.field_element_type() == Element_type::COMPLEX) &&
436  (x.field_element_type() == Element_type::COMPLEX)) {
437  int Nin = y.nin();
438  int Nvol = y.nvol();
439  int Nex = y.nex();
440  assert(x.check_size(Nin, Nvol, Nex));
441 
442  double *RESTRICT yp = y.ptr(0);
443  const double *RESTRICT xp = x.ptr(0);
444 
445  int ith, nth, is, ns;
446  set_threadtask(ith, nth, is, ns, Nvol);
447 
448  double ar = real(a);
449  double ai = imag(a);
450  int Nin2 = Nin / 2;
451 
452  for (int ex = 0; ex < Nex; ++ex) {
453  for (int site = is; site < ns; ++site) {
454  int kv = Nin * (site + Nvol * ex);
455  for (int k = 0; k < Nin2; ++k) {
456  int kr = 2 * k;
457  int ki = 2 * k + 1;
458  yp[kr + kv] += ar * xp[kr + kv] - ai * xp[ki + kv];
459  yp[ki + kv] += ar * xp[ki + kv] + ai * xp[kr + kv];
460  }
461  }
462  }
463  } else {
464  vout.crucial("Error at %s: unsupported types.\n", __func__);
465  exit(EXIT_FAILURE);
466  }
467 }
468 
469 
470 //====================================================================
471 void axpy(Field& y, const int exy,
472  const dcomplex a, const Field& x, const int exx)
473 {
474  if (imag(a) == 0.0) {
475  return axpy(y, exy, real(a), x, exx);
476  } else if ((y.field_element_type() == Element_type::COMPLEX) &&
477  (x.field_element_type() == Element_type::COMPLEX)) {
478  int Nin = y.nin();
479  int Nvol = y.nvol();
480  assert(x.nin() == Nin);
481  assert(x.nvol() == Nvol);
482 
483  double *RESTRICT yp = y.ptr(0, 0, exy);
484  const double *RESTRICT xp = x.ptr(0, 0, exx);
485 
486  int ith, nth, is, ns;
487  set_threadtask(ith, nth, is, ns, Nvol);
488 
489  double ar = real(a);
490  double ai = imag(a);
491  int Nin2 = Nin / 2;
492 
493  for (int site = is; site < ns; ++site) {
494  int kv = Nin * site;
495  for (int k = 0; k < Nin2; ++k) {
496  int kr = 2 * k;
497  int ki = 2 * k + 1;
498  yp[kr + kv] += ar * xp[kr + kv] - ai * xp[ki + kv];
499  yp[ki + kv] += ar * xp[ki + kv] + ai * xp[kr + kv];
500  }
501  }
502  } else {
503  vout.crucial("Error at %s: unsupported types.\n", __func__);
504  exit(EXIT_FAILURE);
505  }
506 }
507 
508 
509 //====================================================================
510 void aypx(const double a, Field& y, const Field& x)
511 {
512  int Nin = y.nin();
513  int Nvol = y.nvol();
514  int Nex = y.nex();
515  assert(x.check_size(Nin, Nvol, Nex));
516 
517  double *RESTRICT yp = y.ptr(0);
518  const double *RESTRICT xp = x.ptr(0);
519 
520  int ith, nth, is, ns;
521  set_threadtask(ith, nth, is, ns, Nvol);
522 
523  for (int ex = 0; ex < Nex; ++ex) {
524  for (int site = is; site < ns; ++site) {
525  int kv = Nin * (site + Nvol * ex);
526  for (int in = 0; in < Nin; ++in) {
527  yp[in + kv] = a * yp[in + kv] + xp[in + kv];
528  }
529  }
530  }
531 }
532 
533 
534 //====================================================================
535 void aypx(const dcomplex a, Field& y, const Field& x)
536 {
537  if (imag(a) == 0.0) {
538  return aypx(real(a), y, x);
539  } else if ((y.field_element_type() == Element_type::COMPLEX) &&
540  (x.field_element_type() == Element_type::COMPLEX)) {
541  int Nin = y.nin();
542  int Nvol = y.nvol();
543  int Nex = y.nex();
544  assert(x.check_size(Nin, Nvol, Nex));
545 
546  double *RESTRICT yp = y.ptr(0);
547  const double *RESTRICT xp = x.ptr(0);
548 
549  int ith, nth, is, ns;
550  set_threadtask(ith, nth, is, ns, Nvol);
551 
552  double ar = real(a);
553  double ai = imag(a);
554  int Nin2 = Nin / 2;
555 
556  for (int ex = 0; ex < Nex; ++ex) {
557  for (int site = is; site < ns; ++site) {
558  int kv = Nin * (site + Nvol * ex);
559  for (int k = 0; k < Nin2; ++k) {
560  int kr = 2 * k;
561  int ki = 2 * k + 1;
562  double yr = yp[kr + kv];
563  double yi = yp[ki + kv];
564  yp[kr + kv] = ar * yr - ai * yi + xp[kr + kv];
565  yp[ki + kv] = ar * yi + ai * yr + xp[ki + kv];
566  }
567  }
568  }
569  } else {
570  vout.crucial("Error at %s: unsupported types.\n", __func__);
571  exit(EXIT_FAILURE);
572  }
573 }
574 
575 
576 //====================================================================
577 double dot(const Field& y, const Field& x)
578 {
579  int Nin = y.nin();
580  int Nvol = y.nvol();
581  int Nex = y.nex();
582  assert(x.check_size(Nin, Nvol, Nex));
583 
584  const double *RESTRICT yp = y.ptr(0);
585  const double *RESTRICT xp = x.ptr(0);
586 
587  int ith, nth, is, ns;
588  set_threadtask(ith, nth, is, ns, Nvol);
589 
590  double a = 0.0;
591 
592  for (int ex = 0; ex < Nex; ++ex) {
593  for (int site = is; site < ns; ++site) {
594  int kv = Nin * (site + Nvol * ex);
595  for (int in = 0; in < Nin; ++in) {
596  a += yp[in + kv] * xp[in + kv];
597  }
598  }
599  }
600 
601  ThreadManager::reduce_sum_global(a, ith, nth);
602 
603  return a;
604 }
605 
606 
607 //====================================================================
608 double dot(const Field& y, const int exy, const Field& x, const int exx)
609 {
610  assert(x.nin() == y.nin());
611  assert(x.nvol() == y.nvol());
612 
613  int Nin = y.nin();
614  int Nvol = y.nvol();
615 
616  const double *RESTRICT yp = y.ptr(0, 0, exy);
617  const double *RESTRICT xp = x.ptr(0, 0, exx);
618 
619  int ith, nth, is, ns;
620  set_threadtask(ith, nth, is, ns, Nvol);
621 
622  double a = 0.0;
623 
624  for (int site = is; site < ns; ++site) {
625  int kv = Nin * site;
626  for (int in = 0; in < Nin; ++in) {
627  a += yp[in + kv] * xp[in + kv];
628  }
629  }
630 
631  ThreadManager::reduce_sum_global(a, ith, nth);
632 
633  return a;
634 }
635 
636 
637 //====================================================================
638 void dot_and_norm2(double& yx, double& y2, double& x2,
639  const Field& y, const int exy,
640  const Field& x, const int exx)
641 {
642  int Nin = y.nin();
643  int Nvol = y.nvol();
644  assert(x.nin() == Nin);
645  assert(x.nvol() == Nvol);
646 
647  const double *RESTRICT yp = y.ptr(0, 0, exy);
648  const double *RESTRICT xp = x.ptr(0, 0, exx);
649 
650  int ith, nth, is, ns;
651  set_threadtask(ith, nth, is, ns, Nvol);
652 
653  double sum_yx = 0.0;
654  double sum_x2 = 0.0;
655  double sum_y2 = 0.0;
656 
657  for (int site = is; site < ns; ++site) {
658  int kv = Nin * site;
659  for (int in = 0; in < Nin; ++in) {
660  sum_yx += yp[in + kv] * xp[in + kv];
661  sum_x2 += xp[in + kv] * xp[in + kv];
662  sum_y2 += yp[in + kv] * yp[in + kv];
663  }
664  }
665 
666  double prd[3] = { sum_yx, sum_x2, sum_y2 };
667  ThreadManager::reduce_sum_global(prd, 3, ith, nth);
668  yx = prd[0];
669  y2 = prd[1];
670  x2 = prd[2];
671 }
672 
673 
674 //====================================================================
675 void dot_and_norm2(double& yx, double& y2, double& x2,
676  const Field& y, const Field& x)
677 {
678  int Nin = y.nin();
679  int Nvol = y.nvol();
680  int Nex = y.nex();
681  assert(x.check_size(Nin, Nvol, Nex));
682 
683  const double *RESTRICT yp = y.ptr(0);
684  const double *RESTRICT xp = x.ptr(0);
685 
686  int ith, nth, is, ns;
687  set_threadtask(ith, nth, is, ns, Nvol);
688 
689  double sum_yx = 0.0;
690  double sum_x2 = 0.0;
691  double sum_y2 = 0.0;
692 
693  for (int ex = 0; ex < Nex; ++ex) {
694  for (int site = is; site < ns; ++site) {
695  int kv = Nin * (site + Nvol * ex);
696  for (int in = 0; in < Nin; ++in) {
697  sum_yx += yp[in + kv] * xp[in + kv];
698  sum_x2 += xp[in + kv] * xp[in + kv];
699  sum_y2 += yp[in + kv] * yp[in + kv];
700  }
701  }
702  }
703 
704  double prd[3] = { sum_yx, sum_x2, sum_y2 };
705  ThreadManager::reduce_sum_global(prd, 3, ith, nth);
706  yx = prd[0];
707  y2 = prd[1];
708  x2 = prd[2];
709 }
710 
711 
712 //====================================================================
713 dcomplex dotc(const Field& y, const Field& x)
714 {
715  int Nin = y.nin();
716  int Nvol = y.nvol();
717  int Nex = y.nex();
718  assert(x.check_size(Nin, Nvol, Nex));
719 
720  const double *RESTRICT yp = y.ptr(0);
721  const double *RESTRICT xp = x.ptr(0);
722 
723  int ith, nth, is, ns;
724  set_threadtask(ith, nth, is, ns, Nvol);
725 
726  if ((y.field_element_type() == Element_type::COMPLEX) &&
727  (x.field_element_type() == Element_type::COMPLEX)) {
728  double prdr = 0.0;
729  double prdi = 0.0;
730  int Nin2 = Nin / 2;
731 
732  for (int ex = 0; ex < Nex; ++ex) {
733  for (int site = is; site < ns; ++site) {
734  int kv = Nin * (site + Nvol * ex);
735  for (int k = 0; k < Nin2; ++k) {
736  int kr = 2 * k;
737  int ki = 2 * k + 1;
738  prdr += yp[kr + kv] * xp[kr + kv] + yp[ki + kv] * xp[ki + kv];
739  prdi += yp[kr + kv] * xp[ki + kv] - yp[ki + kv] * xp[kr + kv];
740  }
741  }
742  }
743 
744  double prd[2] = { prdr, prdi };
745  ThreadManager::reduce_sum_global(prd, 2, ith, nth);
746 
747  return cmplx(prd[0], prd[1]);
748  } else if ((y.field_element_type() == Element_type::REAL) &&
749  (y.field_element_type() == Element_type::REAL)) {
750  return cmplx(dot(y, x), 0.0);
751  } else {
752  vout.crucial("Error at %s: unsupported arg types\n", __func__);
753  exit(EXIT_FAILURE);
754 
755  return cmplx(0.0, 0.0); // never reached.
756  }
757 }
758 
759 
760 //====================================================================
761 dcomplex dotc(const Field& y, const int exy, const Field& x, const int exx)
762 {
763  int Nin = y.nin();
764  int Nvol = y.nvol();
765  assert(x.nin() == Nin);
766  assert(x.nvol() == Nvol);
767 
768  const double *RESTRICT yp = y.ptr(0, 0, exy);
769  const double *RESTRICT xp = x.ptr(0, 0, exx);
770 
771  int ith, nth, is, ns;
772  set_threadtask(ith, nth, is, ns, Nvol);
773 
774  if ((y.field_element_type() == Element_type::COMPLEX) &&
775  (x.field_element_type() == Element_type::COMPLEX)) {
776  double prdr = 0.0;
777  double prdi = 0.0;
778  int Nin2 = Nin / 2;
779 
780  for (int site = is; site < ns; ++site) {
781  int kv = Nin * site;
782  for (int k = 0; k < Nin2; ++k) {
783  int kr = 2 * k;
784  int ki = 2 * k + 1;
785  prdr += yp[kr + kv] * xp[kr + kv] + yp[ki + kv] * xp[ki + kv];
786  prdi += yp[kr + kv] * xp[ki + kv] - yp[ki + kv] * xp[kr + kv];
787  }
788  }
789 
790  double prd[2] = { prdr, prdi };
791  ThreadManager::reduce_sum_global(prd, 2, ith, nth);
792 
793  return cmplx(prd[0], prd[1]);
794  } else if ((y.field_element_type() == Element_type::REAL) &&
795  (y.field_element_type() == Element_type::REAL)) {
796  return cmplx(dot(y, exy, x, exx), 0.0);
797  } else {
798  vout.crucial("Error at %s: unsupported arg types\n", __func__);
799  exit(EXIT_FAILURE);
800 
801  return cmplx(0.0, 0.0); // never reached.
802  }
803 }
804 
805 
806 //====================================================================
807 void dotc_and_norm2(dcomplex& yx, double& y2, double& x2,
808  const Field& y, const Field& x)
809 {
810  int Nin = y.nin();
811  int Nvol = y.nvol();
812  int Nex = y.nex();
813  assert(x.check_size(Nin, Nvol, Nex));
814 
815  const double *RESTRICT yp = y.ptr(0);
816  const double *RESTRICT xp = x.ptr(0);
817 
818  int ith, nth, is, ns;
819  set_threadtask(ith, nth, is, ns, Nvol);
820 
821  if ((y.field_element_type() == Element_type::COMPLEX) &&
822  (x.field_element_type() == Element_type::COMPLEX)) {
823  double prd_r = 0.0;
824  double prd_i = 0.0;
825  double prd_x2 = 0.0;
826  double prd_y2 = 0.0;
827  int Nin2 = Nin / 2;
828 
829  for (int ex = 0; ex < Nex; ++ex) {
830  for (int site = is; site < ns; ++site) {
831  int kv = Nin * (site + Nvol * ex);
832  for (int k = 0; k < Nin2; ++k) {
833  int kr = 2 * k;
834  int ki = 2 * k + 1;
835  prd_r += yp[kr + kv] * xp[kr + kv] + yp[ki + kv] * xp[ki + kv];
836  prd_i += yp[kr + kv] * xp[ki + kv] - yp[ki + kv] * xp[kr + kv];
837  prd_x2 += xp[kr + kv] * xp[kr + kv] + xp[ki + kv] * xp[ki + kv];
838  prd_y2 += yp[kr + kv] * yp[kr + kv] + yp[ki + kv] * yp[ki + kv];
839  }
840  }
841  }
842 
843  double prd[4] = { prd_r, prd_i, prd_x2, prd_y2 };
844  ThreadManager::reduce_sum_global(prd, 4, ith, nth);
845 
846  yx = cmplx(prd[0], prd[1]);
847  y2 = prd[2];
848  x2 = prd[3];
849  } else if ((y.field_element_type() == Element_type::REAL) &&
850  (y.field_element_type() == Element_type::REAL)) {
851  double yx_re = 0.0;
852  dot_and_norm2(yx_re, y2, x2, y, x);
853  yx = cmplx(yx_re, 0.0);
854  } else {
855  vout.crucial("Error at %s: unsupported arg types.\n", __func__);
856  exit(EXIT_FAILURE);
857  }
858 }
859 
860 
861 //====================================================================
862 void dotc_and_norm2(dcomplex& yx, double& y2, double& x2,
863  const Field& y, const int exy,
864  const Field& x, const int exx)
865 {
866  int Nin = y.nin();
867  int Nvol = y.nvol();
868  assert(x.nin() == Nin);
869  assert(x.nvol() == Nvol);
870 
871  const double *RESTRICT yp = y.ptr(0, 0, exy);
872  const double *RESTRICT xp = x.ptr(0, 0, exx);
873 
874  int ith, nth, is, ns;
875  set_threadtask(ith, nth, is, ns, Nvol);
876 
877  if ((y.field_element_type() == Element_type::COMPLEX) &&
878  (x.field_element_type() == Element_type::COMPLEX)) {
879  double prd_r = 0.0;
880  double prd_i = 0.0;
881  double prd_x2 = 0.0;
882  double prd_y2 = 0.0;
883  int Nin2 = Nin / 2;
884 
885  for (int site = is; site < ns; ++site) {
886  int kv = Nin * site;
887  for (int k = 0; k < Nin2; ++k) {
888  int kr = 2 * k;
889  int ki = 2 * k + 1;
890  prd_r += yp[kr + kv] * xp[kr + kv] + yp[ki + kv] * xp[ki + kv];
891  prd_i += yp[kr + kv] * xp[ki + kv] - yp[ki + kv] * xp[kr + kv];
892  prd_x2 += xp[kr + kv] * xp[kr + kv] + xp[ki + kv] * xp[ki + kv];
893  prd_y2 += yp[kr + kv] * yp[kr + kv] + yp[ki + kv] * yp[ki + kv];
894  }
895  }
896 
897  double prd[4] = { prd_r, prd_i, prd_x2, prd_y2 };
898  ThreadManager::reduce_sum_global(prd, 4, ith, nth);
899 
900  yx = cmplx(prd[0], prd[1]);
901  y2 = prd[2];
902  x2 = prd[3];
903  } else if ((y.field_element_type() == Element_type::REAL) &&
904  (y.field_element_type() == Element_type::REAL)) {
905  double yx_re = 0.0;
906  dot_and_norm2(yx_re, y2, x2, y, exy, x, exx);
907  yx = cmplx(yx_re, 0.0);
908  } else {
909  vout.crucial("Error at %s: unsupported arg types.\n", __func__);
910  exit(EXIT_FAILURE);
911  }
912 }
913 
914 
915 //============================================================END=====
Field::myindex
size_t myindex(const int jin, const int site, const int jex) const
Definition: field.h:64
bridgeIO.h
dot_and_norm2
void dot_and_norm2(double &yx, double &y2, double &x2, const Field &y, const int exy, const Field &x, const int exx)
Definition: field.cpp:638
Field::m_Nex
int m_Nex
external d.o.f.
Definition: field.h:57
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
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
aypx
void aypx(const double a, Field &y, const Field &x)
aypx(y, a, x): y := a * y + x
Definition: field.cpp:510
Field::stat
void stat(double &Fave, double &Fmax, double &Fdev) const
determines the statistics of the field. average, maximum value, and deviation is determined over glob...
Definition: field.cpp:170
Field::xI
void xI()
Definition: field.cpp:139
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:381
dot
double dot(const Field &y, const Field &x)
Definition: field.cpp:577
dotc_and_norm2
void dotc_and_norm2(dcomplex &yx, double &y2, double &x2, const Field &y, const Field &x)
Definition: field.cpp:807
Field::nin
int nin() const
Definition: field.h:126
Field::m_element_type
element_type m_element_type
field complex type
Definition: field.h:58
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:213
Field::norm2
double norm2() const
Definition: field.cpp:114
ThreadManager::reduce_max_global
static void reduce_max_global(double *value, const int num, const int i_thread, const int Nthread)
global reduction with max for an array: double values are assumed thread local.
Definition: threadManager.cpp:361
Field::check
void check()
Definition: field.cpp:29
Field::class_name
static const std::string class_name
Definition: field.h:52
ThreadManager::reduce_sum_global
static void reduce_sum_global(dcomplex &value, const int i_thread, const int Nthread)
global reduction with summation: dcomplex values are assumed thread local.
Definition: threadManager.cpp:288
threadManager.h
field.h
Element_type::REAL
@ REAL
Definition: bridge_defs.h:43
Field::nvol
int nvol() const
Definition: field.h:127
dotc
dcomplex dotc(const Field &y, const Field &x)
Definition: field.cpp:713
CommonParameters::NPE
static int NPE()
Definition: commonParameters.h:101
Field::field_element_type
element_type field_element_type() const
Definition: field.h:129
Field::setc
void setc(double a)
Definition: field.cpp:56
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
Element_type::COMPLEX
@ COMPLEX
Definition: bridge_defs.h:43
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:262
field_thread-inc.h
Field::m_Nin
int m_Nin
internal d.o.f.
Definition: field.h:55
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:242
Field
Container of Field-type object.
Definition: field.h:46
bridge_defs.h
Field::m_Nvol
int m_Nvol
lattice volume
Definition: field.h:56
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:572
CommonParameters::epsilon_criterion
static double epsilon_criterion()
Definition: commonParameters.h:119