Bridge++  Ver. 1.2.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
field.cpp
Go to the documentation of this file.
1 
14 #include "field.h"
15 #include <cstring>
16 
17 #include "threadManager_OpenMP.h"
18 
19 #include "bridgeIO.h"
20 using Bridge::vout;
21 
22 using std::string;
23 
24 //====================================================================
25 namespace {
26  inline
27  void set_threadtask(int& ith, int& nth, int& is, int& ns,
28  const int size)
29  {
30  nth = ThreadManager_OpenMP::get_num_threads();
31  ith = ThreadManager_OpenMP::get_thread_id();
32 
33  is = size * ith / nth;
34  ns = size * (ith + 1) / nth;
35  }
36 }
37 
38 //====================================================================
39 void Field::check()
40 {
41  // vout.general("Field was constructed.\n");
42 }
43 
44 
45 //====================================================================
46 double dot(const Field& y, const Field& x)
47 {
48  double *yp = const_cast<Field *>(&y)->ptr(0);
49  double *xp = const_cast<Field *>(&x)->ptr(0);
50 
51  // int size = x.ntot();
52  int ith, nth, is, ns;
53 
54  set_threadtask(ith, nth, is, ns, x.ntot());
55 
56  double a = 0.0;
57 
58  // for (int k = 0; k < size; ++k) {
59  for (int k = is; k < ns; ++k) {
60  a += yp[k] * xp[k];
61  }
62 
63  // double b = Communicator::reduce_sum(a);
64  // return b;
65  ThreadManager_OpenMP::reduce_sum_global(a, ith, nth);
66  return a;
67 }
68 
69 
70 //====================================================================
71 double dot(const Field& y, const int exy, const Field& x, const int exx)
72 {
73  assert(x.nin() == y.nin());
74  assert(x.nvol() == y.nvol());
75 
76  double *yp = const_cast<Field *>(&y)->ptr(0, 0, exy);
77  double *xp = const_cast<Field *>(&x)->ptr(0, 0, exx);
78 
79  // int size = x.nin() * x.nvol();
80  int ith, nth, is, ns;
81  set_threadtask(ith, nth, is, ns, x.ntot());
82 
83  double a = 0.0;
84 
85  // for (int k = 0; k < size; ++k) {
86  for (int k = is; k < ns; ++k) {
87  a += yp[k] * xp[k];
88  }
89 
90  // double b = Communicator::reduce_sum(a);
91  // return b;
92  ThreadManager_OpenMP::reduce_sum_global(a, ith, nth);
93  return a;
94 }
95 
96 
97 //====================================================================
98 dcomplex dotc(const Field& y, const Field& x)
99 {
100  double *yp = const_cast<Field *>(&y)->ptr(0);
101  double *xp = const_cast<Field *>(&x)->ptr(0);
102 
103  assert(x.ntot() == y.ntot());
104 
105  // int ntot = x.ntot();
106  int ith, nth, is, ns;
107  set_threadtask(ith, nth, is, ns, x.ntot());
108 
109  if ((y.field_element_type() == Field::COMPLEX) &&
110  (x.field_element_type() == Field::COMPLEX)) {
111  double prdr = 0.0;
112  double prdi = 0.0;
113 
114  // for (int k = 0; k < ntot; k += 2) {
115  for (int k = is; k < ns; k += 2) {
116  prdr += yp[k] * xp[k] + yp[k + 1] * xp[k + 1];
117  prdi += yp[k] * xp[k + 1] - yp[k + 1] * xp[k];
118  }
119 
120  // double prd[2];
121  // prd[0] = prdr;
122  // prd[1] = prdi;
123 
124  // double prd_g[2];
125  // Communicator::reduce_sum(2, prd_g, prd);
126 
127  // return cmplx(prd_g[0], prd_g[1]);
128 
129  ThreadManager_OpenMP::reduce_sum_global(prdr, ith, nth);
130  ThreadManager_OpenMP::reduce_sum_global(prdi, ith, nth);
131 
132  return cmplx(prdr, prdi);
133  } else if ((y.field_element_type() == Field::REAL) &&
134  (y.field_element_type() == Field::REAL)) {
135  return cmplx(dot(y, x), 0.0);
136  } else {
137  vout.crucial("%s: unsupported arg types.\n", __func__);
138  abort();
139  return cmplx(0.0, 0.0); // never reached.
140  }
141 }
142 
143 
144 //====================================================================
145 dcomplex dotc(const Field& y, const int exy, const Field& x, const int exx)
146 {
147  double *yp = const_cast<Field *>(&y)->ptr(0, 0, exy);
148  double *xp = const_cast<Field *>(&x)->ptr(0, 0, exx);
149 
150  assert(x.nin() == y.nin());
151  assert(x.nvol() == y.nvol());
152 
153  // int size = x.nin() * x.nvol();
154  int ith, nth, is, ns;
155  set_threadtask(ith, nth, is, ns, x.nin() * x.nvol());
156 
157  if ((y.field_element_type() == Field::COMPLEX) &&
158  (x.field_element_type() == Field::COMPLEX)) {
159  double prdr = 0.0;
160  double prdi = 0.0;
161 
162  // for (int k = 0; k < size; k += 2) {
163  for (int k = is; k < ns; k += 2) {
164  prdr += yp[k] * xp[k] + yp[k + 1] * xp[k + 1];
165  prdi += yp[k] * xp[k + 1] - yp[k + 1] * xp[k];
166  }
167 
168  // double prd[2];
169  // prd[0] = prdr;
170  // prd[1] = prdi;
171 
172  // double prd_g[2];
173  // Communicator::reduce_sum(2, prd_g, prd);
174 
175  // return cmplx(prd_g[0], prd_g[1]);
176 
177  ThreadManager_OpenMP::reduce_sum_global(prdr, ith, nth);
178  ThreadManager_OpenMP::reduce_sum_global(prdi, ith, nth);
179 
180  return cmplx(prdr, prdi);
181  } else if ((y.field_element_type() == Field::REAL) &&
182  (y.field_element_type() == Field::REAL)) {
183  return cmplx(dot(y, exy, x, exx), 0.0);
184  } else {
185  vout.crucial("%s: unsupported arg types.\n", __func__);
186  abort();
187  return cmplx(0.0, 0.0); // never reached.
188  }
189 }
190 
191 
192 //====================================================================
193 void axpy(Field& y, const double a, const Field& x)
194 {
195  double *yp = const_cast<Field *>(&y)->ptr(0);
196  double *xp = const_cast<Field *>(&x)->ptr(0);
197 
198  // int size = x.ntot();
199  int ith, nth, is, ns;
200 
201  set_threadtask(ith, nth, is, ns, x.ntot());
202 
203  // for (int k = 0; k < size; ++k) {
204  for (int k = is; k < ns; ++k) {
205  yp[k] += a * xp[k];
206  }
207 }
208 
209 
210 //====================================================================
211 void axpy(Field& y, const int exy, const double a, const Field& x, const int exx)
212 {
213  assert(x.nin() == y.nin());
214  assert(x.nvol() == y.nvol());
215 
216  double *yp = const_cast<Field *>(&y)->ptr(0, 0, exy);
217  double *xp = const_cast<Field *>(&x)->ptr(0, 0, exx);
218 
219  // int size = x.nin() * x.nvol();
220  int ith, nth, is, ns;
221  set_threadtask(ith, nth, is, ns, x.nin() * x.nvol());
222 
223  // for (int k = 0; k < size; ++k) {
224  for (int k = is; k < ns; ++k) {
225  yp[k] += a * xp[k];
226  }
227 }
228 
229 
230 //====================================================================
231 void axpy(Field& y, const dcomplex a, const Field& x)
232 {
233  if (imag(a) == 0.0) {
234  return axpy(y, real(a), x);
235  } else if ((y.field_element_type() == Field::REAL) &&
236  (x.field_element_type() == Field::REAL)) {
237  vout.crucial("%s: real vector and complex parameter.\n", __func__);
238  abort();
239 
240  return axpy(y, real(a), x); // ignore imaginary part of a
241  } else if ((y.field_element_type() == Field::COMPLEX) &&
242  (x.field_element_type() == Field::COMPLEX)) {
243  double *yp = const_cast<Field *>(&y)->ptr(0);
244  double *xp = const_cast<Field *>(&x)->ptr(0);
245 
246  assert(x.ntot() == y.ntot());
247 
248  // int ntot = x.ntot();
249  int ith, nth, is, ns;
250  set_threadtask(ith, nth, is, ns, x.ntot());
251 
252  double ar = real(a);
253  double ai = imag(a);
254 
255  // for (int k = 0; k < ntot; k += 2) {
256  for (int k = is; k < ns; k += 2) {
257  yp[k] += ar * xp[k] - ai * xp[k + 1];
258  yp[k + 1] += ar * xp[k + 1] + ai * xp[k];
259  }
260 
261  return;
262  } else {
263  vout.crucial("%s: unsupported types.\n", __func__);
264  abort();
265  }
266 }
267 
268 
269 //====================================================================
270 void axpy(Field& y, const int exy, const dcomplex a, const Field& x, const int exx)
271 {
272  if (imag(a) == 0.0) {
273  return axpy(y, exy, real(a), x, exx);
274  } else if ((y.field_element_type() == Field::REAL) &&
275  (x.field_element_type() == Field::REAL)) {
276  vout.crucial("%s: real vector and complex parameter.\n", __func__);
277  abort();
278 
279  return axpy(y, real(a), x); // ignore imaginary part of a
280  } else if ((y.field_element_type() == Field::COMPLEX) &&
281  (x.field_element_type() == Field::COMPLEX)) {
282  double *yp = const_cast<Field *>(&y)->ptr(0, 0, exy);
283  double *xp = const_cast<Field *>(&x)->ptr(0, 0, exx);
284 
285  assert(x.nin() == y.nin());
286  assert(x.nvol() == y.nvol());
287 
288  // int size = x.nin() * x.nvol();
289  int ith, nth, is, ns;
290  set_threadtask(ith, nth, is, ns, x.nin() * x.nvol());
291 
292  double ar = real(a);
293  double ai = imag(a);
294 
295  // for (int k = 0; k < size; k += 2) {
296  for (int k = is; k < ns; k += 2) {
297  yp[k] += ar * xp[k] - ai * xp[k + 1];
298  yp[k + 1] += ar * xp[k + 1] + ai * xp[k];
299  }
300 
301  return;
302  } else {
303  vout.crucial("%s: unsupported types.\n", __func__);
304  abort();
305  }
306 }
307 
308 
309 //====================================================================
310 void scal(Field& x, const double a)
311 {
312  // x.field *= a;
313 
314  double *xp = const_cast<Field *>(&x)->ptr(0, 0, 0);
315  int ith, nth, is, ns;
316 
317  set_threadtask(ith, nth, is, ns, x.ntot());
318 
319  for (int k = is; k < ns; ++k) {
320  xp[k] *= a;
321  }
322 }
323 
324 
325 //====================================================================
326 void scal(Field& x, const int exx, const double a)
327 {
328  double *xp = const_cast<Field *>(&x)->ptr(0, 0, exx);
329  // int size = x.nin() * x.nvol();
330  int ith, nth, is, ns;
331 
332  set_threadtask(ith, nth, is, ns, x.nin() * x.nvol());
333 
334  // for (int k = 0; k < size; ++k) {
335  for (int k = is; k < ns; ++k) {
336  // x.field[k] *= a;
337  xp[k] *= a;
338  }
339 }
340 
341 
342 //====================================================================
343 void scal(Field& x, const dcomplex a)
344 {
345  if (x.field_element_type() == Field::REAL) {
346  vout.crucial("%s: real vector and complex parameter.\n", __func__);
347  abort();
348 
349  // x.field *= real(a); // ignore imaginary part of a
350  scal(x, a);
351  } else if (x.field_element_type() == Field::COMPLEX) {
352  double *xp = const_cast<Field *>(&x)->ptr(0);
353  //int ntot = x.ntot();
354  int ith, nth, is, ns;
355  set_threadtask(ith, nth, is, ns, x.ntot());
356 
357  double ar = real(a);
358  double ai = imag(a);
359 
360  // for (int k = 0; k < ntot; k += 2) {
361  for (int k = is; k < ns; k += 2) {
362  double xr = xp[k];
363  double xi = xp[k + 1];
364 
365  xp[k] = ar * xr - ai * xi;
366  xp[k + 1] = ar * xi + ai * xr;
367  }
368  } else {
369  vout.crucial("%s: unsupported field type.\n", __func__);
370  abort();
371  }
372 }
373 
374 
375 //====================================================================
376 void scal(Field& x, const int exx, const dcomplex a)
377 {
378  if (x.field_element_type() == Field::REAL) {
379  vout.crucial("%s: real vector and complex parameter.\n", __func__);
380  abort();
381 
382  // x.field *= real(a); // ignore imaginary part of a
383  scal(x, exx, a);
384  } else if (x.field_element_type() == Field::COMPLEX) {
385  double *xp = const_cast<Field *>(&x)->ptr(0, 0, exx);
386  // int size = x.nin() * x.nvol();
387  int ith, nth, is, ns;
388  set_threadtask(ith, nth, is, ns, x.nin() * x.nvol());
389 
390  double ar = real(a);
391  double ai = imag(a);
392 
393  // for (int k = 0; k < size; k += 2) {
394  for (int k = is; k < ns; k += 2) {
395  double xr = xp[k];
396  double xi = xp[k + 1];
397 
398  xp[k] = ar * xr - ai * xi;
399  xp[k + 1] = ar * xi + ai * xr;
400  }
401  } else {
402  vout.crucial("%s: unsupported field type.\n", __func__);
403  abort();
404  }
405 }
406 
407 
408 //====================================================================
409 void copy(Field& y, const Field& x)
410 {
411  // y.field = x.field;
412 
413  assert(x.nin() == y.nin());
414  assert(x.nvol() == y.nvol());
415  assert(x.nex() == y.nex());
416 
417  double *yp = const_cast<Field *>(&y)->ptr(0, 0, 0);
418  double *xp = const_cast<Field *>(&x)->ptr(0, 0, 0);
419 
420  int ith, nth, is, ns;
421  set_threadtask(ith, nth, is, ns, x.ntot());
422 
423  for (int k = is; k < ns; ++k) {
424  yp[k] = xp[k];
425  }
426 }
427 
428 
429 //====================================================================
430 void copy(Field& y, const int exy, const Field& x, const int exx)
431 {
432  assert(x.nin() == y.nin());
433  assert(x.nvol() == y.nvol());
434 
435  double *yp = const_cast<Field *>(&y)->ptr(0, 0, exy);
436  double *xp = const_cast<Field *>(&x)->ptr(0, 0, exx);
437 
438  // int size = x.nin() * x.nvol();
439  int ith, nth, is, ns;
440  set_threadtask(ith, nth, is, ns, x.nin() * x.nvol());
441 
442  for (int k = is; k < ns; ++k) {
443  yp[k] = xp[k];
444  }
445 
449  // memcpy(yp, xp, sizeof(double) * size);
450 }
451 
452 
453 //====================================================================
454 void Field::set(double a)
455 {
456  int ith, nth, is, ns;
457 
458  set_threadtask(ith, nth, is, ns, ntot());
459 
460  double *yp = const_cast<Field *>(this)->ptr(0);
461 
462  for (int k = is; k < ns; ++k) {
463  yp[k] = a;
464  }
465 }
466 
467 
468 //====================================================================
469 double Field::norm2() const
470 {
471  int ith, nth, is, ns;
472 
473  set_threadtask(ith, nth, is, ns, ntot());
474 
475  double *yp = const_cast<Field *>(this)->ptr(0);
476 
477  double a = 0.0;
478 
479  for (int k = is; k < ns; ++k) {
480  a += yp[k] * yp[k];
481  }
482 
483  ThreadManager_OpenMP::reduce_sum_global(a, ith, nth);
484  return a;
485 }
486 
487 
488 //====================================================================
489 void aypx(const double a, Field& y, const Field& x)
490 {
491  int ith, nth, is, ns;
492 
493  set_threadtask(ith, nth, is, ns, x.ntot());
494 
495  double *yp = const_cast<Field *>(&y)->ptr(0);
496  double *xp = const_cast<Field *>(&x)->ptr(0);
497 
498  for (int k = is; k < ns; ++k) {
499  yp[k] = a * yp[k] + xp[k];
500  }
501 }
502 
503 
504 //====================================================================
505 void aypx(const dcomplex a, Field& y, const Field& x)
506 {
507  if (imag(a) == 0.0) {
508  return aypx(real(a), y, x);
509  } else if ((y.field_element_type() == Field::REAL) &&
510  (x.field_element_type() == Field::REAL)) {
511  vout.crucial("%s: real vector and complex parameter.\n", __func__);
512  abort();
513 
514  return aypx(real(a), y, x); // ignore imaginary part of a
515  } else if ((y.field_element_type() == Field::COMPLEX) &&
516  (x.field_element_type() == Field::COMPLEX)) {
517  double *yp = const_cast<Field *>(&y)->ptr(0);
518  double *xp = const_cast<Field *>(&x)->ptr(0);
519 
520  assert(x.ntot() == y.ntot());
521 
522  int ith, nth, is, ns;
523  set_threadtask(ith, nth, is, ns, x.ntot());
524 
525  double ar = real(a);
526  double ai = imag(a);
527 
528  for (int k = is; k < ns; k += 2) {
529  double ypr = yp[k];
530  double ypi = yp[k + 1];
531  yp[k] = ar * ypr - ai * ypi + xp[k];
532  yp[k + 1] = ar * ypi + ai * ypr + xp[k + 1];
533  }
534 
535  return;
536  } else {
537  vout.crucial("%s: unsupported types.\n", __func__);
538  abort();
539  }
540 }
541 
542 
543 //====================================================================
544 void Field::stat(double& Fave, double& Fmax, double& Fdev) const
545 {
546  assert(ThreadManager_OpenMP::get_num_threads() == 1);
547 
548  double sum = 0.0;
549  double sum2 = 0.0;
550 
551  Fmax = 0.0;
552  for (int ex = 0, nnex = nex(); ex < nnex; ++ex) {
553  for (int site = 0, nnvol = nvol(); site < nnvol; ++site) {
554  double fst = 0.0;
555  for (int in = 0, nnin = nin(); in < nnin; ++in) {
556  double fv = field[myindex(in, site, ex)];
557  fst += fv * fv;
558  }
559  sum2 += fst;
560  fst = sqrt(fst);
561  sum += fst;
562  if (fst > Fmax) Fmax = fst;
563  }
564  }
565 
566  sum = Communicator::reduce_sum(sum);
567  sum2 = Communicator::reduce_sum(sum2);
568  Fmax = Communicator::reduce_max(Fmax);
569 
570  double perdeg = 1.0 / ((double)CommonParameters::Lvol() * nex());
571  Fave = sum * perdeg;
572 
573  double fval = sum2 * perdeg;
574  fval -= Fave * Fave;
575  Fdev = sqrt(fval);
576 }
577 
578 
579 //====================================================================
580 void report_field_stat(const Bridge::VerboseLevel vl,
581  const std::string& msg, const Field& f)
582 {
583  assert(ThreadManager_OpenMP::get_num_threads() == 1);
584 
585  double favg, fmax, fdev;
586  f.stat(favg, fmax, fdev);
587 
588  vout.general(vl,
589  " %s: avg = %12.6f max = %12.6f dev = %12.6f\n",
590  msg.c_str(),
591  favg, fmax, fdev);
592 }
593 
594 
595 //====================================================================
596 //============================================================END=====
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:310
report_field_stat
void report_field_stat(const Bridge::VerboseLevel vl, const std::string &msg, const Field &f)
Definition: field.cpp:580
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:207
threadManager_OpenMP.h
ThreadManager_OpenMP::get_num_threads
static int get_num_threads()
returns available number of threads.
Definition: threadManager_OpenMP.cpp:95
Communicator::reduce_max
static int reduce_max(int count, double *recv_buf, double *send_buf, int pattern=0)
find a global maximum of an array of double over the communicator. pattern specifies the dimensions t...
Definition: communicator.cpp:244
Field::norm2
double norm2() const
Definition: field.cpp:469
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:128
dot
double dot(const Field &y, const Field &x)
Definition: field.cpp:46
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:38
Field::ptr
double * ptr(const int jin, const int site, const int jex)
Definition: field.h:118
Field
Container of Field-type object.
Definition: field.h:37
Field::nvol
int nvol() const
Definition: field.h:101
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:409
CommonParameters::Lvol
static int Lvol()
Definition: commonParameters.h:86
ThreadManager_OpenMP::reduce_sum_global
static void reduce_sum_global(double &value, const int ith, const int nth)
global reduction with summation: value is assumed thread local.
Definition: threadManager_OpenMP.cpp:155
ThreadManager_OpenMP::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager_OpenMP.cpp:107
Field::nin
int nin() const
Definition: field.h:100
dotc
dcomplex dotc(const Field &y, const Field &x)
Definition: field.cpp:98
aypx
void aypx(const double a, Field &y, const Field &x)
aypx(y, a, x): y := a * y + x
Definition: field.cpp:489
Field::check
void check()
Definition: field.cpp:39
Field::nex
int nex() const
Definition: field.h:102
Field::field
std::valarray< double > field
Definition: field.h:49
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:193
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:26
ParameterCheck::vl
Bridge::VerboseLevel vl
Definition: checker.cpp:18
Field::ntot
int ntot() const
Definition: field.h:105
Bridge::VerboseLevel
VerboseLevel
Definition: bridgeIO.h:25
Field::field_element_type
element_type field_element_type() const
Definition: field.h:103
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:544
Communicator::reduce_sum
static int reduce_sum(int count, double *recv_buf, double *send_buf, int pattern=0)
make a global sum of an array of double over the communicator. pattern specifies the dimensions to be...
Definition: communicator.cpp:224
Field::myindex
int myindex(const int jin, const int site, const int jex) const
Definition: field.h:52