LatticeYangMills
lattice.h
Go to the documentation of this file.
1 /******************************************************************************
2 *
3 * MIT License
4 *
5 * Copyright (c) 2018 Giovanni Pederiva
6 *
7 * Permission is hereby granted, free of charge, to any person obtaining a copy
8 * of this software and associated documentation files (the "Software"), to deal
9 * in the Software without restriction, including without limitation the rights
10 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
11 * copies of the Software, and to permit persons to whom the Software is
12 * furnished to do so, subject to the following conditions:
13 *
14 * The above copyright notice and this permission notice shall be included in
15 * all copies or substantial portions of the Software.
16 *
17 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
18 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
19 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
20 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
21 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
22 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
23 * SOFTWARE.
24 ******************************************************************************/
25 
34 #pragma once
35 #include <vector>
36 #include <array>
37 #include <utility>
38 #include "su3.h"
39 #include "Utils/clusterspecifier.h"
40 #ifndef LACONIA
41 #include <mpi/mpi.h>
42 #else
43 #include <mpi.h>
44 #endif
45 #include "ParallelTools/parallel.h"
46 
61 template <typename T>
62 class Lattice{
63 public:
64  // Default Constructors
65  Lattice();
66  // Copy Constructors
67  Lattice(const Lattice& other) noexcept{
68  allocate(other.size);
69  for(int i = 0; i < sites; i++)
70  lattice[i] = other.lattice[i];
71  }
72  Lattice(Lattice&& other) noexcept :
73  lattice(std::move(other.lattice)),
74  sites(std::move(other.sites)),
75  size(std::move(other.size)) { }
76 
77  // Array-based constructor
78  Lattice(std::array<int, 4> sizeArray);
79  void allocate(std::array<int, 4> sizeArray);
80 
81  // Data
82  std::vector<T> lattice;
83  std::array<int,4> size;
84  int sites;
85 
86  // Element accessors
87  inline T& at(int x, int y, int z, int t) {
88  return lattice[x + size[0]*(y + size[1]*(z + size[2]*t))];
89  }
90 
91  inline T& at(const std::vector<int>& site){
92  return lattice[site[0] + size[0]*(site[1] + size[1]*(site[2] + size[2]*site[3]))];
93  }
94  inline T& at(const std::array<int, 4>& site){
95  return lattice[site[0] + size[0]*(site[1] + size[1]*(site[2] + size[2]*site[3]))];
96  }
97  inline T& at(int i){
98  return lattice[i];
99  }
100  inline const T& at(int x, int y, int z, int t) const {
101  return lattice[x + size[0]*(y + size[1]*(z + size[2]*t))];
102  }
103 
104  inline const T& at(const std::vector<int>& site) const {
105  return lattice[site[0] + size[0]*(site[1] + size[1]*(site[2] + size[2]*site[3]))];
106  }
107  inline const T& at(const std::array<int, 4>& site) const {
108  return lattice[site[0] + size[0]*(site[1] + size[1]*(site[2] + size[2]*site[3]))];
109  }
110  inline const T& at(int i) const {
111  return lattice[i];
112  }
113 
114  // Copy-assignment
115  inline Lattice& operator=(const Lattice& other) noexcept{
116  lattice = other.lattice;
117  sites = other.sites;
118  size = other.size;
119  return *this;
120  }
121  inline Lattice& operator=(Lattice&& other) noexcept{
122  std::swap(lattice, other.lattice);
123  std::swap(sites, other.sites);
124  std::swap(size, other.size);
125  return *this;
126  }
127 
128  // OPERATIOSN BETWEEN LATTICES
129  // Sum
130  inline Lattice& operator+=( const Lattice& other ) noexcept{
131  for(int i = 0; i < sites; i++)
132  lattice[i] += other.lattice[i];
133  return *this;
134  }
135  inline Lattice& operator+=( Lattice&& other ) noexcept{
136  for(int i = 0; i < sites; i++)
137  lattice[i] += other.lattice[i];
138  return *this;
139  }
140 
141  friend Lattice operator+(Lattice lhs, const Lattice& rhs) noexcept{
142  lhs += rhs;
143  return lhs;
144  }
145 
146  friend Lattice operator+(Lattice lhs, Lattice&& rhs) noexcept{
147  lhs += rhs;
148  return lhs;
149  }
150 
151  // Subtraction
152  inline Lattice& operator-=( const Lattice& other ) noexcept{
153  for(int i = 0; i < sites; i++)
154  lattice[i] -= other.lattice[i];
155  return *this;
156  }
157  inline Lattice& operator-=( Lattice&& other ) noexcept{
158  for(int i = 0; i < sites; i++)
159  lattice[i] -= other.lattice[i];
160  return *this;
161  }
162 
163  friend Lattice operator-(Lattice lhs, const Lattice& rhs) noexcept{
164  lhs -= rhs;
165  return lhs;
166  }
167 
168  friend Lattice operator-(Lattice lhs, Lattice&& rhs) noexcept{
169  lhs -= rhs;
170  return lhs;
171  }
172 
173  // Multiplication
174  inline Lattice& operator*=( const Lattice& other ) noexcept{
175  for(int i = 0; i < sites; i++)
176  lattice[i] *= other.lattice[i];
177  return *this;
178  }
179  inline Lattice& operator*=( Lattice&& other ) noexcept{
180  for(int i = 0; i < sites; i++)
181  lattice[i] *= other.lattice[i];
182  return *this;
183  }
184 
185  friend Lattice operator*(Lattice lhs, const Lattice& rhs) noexcept{
186  lhs *= rhs;
187  return lhs;
188  }
189 
190  friend Lattice operator*(Lattice lhs, Lattice&& rhs) noexcept{
191  lhs *= rhs;
192  return lhs;
193  }
194 
195 
196  // OPERATIONS WITH A SCALAR
197  // Sum
198  inline Lattice& operator+=( double scalar ) noexcept {
199  for(int i = 0; i < sites; i++)
200  lattice[i] += scalar;
201  return *this;
202  }
203  friend Lattice operator+(Lattice lhs, double scalar) noexcept {
204  lhs += scalar;
205  return lhs;
206  }
207 
208  // Subtraction
209  inline Lattice& operator-=( double scalar ) noexcept {
210  for(int i = 0; i < sites; i++)
211  lattice[i] -= scalar;
212  return *this;
213  }
214  friend Lattice operator-(Lattice lhs, double scalar) noexcept {
215  lhs -= scalar;
216  return lhs;
217  }
218 
219  // Multiplication
220  inline Lattice& operator*=( double scalar ) noexcept {
221  for(int i = 0; i < sites; i++)
222  lattice[i] *= scalar;
223  return *this;
224  }
225  friend Lattice operator*(Lattice lhs, double scalar) noexcept {
226  lhs *= scalar;
227  return lhs;
228  }
229 };
230 
231 template <typename T>
232 inline Lattice<T> adj(const Lattice<T> &base){
233  Lattice<T> result(base.size);
234  for(int i = 0; i < result.sites; i++)
235  result.at(i) = ~base.at(i);
236  return std::move(result);
237 }
238 
239 template <typename T>
240 inline Lattice<T> adj(Lattice<T>&& base){
241  Lattice<T> result(base.size);
242  for(int i = 0; i < result.sites; i++)
243  result.at(i) = ~base.at(i);
244  return std::move(result);
245 }
246 
247 template <typename T>
248 inline Lattice<T> exp(const Lattice<T> &base){
249  Lattice<T> result(base.size);
250  for(int i = 0; i < result.sites; i++)
251  result.at(i) = exp(base.at(i));
252  return std::move(result);
253 }
254 
255 template <typename T>
256 inline Lattice<T> exp(Lattice<T>&& base){
257  Lattice<T> result(base.size);
258  for(int i = 0; i < result.sites; i++)
259  result.at(i) = exp(base.at(i));
260  return std::move(result);
261 }
262 
263 
264 template <typename T>
265 inline T sum(Lattice<T>& base){
266  T result = 0;
267  for(int i = 0; i < base.sites; i++)
268  result += base.at(i);
269  return result;
270 }
271 
272 template <>
273 inline SU3 sum(Lattice<SU3>& base){
274  SU3 result;
275  result.setSU3Zero();
276  for(int i = 0; i < base.sites; i++)
277  result += base.at(i);
278  return result;
279 }
280 
281 template <typename T>
282 inline T sum(Lattice<T>&& base){
283  T result = 0;
284  for(int i = 0; i < base.sites; i++)
285  result += base.at(i);
286  return result;
287 }
288 
289 template <>
290 inline SU3 sum(Lattice<SU3>&& base){
291  SU3 result;
292  result.setSU3Zero();
293  for(int i = 0; i < base.sites; i++)
294  result += base.at(i);
295  return std::move(result);
296 }
297 
298 inline Lattice<double> realTrace(Lattice<SU3> &base){
299  Lattice<double> result(base.size);
300  for(int i = 0; i < base.sites; i++)
301  result.at(i) = base.at(i).realTrace();
302  return std::move(result);
303 }
304 
305 inline Lattice<double> imagTrace(Lattice<SU3> &base){
306  Lattice<double> result(base.size);
307  for(int i = 0; i < base.sites; i++)
308  result.at(i) = base.at(i).imagTrace();
309  return std::move(result);
310 }
311 
312 inline Lattice<double> realTrace(Lattice<SU3> &&base){
313  Lattice<double> result(base.size);
314  for(int i = 0; i < base.sites; i++)
315  result.at(i) = base.at(i).realTrace();
316  return std::move(result);
317 }
318 
319 inline Lattice<double> imagTrace(Lattice<SU3> &&base){
320  Lattice<double> result(base.size);
321  for(int i = 0; i < base.sites; i++)
322  result.at(i) = base.at(i).imagTrace();
323  return std::move(result);
324 }
325 
326 
327 template <typename T>
328 Lattice<T>::Lattice(){
329 
330 }
331 
332 
333 template <typename T>
334 Lattice<T>::Lattice(std::array<int, 4> sizeArray){
335  allocate(sizeArray);
336 }
337 
338 template <typename T>
339 void Lattice<T>::allocate(std::array<int, 4> sizeArray){
340  // Allocate memory
341  size = sizeArray;
342  sites = size[0] * size[1] * size[2] * size[3];
343  lattice.resize(sites);
344 }
345 
346 
347 void setLatticeImagIdentityValue(Lattice<SU3>& su3lat, const Lattice<double>& lat);
348 void setLatticeImagIdentityValue(Lattice<SU3>& su3lat, Lattice<double>&& lat);
349 void setToZero(Lattice<SU3>& su3lat);
350 void setToIdentity(Lattice<SU3>& su3lat);
351 
352 
353 
354 
355 
356 
357 
358 
359 
360 template<typename T>
361 Lattice<T> shift(const Lattice<T>& lat, int shiftDir, int shiftStep){
362  // Create Shifted Matrix
363  Lattice<T> shifted;
364  shifted.allocate(lat.size);
365 
366  // Create send and recv buffers
367  std::array<int,4> bufferSize;
368  Lattice<T> sendBuffer, recvBuffer;
369  for(int i = 0; i < 4; i++)
370  i != shiftDir ? bufferSize[i] = shifted.size[i] : bufferSize[i] = 1;
371  sendBuffer.allocate(bufferSize);
372  recvBuffer.allocate(bufferSize);
373 
374  // Fill send buffer
375  // X
376  if(shiftDir == 0 && shiftStep == 1){
377  for(int y = 0; y < lat.size[1]; y++){
378  for(int z = 0; z < lat.size[2]; z++){
379  for(int t = 0; t < lat.size[3]; t++){
380  sendBuffer.at(0,y,z,t) = lat.at(0,y,z,t);
381  }}}
382  }
383  else if(shiftDir == 0 && shiftStep == -1){
384  for(int y = 0; y < lat.size[1]; y++){
385  for(int z = 0; z < lat.size[2]; z++){
386  for(int t = 0; t < lat.size[3]; t++){
387  sendBuffer.at(0,y,z,t) = lat.at(lat.size[0] - 1,y,z,t);
388  }}}
389  }
390 
391  // Y
392  else if(shiftDir == 1 && shiftStep == 1){
393  for(int x = 0; x < lat.size[0]; x++){
394  for(int z = 0; z < lat.size[2]; z++){
395  for(int t = 0; t < lat.size[3]; t++){
396  sendBuffer.at(x,0,z,t) = lat.at(x,0,z,t);
397  }}}
398  }
399  else if(shiftDir == 1 && shiftStep == -1){
400  for(int x = 0; x < lat.size[0]; x++){
401  for(int z = 0; z < lat.size[2]; z++){
402  for(int t = 0; t < lat.size[3]; t++){
403  sendBuffer.at(x,0,z,t) = lat.at(x,lat.size[1] - 1,z,t);
404  }}}
405  }
406 
407  // Z
408  else if(shiftDir == 2 && shiftStep == 1){
409  for(int x = 0; x < lat.size[0]; x++){
410  for(int y = 0; y < lat.size[1]; y++){
411  for(int t = 0; t < lat.size[3]; t++){
412  sendBuffer.at(x,y,0,t) = lat.at(x,y,0,t);
413  }}}
414  }
415  else if(shiftDir == 2 && shiftStep == -1){
416  for(int x = 0; x < lat.size[0]; x++){
417  for(int y = 0; y < lat.size[1]; y++){
418  for(int t = 0; t < lat.size[3]; t++){
419  sendBuffer.at(x,y,0,t) = lat.at(x,y,lat.size[2] - 1,t);
420  }}}
421  }
422  else if(shiftDir == 3 && shiftStep == 1){
423  for(int x = 0; x < lat.size[0]; x++){
424  for(int y = 0; y < lat.size[1]; y++){
425  for(int z = 0; z < lat.size[2]; z++){
426  sendBuffer.at(x,y,z,0) = lat.at(x,y,z,0);
427  }}}
428  }
429  else if(shiftDir == 3 && shiftStep == -1){
430  for(int x = 0; x < lat.size[0]; x++){
431  for(int y = 0; y < lat.size[1]; y++){
432  for(int z = 0; z < lat.size[2]; z++){
433  sendBuffer.at(x,y,z,0) = lat.at(x,y,z,lat.size[3] - 1);
434  }}}
435  }
436 
437  // Send non-blockingly the msg
438  MPI_Request sendRequest, recvRequest;
439  MPI_Status status;
440  int bufferCount = 18*bufferSize[0]*bufferSize[1]*bufferSize[2]*bufferSize[3];
441  int shiftRank;
442  shiftStep == -1 ? shiftRank = 0 : shiftRank = 1;
443 
444  MPI_Isend(sendBuffer.lattice.data(), bufferCount, MPI_DOUBLE,
445  Parallel::getNeighbor(shiftDir, abs(shiftRank-1)), 0,
446  MPI_COMM_WORLD, &sendRequest);
447  MPI_Irecv(recvBuffer.lattice.data(), bufferCount, MPI_DOUBLE,
448  Parallel::getNeighbor(shiftDir, shiftRank), 0,
449  MPI_COMM_WORLD, &recvRequest);
450 
451 
452  // Shift the internal points
453  // X
454  if(shiftDir == 0 && shiftStep == 1){
455  for(int x = 0; x < lat.size[0]-1; x++){
456  for(int y = 0; y < lat.size[1]; y++){
457  for(int z = 0; z < lat.size[2]; z++){
458  for(int t = 0; t < lat.size[3]; t++){
459  shifted.at(x,y,z,t) = lat.at(x+1,y,z,t);
460  }}}}
461  }
462  else if(shiftDir == 0 && shiftStep == -1){
463  for(int x = 1; x < lat.size[0]; x++){
464  for(int y = 0; y < lat.size[1]; y++){
465  for(int z = 0; z < lat.size[2]; z++){
466  for(int t = 0; t < lat.size[3]; t++){
467  shifted.at(x,y,z,t) = lat.at(x-1,y,z,t);
468  }}}}
469  }
470  // Y
471  else if(shiftDir == 1 && shiftStep == 1){
472  for(int x = 0; x < lat.size[0]; x++){
473  for(int y = 0; y < lat.size[1]-1; y++){
474  for(int z = 0; z < lat.size[2]; z++){
475  for(int t = 0; t < lat.size[3]; t++){
476  shifted.at(x,y,z,t) = lat.at(x,y+1,z,t);
477  }}}}
478  }
479  else if(shiftDir == 1 && shiftStep == -1){
480  for(int x = 0; x < lat.size[0]; x++){
481  for(int y = 1; y < lat.size[1]; y++){
482  for(int z = 0; z < lat.size[2]; z++){
483  for(int t = 0; t < lat.size[3]; t++){
484  shifted.at(x,y,z,t) = lat.at(x,y-1,z,t);
485  }}}}
486  }
487  // Z
488  else if(shiftDir == 2 && shiftStep == 1){
489  for(int x = 0; x < lat.size[0]; x++){
490  for(int y = 0; y < lat.size[1]; y++){
491  for(int z = 0; z < lat.size[2]-1; z++){
492  for(int t = 0; t < lat.size[3]; t++){
493  shifted.at(x,y,z,t) = lat.at(x,y,z+1,t);
494  }}}}
495  }
496  else if(shiftDir == 2 && shiftStep == -1){
497  for(int x = 0; x < lat.size[0]; x++){
498  for(int y = 0; y < lat.size[1]; y++){
499  for(int z = 1; z < lat.size[2]; z++){
500  for(int t = 0; t < lat.size[3]; t++){
501  shifted.at(x,y,z,t) = lat.at(x,y,z-1,t);
502  }}}}
503  }
504  // T
505  else if(shiftDir == 3 && shiftStep == 1){
506  for(int x = 0; x < lat.size[0]; x++){
507  for(int y = 0; y < lat.size[1]; y++){
508  for(int z = 0; z < lat.size[2]; z++){
509  for(int t = 0; t < lat.size[3]-1; t++){
510  shifted.at(x,y,z,t) = lat.at(x,y,z,t+1);
511  }}}}
512  }
513  else if(shiftDir == 3 && shiftStep == -1){
514  for(int x = 0; x < lat.size[0]; x++){
515  for(int y = 0; y < lat.size[1]; y++){
516  for(int z = 0; z < lat.size[2]; z++){
517  for(int t = 1; t < lat.size[3]; t++){
518  shifted.at(x,y,z,t) = lat.at(x,y,z,t-1);
519  }}}}
520  }
521 
522 
523 
524  // Wait for the message to arrive
525  MPI_Wait(&sendRequest, &status);
526  MPI_Wait(&recvRequest, &status);
527 
528 
529  // Complete the result with the received data
530  // X
531  if(shiftDir == 0 && shiftStep == 1){
532  for(int y = 0; y < lat.size[1]; y++){
533  for(int z = 0; z < lat.size[2]; z++){
534  for(int t = 0; t < lat.size[3]; t++){
535  shifted.at(lat.size[0] - 1,y,z,t) = recvBuffer.at(0,y,z,t);
536  }}}
537  }
538  else if(shiftDir == 0 && shiftStep == -1){
539  for(int y = 0; y < lat.size[1]; y++){
540  for(int z = 0; z < lat.size[2]; z++){
541  for(int t = 0; t < lat.size[3]; t++){
542  shifted.at(0,y,z,t) = recvBuffer.at(0,y,z,t);
543  }}}
544  }
545 
546  // Y
547  else if(shiftDir == 1 && shiftStep == 1){
548  for(int x = 0; x < lat.size[0]; x++){
549  for(int z = 0; z < lat.size[2]; z++){
550  for(int t = 0; t < lat.size[3]; t++){
551  shifted.at(x,lat.size[1] - 1,z,t) = recvBuffer.at(x,0,z,t);
552  }}}
553  }
554  else if(shiftDir == 1 && shiftStep == -1){
555  for(int x = 0; x < lat.size[0]; x++){
556  for(int z = 0; z < lat.size[2]; z++){
557  for(int t = 0; t < lat.size[3]; t++){
558  shifted.at(x,0,z,t) = recvBuffer.at(x,0,z,t);
559  }}}
560  }
561 
562  // Z
563  else if(shiftDir == 2 && shiftStep == 1){
564  for(int x = 0; x < lat.size[0]; x++){
565  for(int y = 0; y < lat.size[1]; y++){
566  for(int t = 0; t < lat.size[3]; t++){
567  shifted.at(x,y,lat.size[2] - 1,t) = recvBuffer.at(x,y,0,t);
568  }}}
569  }
570  else if(shiftDir == 2 && shiftStep == -1){
571  for(int x = 0; x < lat.size[0]; x++){
572  for(int y = 0; y < lat.size[1]; y++){
573  for(int t = 0; t < lat.size[3]; t++){
574  shifted.at(x,y,0,t) = recvBuffer.at(x,y,0,t);
575  }}}
576  }
577  else if(shiftDir == 3 && shiftStep == 1){
578  for(int x = 0; x < lat.size[0]; x++){
579  for(int y = 0; y < lat.size[1]; y++){
580  for(int z = 0; z < lat.size[2]; z++){
581  shifted.at(x,y,z,lat.size[3] - 1) = recvBuffer.at(x,y,z,0);
582  }}}
583  }
584  else if(shiftDir == 3 && shiftStep == -1){
585  for(int x = 0; x < lat.size[0]; x++){
586  for(int y = 0; y < lat.size[1]; y++){
587  for(int z = 0; z < lat.size[2]; z++){
588  shifted.at(x,y,z,0) = recvBuffer.at(x,y,z,0);
589  }}}
590  }
591 
592  // Return
593  return shifted;
594 }
595 
596 
597 template<typename T>
598 Lattice<T> shift(Lattice<T>&& lat, int shiftDir, int shiftStep){
599  // Create Shifted Matrix
600  Lattice<T> shifted;
601  shifted.allocate(lat.size);
602 
603  // Create send and recv buffers
604  std::array<int,4> bufferSize;
605  Lattice<T> sendBuffer, recvBuffer;
606  for(int i = 0; i < 4; i++)
607  i != shiftDir ? bufferSize[i] = shifted.size[i] : bufferSize[i] = 1;
608  sendBuffer.allocate(bufferSize);
609  recvBuffer.allocate(bufferSize);
610 
611  // Fill send buffer
612  // X
613  if(shiftDir == 0 && shiftStep == 1){
614  for(int y = 0; y < lat.size[1]; y++){
615  for(int z = 0; z < lat.size[2]; z++){
616  for(int t = 0; t < lat.size[3]; t++){
617  sendBuffer.at(0,y,z,t) = lat.at(0,y,z,t);
618  }}}
619  }
620  else if(shiftDir == 0 && shiftStep == -1){
621  for(int y = 0; y < lat.size[1]; y++){
622  for(int z = 0; z < lat.size[2]; z++){
623  for(int t = 0; t < lat.size[3]; t++){
624  sendBuffer.at(0,y,z,t) = lat.at(lat.size[0] - 1,y,z,t);
625  }}}
626  }
627 
628  // Y
629  else if(shiftDir == 1 && shiftStep == 1){
630  for(int x = 0; x < lat.size[0]; x++){
631  for(int z = 0; z < lat.size[2]; z++){
632  for(int t = 0; t < lat.size[3]; t++){
633  sendBuffer.at(x,0,z,t) = lat.at(x,0,z,t);
634  }}}
635  }
636  else if(shiftDir == 1 && shiftStep == -1){
637  for(int x = 0; x < lat.size[0]; x++){
638  for(int z = 0; z < lat.size[2]; z++){
639  for(int t = 0; t < lat.size[3]; t++){
640  sendBuffer.at(x,0,z,t) = lat.at(x,lat.size[1] - 1,z,t);
641  }}}
642  }
643 
644  // Z
645  else if(shiftDir == 2 && shiftStep == 1){
646  for(int x = 0; x < lat.size[0]; x++){
647  for(int y = 0; y < lat.size[1]; y++){
648  for(int t = 0; t < lat.size[3]; t++){
649  sendBuffer.at(x,y,0,t) = lat.at(x,y,0,t);
650  }}}
651  }
652  else if(shiftDir == 2 && shiftStep == -1){
653  for(int x = 0; x < lat.size[0]; x++){
654  for(int y = 0; y < lat.size[1]; y++){
655  for(int t = 0; t < lat.size[3]; t++){
656  sendBuffer.at(x,y,0,t) = lat.at(x,y,lat.size[2] - 1,t);
657  }}}
658  }
659  else if(shiftDir == 3 && shiftStep == 1){
660  for(int x = 0; x < lat.size[0]; x++){
661  for(int y = 0; y < lat.size[1]; y++){
662  for(int z = 0; z < lat.size[2]; z++){
663  sendBuffer.at(x,y,z,0) = lat.at(x,y,z,0);
664  }}}
665  }
666  else if(shiftDir == 3 && shiftStep == -1){
667  for(int x = 0; x < lat.size[0]; x++){
668  for(int y = 0; y < lat.size[1]; y++){
669  for(int z = 0; z < lat.size[2]; z++){
670  sendBuffer.at(x,y,z,0) = lat.at(x,y,z,lat.size[3] - 1);
671  }}}
672  }
673 
674  // Send non-blockingly the msg
675  MPI_Request sendRequest, recvRequest;
676  MPI_Status status;
677  int bufferCount = 18*bufferSize[0]*bufferSize[1]*bufferSize[2]*bufferSize[3];
678  int shiftRank;
679  shiftStep == -1 ? shiftRank = 0 : shiftRank = 1;
680 
681  MPI_Isend(sendBuffer.lattice.data(), bufferCount, MPI_DOUBLE,
682  Parallel::getNeighbor(shiftDir, abs(shiftRank-1)), 0,
683  MPI_COMM_WORLD, &sendRequest);
684  MPI_Irecv(recvBuffer.lattice.data(), bufferCount, MPI_DOUBLE,
685  Parallel::getNeighbor(shiftDir, shiftRank), 0,
686  MPI_COMM_WORLD, &recvRequest);
687 
688 
689  // Shift the internal points
690  // X
691  if(shiftDir == 0 && shiftStep == 1){
692  for(int x = 0; x < lat.size[0]-1; x++){
693  for(int y = 0; y < lat.size[1]; y++){
694  for(int z = 0; z < lat.size[2]; z++){
695  for(int t = 0; t < lat.size[3]; t++){
696  shifted.at(x,y,z,t) = lat.at(x+1,y,z,t);
697  }}}}
698  }
699  else if(shiftDir == 0 && shiftStep == -1){
700  for(int x = 1; x < lat.size[0]; x++){
701  for(int y = 0; y < lat.size[1]; y++){
702  for(int z = 0; z < lat.size[2]; z++){
703  for(int t = 0; t < lat.size[3]; t++){
704  shifted.at(x,y,z,t) = lat.at(x-1,y,z,t);
705  }}}}
706  }
707  // Y
708  else if(shiftDir == 1 && shiftStep == 1){
709  for(int x = 0; x < lat.size[0]; x++){
710  for(int y = 0; y < lat.size[1]-1; y++){
711  for(int z = 0; z < lat.size[2]; z++){
712  for(int t = 0; t < lat.size[3]; t++){
713  shifted.at(x,y,z,t) = lat.at(x,y+1,z,t);
714  }}}}
715  }
716  else if(shiftDir == 1 && shiftStep == -1){
717  for(int x = 0; x < lat.size[0]; x++){
718  for(int y = 1; y < lat.size[1]; y++){
719  for(int z = 0; z < lat.size[2]; z++){
720  for(int t = 0; t < lat.size[3]; t++){
721  shifted.at(x,y,z,t) = lat.at(x,y-1,z,t);
722  }}}}
723  }
724  // Z
725  else if(shiftDir == 2 && shiftStep == 1){
726  for(int x = 0; x < lat.size[0]; x++){
727  for(int y = 0; y < lat.size[1]; y++){
728  for(int z = 0; z < lat.size[2]-1; z++){
729  for(int t = 0; t < lat.size[3]; t++){
730  shifted.at(x,y,z,t) = lat.at(x,y,z+1,t);
731  }}}}
732  }
733  else if(shiftDir == 2 && shiftStep == -1){
734  for(int x = 0; x < lat.size[0]; x++){
735  for(int y = 0; y < lat.size[1]; y++){
736  for(int z = 1; z < lat.size[2]; z++){
737  for(int t = 0; t < lat.size[3]; t++){
738  shifted.at(x,y,z,t) = lat.at(x,y,z-1,t);
739  }}}}
740  }
741  // T
742  else if(shiftDir == 3 && shiftStep == 1){
743  for(int x = 0; x < lat.size[0]; x++){
744  for(int y = 0; y < lat.size[1]; y++){
745  for(int z = 0; z < lat.size[2]; z++){
746  for(int t = 0; t < lat.size[3]-1; t++){
747  shifted.at(x,y,z,t) = lat.at(x,y,z,t+1);
748  }}}}
749  }
750  else if(shiftDir == 3 && shiftStep == -1){
751  for(int x = 0; x < lat.size[0]; x++){
752  for(int y = 0; y < lat.size[1]; y++){
753  for(int z = 0; z < lat.size[2]; z++){
754  for(int t = 1; t < lat.size[3]; t++){
755  shifted.at(x,y,z,t) = lat.at(x,y,z,t-1);
756  }}}}
757  }
758 
759 
760 
761  // Wait for the message to arrive
762  MPI_Wait(&sendRequest, &status);
763  MPI_Wait(&recvRequest, &status);
764 
765 
766  // Complete the result with the received data
767  // X
768  if(shiftDir == 0 && shiftStep == 1){
769  for(int y = 0; y < lat.size[1]; y++){
770  for(int z = 0; z < lat.size[2]; z++){
771  for(int t = 0; t < lat.size[3]; t++){
772  shifted.at(lat.size[0] - 1,y,z,t) = recvBuffer.at(0,y,z,t);
773  }}}
774  }
775  else if(shiftDir == 0 && shiftStep == -1){
776  for(int y = 0; y < lat.size[1]; y++){
777  for(int z = 0; z < lat.size[2]; z++){
778  for(int t = 0; t < lat.size[3]; t++){
779  shifted.at(0,y,z,t) = recvBuffer.at(0,y,z,t);
780  }}}
781  }
782 
783  // Y
784  else if(shiftDir == 1 && shiftStep == 1){
785  for(int x = 0; x < lat.size[0]; x++){
786  for(int z = 0; z < lat.size[2]; z++){
787  for(int t = 0; t < lat.size[3]; t++){
788  shifted.at(x,lat.size[1] - 1,z,t) = recvBuffer.at(x,0,z,t);
789  }}}
790  }
791  else if(shiftDir == 1 && shiftStep == -1){
792  for(int x = 0; x < lat.size[0]; x++){
793  for(int z = 0; z < lat.size[2]; z++){
794  for(int t = 0; t < lat.size[3]; t++){
795  shifted.at(x,0,z,t) = recvBuffer.at(x,0,z,t);
796  }}}
797  }
798 
799  // Z
800  else if(shiftDir == 2 && shiftStep == 1){
801  for(int x = 0; x < lat.size[0]; x++){
802  for(int y = 0; y < lat.size[1]; y++){
803  for(int t = 0; t < lat.size[3]; t++){
804  shifted.at(x,y,lat.size[2] - 1,t) = recvBuffer.at(x,y,0,t);
805  }}}
806  }
807  else if(shiftDir == 2 && shiftStep == -1){
808  for(int x = 0; x < lat.size[0]; x++){
809  for(int y = 0; y < lat.size[1]; y++){
810  for(int t = 0; t < lat.size[3]; t++){
811  shifted.at(x,y,0,t) = recvBuffer.at(x,y,0,t);
812  }}}
813  }
814  else if(shiftDir == 3 && shiftStep == 1){
815  for(int x = 0; x < lat.size[0]; x++){
816  for(int y = 0; y < lat.size[1]; y++){
817  for(int z = 0; z < lat.size[2]; z++){
818  shifted.at(x,y,z,lat.size[3] - 1) = recvBuffer.at(x,y,z,0);
819  }}}
820  }
821  else if(shiftDir == 3 && shiftStep == -1){
822  for(int x = 0; x < lat.size[0]; x++){
823  for(int y = 0; y < lat.size[1]; y++){
824  for(int z = 0; z < lat.size[2]; z++){
825  shifted.at(x,y,z,0) = recvBuffer.at(x,y,z,0);
826  }}}
827  }
828 
829  // Return
830  return shifted;
831 }
SU3
Implementation of a class to perform arithmetics between links.
Definition: su3.h:53
setLatticeImagIdentityValue
void setLatticeImagIdentityValue(Lattice< SU3 > &su3lat, const Lattice< double > &lat)
sets a lattice object to only diagonal matrices with complex values
Definition: lattice.cpp:50
setToIdentity
void setToIdentity(Lattice< SU3 > &su3lat)
sets a lattice object to all identity matrices
Definition: lattice.cpp:79
setToZero
void setToZero(Lattice< SU3 > &su3lat)
sets a lattice object to all zero matrices
Definition: lattice.cpp:71
parallel.h
Utilities for parallelization.
su3.h
Basic library to implement SU3 matrix arithmetics and functions.
exp
SU3 exp(const SU3 &Q)
exponential of a SU3 algebra matrix
Definition: su3.cpp:136
Parallel::getNeighbor
static int getNeighbor(int direction, int sign)
return the rank of the neighbor along the given direction and sign
Definition: parallel.cpp:203
Lattice
Template class to store an array with 4 dimensional indices of a given datatype. Includes functionali...
Definition: action.h:38
Generated by   doxygen 1.8.11