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