Bridge++  Ver. 2.0.2
MultiGrid_Clover-tmpl.h
Go to the documentation of this file.
1 
10 #include "lib/Tools/randomNumberManager.h"
11 #include "lib/Tools/randomNumbers_Mseries.h"
12 #include "lib/Tools/timer.h"
13 
14 //====================================================================
15 namespace {
16  // memory saving implementation: no work array of fields are needed
17  template<typename INDEX, typename AFIELD>
18  void gramschmidt_in_block_impl(std::vector<AFIELD>& v,
19  AFopr<AFIELD> *afopr,
20  const INDEX& block_index,
21  AFIELD& tmp1, AFIELD& tmp2,
22  typename AFIELD::real_t *block_norm,
23  typename AFIELD::complex_t *block_prod)
24  {
25  typedef typename AFIELD::real_t real_t;
26  typedef typename AFIELD::complex_t complex_t;
27 
28  vout.detailed("gramschimdt_in_block, called (general version with smaller memory): size of(real_t)=%d\n", sizeof(real_t));
29 
30 #ifdef DEBUG_GS
31  {
32  double norm2 = v[0].norm2();
33  vout.general("initial: n2[0]=%e\n", norm2);
34  }
35 #endif
36 
37  int nvector = v.size();
38 
39  int coarse_nvol = block_index.coarse_nvol();
40  int ith, nth, is, ns;
41  set_threadtask(ith, nth, is, ns, coarse_nvol);
42 
43 
44  // make <v[j]|v[i]> = 0 for j < i
45  for (int i = 0; i < nvector; i++) {
46  // vout.general("i = %d norm2 = %f\n", i, v[i].norm2());
47 
48 #pragma omp barrier
49 
50  // linear algebras on the fine lattice
51  copy(tmp1, v[i]);
52 #pragma omp barrier
53 
54  afopr->mult_gm5(tmp2, v[i]);
55 
56  axpy(tmp1, real_t(1.0), tmp2); // tmp1 = (1 + gm5) v[i]
57  axpy(v[i], real_t(-1.0), tmp2); // v[i] := (1 - gm5) v[i]
58 #pragma omp barrier // work on the fine lattice, done
59 
60  // linear algebras with block
61  // ch=-1
62  for (int j = 0; j < i; j++) {
63  // P|v[i]> := P|v[i]> - |v[j]> <v[j]|P|v[i]>
64  block_dotc(block_prod, v[j], v[i], block_index);
65  block_axpy(v[i], block_prod, v[j], real_t(-1.0), block_index);
66  }
67 #pragma omp barrier // work on the blocks, done
68 
69  // normalize
70  // linear algebras on the fine lattice
71  afopr->mult_gm5(tmp2, v[i]);
72  axpy(v[i], real_t(-1.0), tmp2); // apply P(-)
73 #pragma omp barrier // work on the fine lattice, done
74 
75  // linear algebras with block
76  block_norm2(block_norm, v[i], block_index);
77 #pragma omp barrier // work on the blocks, done
78 
79  for (int k = is; k < ns; k++) {
80  double n2 = block_norm[k];
81  block_norm[k] = 1.0 / sqrt(n2);
82  }
83 #pragma omp barrier
84 
85  block_scal(v[i], block_norm, block_index);
86 
87 
88  // ch=+1
89  for (int j = 0; j < i; ++j) {
90  // |Pv[i]> := |Pv[i]> - |v[j]> <v[j]|Pv[i]>
91  block_dotc(block_prod, v[j], tmp1, block_index);
92  block_axpy(tmp1, block_prod, v[j], real_t(-1.0), block_index);
93  }
94 #pragma omp barrier // work on blocks, done
95 
96 
97  // normalize
98 
99  // linear algebras on the fine lattice
100  afopr->mult_gm5(tmp2, tmp1);
101  axpy(tmp1, real_t(1.0), tmp2); // apply P
102 #pragma omp barrier // work on the fine lattice, done
103 
104  // linear algebra with blocks
105  block_norm2(block_norm, tmp1, block_index);
106 #pragma omp barrier // work on the blocks, done
107 
108  for (int k = is; k < ns; ++k) {
109  double n2 = block_norm[k];
110  block_norm[k] = 1.0 / sqrt(n2);
111  }
112 #pragma omp barrier
113 
114  // merge both chiralities
115  block_axpy(v[i], block_norm, tmp1, real_t(1.0), block_index);
116  } // i
117 
118 #pragma omp barrier
119 
120 #ifdef DEBUG
121  {
122  double norm2 = v[0].norm2();
123  vout.general(" after Gramschmidt n2[0]=%e\n", norm2);
124  }
125 #endif
126  }
127 
128 
129  //====================================================================
130  template<typename INDEX, class AFIELD1, class AFIELD2>
131  void make_coarse_vector_impl(AFIELD1& coarse_vector, const AFIELD2& fine_vector,
132  const std::vector<AFIELD2>& testvectors,
133  AFopr<AFIELD2> *afopr, const INDEX& block_index,
134  AFIELD2& tmp1, AFIELD2& tmp2, typename AFIELD1::complex_t *inner_prod)
135  {
136 #ifdef DEBUG
137  vout.detailed("%s is called\n", __func__);
138 #endif
139  typedef typename AFIELD1::complex_t complex_t;
140  typedef typename AFIELD1::real_t real_t;
141 
142  const int num_vectors = testvectors.size();
143  assert(coarse_vector.nin() == 4 * num_vectors); // 4: 2 for complex, 2 for chirality
144 
145  int ith, nth, is, ns;
146  set_threadtask(ith, nth, is, ns, block_index.coarse_nvol());
147  int coarse_Nx = block_index.coarse_lattice_size(0);
148  int coarse_Ny = block_index.coarse_lattice_size(1);
149  int coarse_Nz = block_index.coarse_lattice_size(2);
150  int coarse_Nt = block_index.coarse_lattice_size(3);
151  AIndex_coarse_lex<real_t, AFIELD1::IMPL> index(coarse_Nx, coarse_Ny, coarse_Nz, coarse_Nt, num_vectors, 2);
152 
153 
154 
155  for (int i = 0; i < num_vectors; i++) {
156 #pragma omp barrier
157 
158  //linear algebra on fine lattice
159  copy(tmp2, testvectors[i]);
160  afopr->mult_gm5(tmp1, tmp2);
161  axpy(tmp2, -1.0, tmp1); // tmp2=(1-gm5)v[i]
162  axpy(tmp1, +1.0, testvectors[i]); // tmp1=(1+gm5)v[i]
163 
164 #pragma omp barrier
165 
166  // linear algerbar with blocks
167 
168  // ch=-1
169  block_dotc(inner_prod, tmp2, fine_vector, block_index);
170  for (int block_id = is; block_id < ns; block_id++) {
171  int index_re = index.idx_SPr(i, 0, block_id, 0);
172  int index_im = index.idx_SPi(i, 0, block_id, 0);
173  coarse_vector.set(index_re, 0.5 * real(inner_prod[block_id]));
174  coarse_vector.set(index_im, 0.5 * imag(inner_prod[block_id]));
175  }
176 
177  // ch=+1
178  block_dotc(inner_prod, tmp1, fine_vector, block_index);
179  for (int block_id = is; block_id < ns; block_id++) {
180  int index_re = index.idx_SPr(i, 1, block_id, 0);
181  int index_im = index.idx_SPi(i, 1, block_id, 0);
182  coarse_vector.set(index_re, 0.5 * real(inner_prod[block_id]));
183  coarse_vector.set(index_im, 0.5 * imag(inner_prod[block_id]));
184  }
185  } // i
186 #pragma omp barrier
187  }
188 
189 
190  //====================================================================
191  template<typename INDEX, class AFIELD1, class AFIELD2>
192  void make_fine_vector_impl(AFIELD2& fine_vector, const AFIELD1& coarse_vector,
193  const std::vector<AFIELD2>& testvectors,
194  AFopr<AFIELD2> *afopr, const INDEX& block_index,
195  AFIELD2& tmp1, AFIELD2& tmp2, typename AFIELD1::complex_t *complex_array)
196  {
197 #ifdef DEBUG
198  vout.detailed("%s is called\n", __func__);
199 #endif
200 
201  typedef typename AFIELD1::complex_t complex_t;
202  typedef typename AFIELD1::real_t real_t;
203 
204  const int num_vectors = testvectors.size();
205  assert(coarse_vector.nin() == 4 * num_vectors); // 4: 2 for complex, 2 for chirality
206 
207  int ith, nth, is, ns;
208  set_threadtask(ith, nth, is, ns, block_index.coarse_nvol());
209  int coarse_Nx = block_index.coarse_lattice_size(0);
210  int coarse_Ny = block_index.coarse_lattice_size(1);
211  int coarse_Nz = block_index.coarse_lattice_size(2);
212  int coarse_Nt = block_index.coarse_lattice_size(3);
213  AIndex_coarse_lex<real_t, AFIELD1::IMPL> index(coarse_Nx, coarse_Ny, coarse_Nz, coarse_Nt, num_vectors, 2);
214 
215 #pragma omp barrier
216 
217  // linear algebra on fine lattice
218 
219  fine_vector.set(0.0);
220 
221  for (int i = 0; i < num_vectors; i++) {
222  copy(tmp2, testvectors[i]);
223  afopr->mult_gm5(tmp1, tmp2);
224  axpy(tmp2, -1.0, tmp1); // tmp2=(1-gm5)v[i]
225  axpy(tmp1, +1.0, testvectors[i]); // tmp1=(1+gm5)v[i]
226 
227 #pragma omp barrier
228 
229  // linear algebra with bolcks
230  for (int block_id = is; block_id < ns; block_id++) { // ch=-1
231  int index_re = index.idx_SPr(i, 0, block_id, 0);
232  int index_im = index.idx_SPi(i, 0, block_id, 0);
233  real_t re = coarse_vector.cmp(index_re);
234  real_t im = coarse_vector.cmp(index_im);
235  complex_array[block_id] = complex_t(re, im);
236  }
237  block_axpy(fine_vector, &(complex_array[0]), tmp2, real_t(0.5),
238  block_index);
239 
240  for (int block_id = is; block_id < ns; block_id++) { // ch=+1
241  int index_re = index.idx_SPr(i, 1, block_id, 0);
242  int index_im = index.idx_SPi(i, 1, block_id, 0);
243  real_t re = coarse_vector.cmp(index_re);
244  real_t im = coarse_vector.cmp(index_im);
245  complex_array[block_id] = complex_t(re, im);
246  }
247  block_axpy(fine_vector, &(complex_array[0]), tmp1, real_t(0.5),
248  block_index);
249  }//i
250  }
251 
252 
253  //====================================================================
254  template<class AFIELD>
255  void set_testvectors_impl_global(std::vector<AFIELD>& testvectors, AFopr<AFIELD> *afopr, Field& vec_tmp)
256  {
257  vout.general("using global random numbers for the initial testvectors\n");
258 
259  typedef typename AFIELD::real_t real_t;
260  int nvec = testvectors.size();
261 
262  RandomNumbers *random = RandomNumberManager::getInstance();
263  for (int i = 0; i < nvec; i++) {
264  random->uniform_lex_global(vec_tmp);
265  //#pragma omp parallel
266  {
267  if (afopr->needs_convert()) {
268  if (i == 0) {
269  vout.detailed("convert in rep. required.\n");
270  }
271  afopr->convert(testvectors[i], vec_tmp);
272  } else {
273  if (i == 0) {
274  vout.detailed("convert in rep. not required.\n");
275  }
276  AIndex_lex<real_t, AFIELD::IMPL> index_f;
277  convert(index_f, testvectors[i], vec_tmp);
278  }
279  } // omp parallel
280 #pragma omp barrier
281  } // i
282  }
283 
284 
285  //====================================================================
286  template<class AFIELD>
287  void set_testvectors_impl_global(std::vector<AFIELD>& testvectors, AFopr<AFIELD> *afopr)
288  {
289  std::string class_name = __func__;
290  ThreadManager::assert_single_thread(class_name);
291 
292  int nin = testvectors[0].nin();
293  size_t nvol = testvectors[0].nvol();
294  Field vec_tmp(nin, nvol, 1);
295 #pragma omp parallel
296  {
297  set_testvectors_impl_global(testvectors, afopr, vec_tmp);
298  }
299  }
300 
301 
302  //====================================================================
303  template<class AFIELD>
304  void set_testvectors_impl_local(std::vector<AFIELD>& testvectors, AFopr<AFIELD> *afopr, Field& vec_tmp)
305  {
306  vout.general("using local random numbers for the initial testvectors\n");
307 
308  typedef typename AFIELD::real_t real_t;
309 
310  int nvec = testvectors.size();
311  int nin = testvectors[0].nin();
312  size_t nvol = testvectors[0].nvol();
313  assert(nin == vec_tmp.nin() && nvol == vec_tmp.nvol());
314  {
315  // const int seed_base=13957039; // charged pion: 139.57039 MeV
316  const int seed_base = 1395704; // charged pion: 139.57039 MeV
317  // const int seed_base=493677; // charged Kaon: 493.677 MeV
318  // const int seed_base=1234567;
319  int myrank = Communicator::nodeid();
320  int is0, ns0, ith, nth;
321  set_threadtask(ith, nth, is0, ns0, nvol);
322  int seed = seed_base + nth * myrank + ith;
323  RandomNumbers_Mseries random(seed);
324  size_t is = is0 * nin;
325  size_t ns = ns0 * nin;
326  for (int i = 0; i < nvec; i++) {
327  double *pv = vec_tmp.ptr(0);
328  for (size_t s = is; s < ns; s++) {
329  pv[s] = random.get();
330  }
331 #pragma omp barrier
332  if (afopr->needs_convert()) {
333  if (i == 0) {
334  vout.detailed("convert in rep. required.\n");
335  }
336  afopr->convert(testvectors[i], vec_tmp);
337  } else {
338  if (i == 0) {
339  vout.detailed("convert in rep. not required.\n");
340  }
341  AIndex_lex<real_t, AFIELD::IMPL> index_f;
342  convert(index_f, testvectors[i], vec_tmp);
343  }
344 #pragma omp barrier
345  } // i
346  }
347  }
348 
349 
350  //====================================================================
351  template<class AFIELD>
352  void set_testvectors_impl_local(std::vector<AFIELD>& testvectors, AFopr<AFIELD> *afopr)
353  {
354  std::string class_name = __func__;
355  ThreadManager::assert_single_thread(class_name);
356 
357  int nin = testvectors[0].nin();
358  size_t nvol = testvectors[0].nvol();
359  Field vec_tmp(nin, nvol, 1);
360 
361 #pragma omp parallel
362  {
363  set_testvectors_impl_local(testvectors, afopr, vec_tmp);
364  }
365  }
366 } // end of namespace
367 
368 //====================================================================
369 template<class AFIELD1, class AFIELD2>
370 void MultiGrid_Clover<AFIELD1, AFIELD2>::init_resources()
371 {
372  m_tmp1.reset(m_nin, m_fine_nvol, 1);
373  m_tmp2.reset(m_nin, m_fine_nvol, 1);
374  m_field_tmp.reset(m_nin, m_fine_nvol, 1);
375 
376  int coarse_nvol = m_block_index.coarse_nvol();
377  m_real_array.resize(coarse_nvol);
378  m_complex_array.resize(coarse_nvol);
379  m_coarse_array.resize(coarse_nvol * 2 * m_nvec);
380 }
381 
382 
383 //====================================================================
384 template<class AFIELD1, class AFIELD2>
385 void MultiGrid_Clover<AFIELD1, AFIELD2>::set_testvectors(const std::vector<AFIELD2>& w)
386 {
387  assert(w.size() == m_nvec);
388  for (int ivec = 0; ivec < m_nvec; ++ivec) {
389  copy(m_testvectors[ivec], w[ivec]);
390  }
391 }
392 
393 
394 template<class AFIELD1, class AFIELD2>
395 void MultiGrid_Clover<AFIELD1, AFIELD2>::set_testvectors()
396 {
397  Timer timer;
398  timer.start();
399 
400  //set_testvectors_impl_global(m_testvectors, m_afopr_fine);
401  set_testvectors_impl_local(m_testvectors, m_afopr_fine, m_field_tmp);
402 
403  timer.stop();
404  double elapsed_time = timer.elapsed_sec();
405  vout.detailed(m_vl, "%s: time for set_testvector = %e sec\n",
406  class_name.c_str(), elapsed_time);
407 }
408 
409 
410 //====================================================================
411 template<class AFIELD1, class AFIELD2>
412 void MultiGrid_Clover<AFIELD1, AFIELD2>::set_coarse_array(
413  const AFIELD1& coarse_vector) const
414 
415 {
416  typedef typename AFIELD1::real_t real_t;
417  typedef typename AFIELD1::complex_t complex_t;
418 
419  const int coarse_nvol = m_block_index.coarse_nvol();
420  const int num_vectors = m_testvectors.size();
421 
422  const int Nxc = m_block_index.coarse_lattice_size(0);
423  const int Nyc = m_block_index.coarse_lattice_size(1);
424  const int Nzc = m_block_index.coarse_lattice_size(2);
425  const int Ntc = m_block_index.coarse_lattice_size(3);
426 
427  AIndex_coarse_lex<real_t, AFIELD1::IMPL> index_c(
428  Nxc, Nyc, Nzc, Ntc, num_vectors, 2);
429 
430  int ith, nth, is_coarse, ns_coarse;
431  set_threadtask(ith, nth, is_coarse, ns_coarse, coarse_nvol);
432 
433 #pragma omp barrier
434 
435  for (int ivec = 0; ivec < num_vectors; ++ivec) {
436  for (int block = is_coarse; block < ns_coarse; ++block) {
437  int index_re1 = index_c.idx_SPr(ivec, 0, block, 0);
438  int index_im1 = index_c.idx_SPi(ivec, 0, block, 0);
439  real_t re1 = coarse_vector.cmp(index_re1);
440  real_t im1 = coarse_vector.cmp(index_im1);
441 
442  int index_re2 = index_c.idx_SPr(ivec, 1, block, 0);
443  int index_im2 = index_c.idx_SPi(ivec, 1, block, 0);
444  real_t re2 = coarse_vector.cmp(index_re2);
445  real_t im2 = coarse_vector.cmp(index_im2);
446 
447  // impl-2
448  // |f> += cm (1-gm5)/2 |i> + cp (1+gm5)/2 |i>
449  // = (cm+cp) |i>/2 + (-cm+cp) gm5|i>/2
450  // cm=(re1,im1), cp=(re2,im2)
451  m_coarse_array[block + coarse_nvol * (2 * ivec)]
452  = complex_t(re2 + re1, im2 + im1);
453  m_coarse_array[block + coarse_nvol * (2 * ivec + 1)]
454  = complex_t(re2 - re1, im2 - im1);
455 
456  /*
457  // impl-1
458  m_coarse_array[block + coarse_nvol * (2*ivec)]
459  = complex_t(re1, im1);
460  m_coarse_array[block + coarse_nvol * (2*ivec+1)]
461  = complex_t(re2, im2);
462  */
463  }
464  }
465 
466 #pragma omp barrier
467 }
468 
469 
470 //====================================================================
471 template<class AFIELD1, class AFIELD2>
472 void MultiGrid_Clover<AFIELD1, AFIELD2>::set_coarse_vector(
473  AFIELD1& coarse_vector) const
474 {
475  typedef typename AFIELD1::real_t real_t;
476  typedef typename AFIELD1::complex_t complex_t;
477 
478  const int coarse_nvol = m_block_index.coarse_nvol();
479  const int num_vectors = m_testvectors.size();
480 
481  const int Nxc = m_block_index.coarse_lattice_size(0);
482  const int Nyc = m_block_index.coarse_lattice_size(1);
483  const int Nzc = m_block_index.coarse_lattice_size(2);
484  const int Ntc = m_block_index.coarse_lattice_size(3);
485 
486  AIndex_coarse_lex<real_t, AFIELD1::IMPL> index_c(
487  Nxc, Nyc, Nzc, Ntc, num_vectors, 2);
488 
489  complex_t *array1 = &m_coarse_array[0];
490 
491  int ith, nth, is_coarse, ns_coarse;
492  set_threadtask(ith, nth, is_coarse, ns_coarse, coarse_nvol);
493 
494 #pragma omp barrier
495 
496  for (int ivec = 0; ivec < num_vectors; ++ivec) {
497  for (int block = is_coarse; block < ns_coarse; ++block) {
498  real_t re1 = real(array1[block + coarse_nvol * (2 * ivec)]);
499  real_t im1 = imag(array1[block + coarse_nvol * (2 * ivec)]);
500  int index_re1 = index_c.idx_SPr(ivec, 0, block, 0);
501  int index_im1 = index_c.idx_SPi(ivec, 0, block, 0);
502  coarse_vector.set(index_re1, re1);
503  coarse_vector.set(index_im1, im1);
504 
505  real_t re2 = real(array1[block + coarse_nvol * (2 * ivec + 1)]);
506  real_t im2 = imag(array1[block + coarse_nvol * (2 * ivec + 1)]);
507  int index_re2 = index_c.idx_SPr(ivec, 1, block, 0);
508  int index_im2 = index_c.idx_SPi(ivec, 1, block, 0);
509  coarse_vector.set(index_re2, re2);
510  coarse_vector.set(index_im2, im2);
511  }
512  }
513 #pragma omp barrier
514  coarse_vector.scal(real_t(0.5));
515 #pragma omp barrier
516 }
517 
518 
519 //====================================================================
520 template<class AFIELD1, class AFIELD2>
521 void MultiGrid_Clover<AFIELD1, AFIELD2>::make_fine_vector(AFIELD2& fine_vector, const AFIELD1& coarse_vector) const
522 {
523  vout.general("%s: make_fine_vector is called\n", class_name.c_str());
524 
525  // fine_vec
526  // = coarse[ic, ch=-1, block] * P_- testvecr[ic, block]
527  // + coarse[ic, ch=+1, block]* P_+ testvec[ic, block]
528  // = (coarse[ic, ch=-1, block] + coarse[ic, ch=+1, block])
529  // *testvec[ic, block] / 2
530  // + (- coarse[ic, ch=-1, block] + coarse[ic, ch=+1, block])
531  // * gamma5 * testvec[ic, block] /2
532 
533  // unique_ptr<Timer> timer(new Timer);
534  // timer->start();
535 
536  // make_fine_vector_impl(fine_vector, coarse_vector,
537  // m_testvectors, m_afopr_fine, m_block_index,
538  // m_tmp1, m_tmp2, &(m_complex_array[0]));
539 
540  typedef typename AFIELD1::complex_t complex_t;
541  typedef typename AFIELD1::real_t real_t;
542 
543  const int num_vectors = m_testvectors.size();
544  assert(coarse_vector.nin() == 4 * num_vectors);
545  // 4: 2 for complex * 2 for chirality
546 
547  const int coarse_nvol = m_block_index.coarse_nvol();
548 
549  set_coarse_array(coarse_vector);
550 
551 #pragma omp barrier
552 
553  fine_vector.set(0.0);
554 
555  for (int ivec = 0; ivec < num_vectors; ++ivec) {
556  // imple-2
557  m_afopr_fine->mult_gm5(m_tmp1, m_testvectors[ivec]);
558 
559  // impl-1
560  //copy(m_tmp2, m_testvectors[ivec]);
561  //m_afopr_fine->mult_gm5(m_tmp1, m_tmp2);
562  //axpy(m_tmp2, -1.0, m_tmp1); // tmp2 = (1 - gm5) v[i]
563  //axpy(m_tmp1, +1.0, m_testvectors[ivec]); // tmp1 = (1 + gm5) v[i]
564 
565 #pragma omp barrier
566 
567  complex_t *array1 = &m_coarse_array[coarse_nvol * (2 * ivec)];
568 
569  // impl-2
570  block_axpy(fine_vector, array1, m_testvectors[ivec], real_t(0.5),
571  m_block_index);
572  // impl-1
573  //block_axpy(fine_vector, array1, m_tmp2, real_t(0.5),
574  // m_block_index);
575 
576  complex_t *array2 = &m_coarse_array[coarse_nvol * (2 * ivec + 1)];
577  block_axpy(fine_vector, array2, m_tmp1, real_t(0.5),
578  m_block_index);
579  }
580 
581  // timer->stop();
582  // double elapsed_time = timer->elapsed_sec();
583  // vout.general(" time for make_fine_vector = %e sec\n", elapsed_time);
584 }
585 
586 
587 //====================================================================
588 template<class AFIELD1, class AFIELD2>
589 void MultiGrid_Clover<AFIELD1, AFIELD2>::make_coarse_vector(AFIELD1& coarse_vector, const AFIELD2& fine_vector) const
590 {
591  vout.general("%s: make_coarse_vector is called\n", class_name.c_str());
592 
593  typedef typename AFIELD1::complex_t complex_t;
594  typedef typename AFIELD1::real_t real_t;
595 
596  // unique_ptr<Timer> timer(new Timer);
597  // timer->start();
598 
599  // make_coarse_vector_impl(coarse_vector, fine_vector,
600  // m_testvectors, m_afopr_fine, m_block_index,
601  // m_tmp1, m_tmp2, &(m_complex_array[0]));
602 
603 
604  const int num_vectors = m_testvectors.size();
605  assert(coarse_vector.nin() == 4 * num_vectors);
606  // 4: 2 for complex, 2 for chirality
607 
608  const int coarse_nvol = m_block_index.coarse_nvol();
609 
610  for (int ivec = 0; ivec < num_vectors; ++ivec) {
611  copy(m_tmp1, m_testvectors[ivec]);
612  m_afopr_fine->mult_gm5(m_tmp2, m_testvectors[ivec]);
613 #pragma omp barrier
614  axpy(m_tmp1, -1.0, m_tmp2); // tmp1 = (1 - gm5) v[i]
615  axpy(m_tmp2, +1.0, m_testvectors[ivec]); // tmp2 = (1 + gm5) v[i]
616 #pragma omp barrier
617 
618  // ch=-1
619  complex_t *array1 = &m_coarse_array[coarse_nvol * (2 * ivec)];
620  block_dotc(array1, m_tmp1, fine_vector, m_block_index);
621 
622  // ch=+1
623  complex_t *array2 = &m_coarse_array[coarse_nvol * (2 * ivec + 1)];
624  block_dotc(array2, m_tmp2, fine_vector, m_block_index);
625 
626 #pragma omp barrier
627  } // ivec
628 
629  set_coarse_vector(coarse_vector);
630 
631  // timer->stop();
632  // double elapsed_time = timer->elapsed_sec();
633  // vout.general(" time for make_coarse_vector = %e sec\n", elapsed_time);
634 }
635 
636 
637 //====================================================================
638 template<class AFIELD1, class AFIELD2>
639 void MultiGrid_Clover<AFIELD1, AFIELD2>::gramschmidt()
640 {
641  gramschmidt(m_testvectors);
642 }
643 
644 
645 //====================================================================
646 template<class AFIELD1, class AFIELD2>
647 void MultiGrid_Clover<AFIELD1, AFIELD2>::gramschmidt(
648  std::vector<AFIELD2>& vectors) const
649 {
650  vout.general("%s: ::gramschmidt is called\n", class_name.c_str());
651 
652  gramschmidt_in_block_impl(vectors, m_afopr_fine, m_block_index,
653  m_tmp1, m_tmp2,
654  &(m_real_array[0]), &(m_complex_array[0]));
655 
656  gramschmidt_in_block_impl(vectors, m_afopr_fine, m_block_index,
657  m_tmp1, m_tmp2,
658  &(m_real_array[0]), &(m_complex_array[0]));
659 }
660 
661 
662 //====================================================================
663 template<class AFIELD1, class AFIELD2>
664 void MultiGrid_Clover<AFIELD1, AFIELD2>::gramschmidt_double(
665  std::vector<AFIELD2>& vectors) const
666 {
667  vout.crucial("%s: gramschmed_double is called\n", class_name.c_str());
668  abort();
669 }
670 
671 
672 //============================================================END=====
MultiGrid_Clover::set_coarse_array
void set_coarse_array(const AFIELD1 &coarse_vector) const
array <- vector (function for optimization)
Definition: MultiGrid_Clover-tmpl.h:412
randomNumbers_Mseries.h
AFopr
Definition: afopr.h:48
AIndex_lex
Definition: aindex_lex_base.h:17
convert
void convert(INDEX &index, AFIELD &v, const Field &w)
Definition: afield-inc.h:129
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:219
RandomNumbers
Base class of random number generators.
Definition: randomNumbers.h:43
AFopr::convert
virtual void convert(AFIELD &, const Field &)
converts a Field object into other format if necessary.
Definition: afopr.h:177
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:380
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
AFopr::mult_gm5
virtual void mult_gm5(AFIELD &, const AFIELD &)
multiplies gamma_5 matrix.
Definition: afopr.h:123
Timer
Definition: timer.h:31
MultiGrid_Clover::make_fine_vector
void make_fine_vector(AFIELD2 &fine_vector, const AFIELD1 &coarse_vector) const
Definition: MultiGrid_Clover-tmpl.h:521
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:212
Field::real_t
double real_t
Definition: field.h:51
AIndex_coarse_lex
Definition: aindex_coarse_lex_base.h:17
block_dotc
void block_dotc(typename AFIELD::complex_t *out, const AFIELD &v, const AFIELD &w, const INDEX &block_index)
Definition: afield_dd-inc.h:24
Timer::start
void start()
Definition: timer.cpp:44
AFopr::needs_convert
virtual bool needs_convert()
returns true if additional field conversion is needed.
Definition: afopr.h:174
MultiGrid_Clover::set_testvectors
void set_testvectors()
Definition: MultiGrid_Clover-tmpl.h:395
timer.h
MultiGrid_Clover::gramschmidt_double
void gramschmidt_double(std::vector< AFIELD2 > &fine_vectors) const
Definition: MultiGrid_Clover-tmpl.h:664
RandomNumbers::get
virtual double get()=0
RandomNumbers_Mseries
Random number generator base on M-series.
Definition: randomNumbers_Mseries.h:46
MultiGrid_Clover::gramschmidt
void gramschmidt()
Definition: MultiGrid_Clover-tmpl.h:639
Field::nvol
int nvol() const
Definition: field.h:127
MultiGrid_Clover::set_coarse_vector
void set_coarse_vector(AFIELD1 &coarse_vector) const
vector <- array (function for optimization)
Definition: MultiGrid_Clover-tmpl.h:472
RandomNumberManager::getInstance
static RandomNumbers * getInstance()
Definition: randomNumberManager.cpp:45
MultiGrid_Clover::init_resources
void init_resources()
Definition: MultiGrid_Clover-tmpl.h:370
block_scal
void block_scal(AFIELD &v, const typename AFIELD::real_t *a, const INDEX &block_index)
Definition: afield_dd-inc.h:221
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
randomNumberManager.h
RandomNumbers::uniform_lex_global
virtual void uniform_lex_global(Field &)
uniform random number defined on global lattice.
Definition: randomNumbers.cpp:106
complex_t
ComplexTraits< double >::complex_t complex_t
Definition: afopr_Clover_coarse_double.cpp:23
block_norm2
void block_norm2(typename AFIELD::real_t *out, const AFIELD &v, const INDEX &block_index)
Definition: afield_dd-inc.h:132
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:180
MultiGrid_Clover::make_coarse_vector
void make_coarse_vector(AFIELD1 &coarse_vector, const AFIELD2 &fine_vector) const
Definition: MultiGrid_Clover-tmpl.h:589
Field
Container of Field-type object.
Definition: field.h:46
Timer::elapsed_sec
double elapsed_sec() const
Definition: timer.cpp:107
block_axpy
void block_axpy(AFIELD &v, const typename AFIELD::real_t *a, const AFIELD &w, const typename AFIELD::real_t fac, const INDEX &block_index)
Definition: afield_dd-inc.h:388
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
Timer::stop
void stop()
Definition: timer.cpp:69
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:512