CHROMA
readcern.cc
Go to the documentation of this file.
1 
2 /*! \file
3  * \brief Read a CERN gauge configuration
4  */
5 
6 #include "chromabase.h"
7 #include "meas/glue/mesplq.h"
8 #include <vector>
9 using namespace std;
10 
11 namespace Chroma
12 {
13 
14 
15  // Read gauge field in CERN format into "u".
16  // Details about CERN format:
17  // - little endian, ints are 4 bytes, doubles are 8 bytes
18  // - 4 ints NT,NX,NY,NZ
19  // - 1 double for average plaquette
20  // - links as SU3 matrices (row major, double precision complex)
21  // in the following order:
22  // - the 8 links in directions +0,-0,...,+3,-3 at the first
23  // odd point, the second odd point, and so on.
24  // - The order of the point (x0,x1,x2,x3) with
25  // Cartesian coordinates in the range
26  // 0<=x0<N0,...,0<=x3<N3 is determined by the index
27  // ix=x3+N3*x2+N2*N3*x1+N1*N2*N3*x0,
28  // where N0,N1,N2,N3 are the global lattice sizes
29 
30 
31 
32 
33 void localPokeCernLink(const multi1d<int>& coord, const double *buf,
34  LatticeColorMatrix& Umufield)
35 {
36  int localsiteindex=Layout::linearSiteIndex(coord);
37  ColorMatrix locColMat;
38  int count=0;
39  for (int row=0;row<QDP::Nc;row++)
40  for (int col=0;col<QDP::Nc;col++){
41  Complex sitecomp=cmplx(Real64(buf[count]),Real64(buf[count+1]));
42  pokeColor(locColMat,sitecomp,row,col);
43  count+=2;}
44  Umufield.elem(localsiteindex) = locColMat.elem();
45 }
46 
47 
48 
49 void pokeCernLinks(const double *rdbuf, int lcbz, int linkdbles, multi1d<int>& coord,
50  int NX, int NY, int NZ, int NT, multi1d<LatticeColorMatrix>& u,
51  StopWatch& cmtimer)
52 {
53  // reorder the sites in the buffer
54  vector<double> cmbuf(8*lcbz*linkdbles);
55  int count=0;
56 #if defined(ARCH_PARSCALAR) || defined(ARCH_PARSCALARVEC)
57  if (Layout::primaryNode()){
58 #endif
59  for (int dir=0;dir<8;dir+=2)
60  for (int lz=0;lz<lcbz;lz++){
61  int indstart=(8*lz+dir)*linkdbles;
62  for (int k=0;k<linkdbles;k++)
63  cmbuf[count++]=rdbuf[indstart+k];}
64  for (int dir=7;dir>0;dir-=2)
65  for (int lz=lcbz-1;lz>=0;lz--){
66  int indstart=(8*lz+dir)*linkdbles;
67  for (int k=0;k<linkdbles;k++)
68  cmbuf[count++]=rdbuf[indstart+k];}
69 #if defined(ARCH_PARSCALAR) || defined(ARCH_PARSCALARVEC)
70  }
71 #endif
72 
73  int zkeep,ykeep,xkeep,tkeep,xshift,yshift,zshift,shift;
74 
75 #if defined(ARCH_PARSCALAR) || defined(ARCH_PARSCALARVEC)
76 
77  // get nodes involved
78  int cnode=Layout::nodeNumber(coord);
79  tkeep=coord[3];
80  if (coord[3]==0) coord[3]=NT-1; else coord[3]--;
81  int tmnode=Layout::nodeNumber(coord);
82  coord[3]=tkeep;
83  xkeep=coord[0];
84  if (coord[0]==0) coord[0]=NX-1; else coord[0]--;
85  int xmnode=Layout::nodeNumber(coord);
86  coord[0]=xkeep;
87  ykeep=coord[1];
88  if (coord[1]==0) coord[1]=NY-1; else coord[1]--;
89  int ymnode=Layout::nodeNumber(coord);
90  coord[1]=ykeep;
91  zkeep=coord[2];
92  if (coord[2]==0) coord[2]=NZ-1; else coord[2]--;
93  int zmnode=Layout::nodeNumber(coord);
94  coord[2]=zkeep;
95 
96  cmtimer.start();
97  size_t dbsize=sizeof(double);
98  if ((cnode==tmnode)&&(cnode==xmnode)&&(cnode==ymnode)&&(cnode==zmnode)){
99  size_t bytes=8*lcbz*linkdbles*dbsize;
100  if (cnode!=0){
101  if (Layout::primaryNode())
102  QDPInternal::sendToWait((void *)&cmbuf[0],cnode,bytes);
103  if (Layout::nodeNumber() == cnode)
104  QDPInternal::recvFromWait((void *)&cmbuf[0],0,bytes);}}
105  else{
106  if (cnode!=0){
107  size_t bytes=(5*lcbz-((cnode!=zmnode)?1:0))*linkdbles*dbsize;
108  if (Layout::primaryNode())
109  QDPInternal::sendToWait((void *)&cmbuf[0],cnode,bytes);
110  if (Layout::nodeNumber() == cnode)
111  QDPInternal::recvFromWait((void *)&cmbuf[0],0,bytes);}
112  count=(5*lcbz-1)*linkdbles;
113  if ((zmnode!=cnode)&&(zmnode!=0)){
114  size_t bytes=linkdbles*dbsize;
115  if (Layout::primaryNode())
116  QDPInternal::sendToWait((void *)&cmbuf[count],zmnode,bytes);
117  if (Layout::nodeNumber() == zmnode)
118  QDPInternal::recvFromWait((void *)&cmbuf[count],0,bytes);}
119  count=5*lcbz*linkdbles;
120  size_t bytes=lcbz*linkdbles*dbsize;
121  if (ymnode!=0){
122  if (Layout::primaryNode())
123  QDPInternal::sendToWait((void *)&cmbuf[count],ymnode,bytes);
124  if (Layout::nodeNumber() == ymnode)
125  QDPInternal::recvFromWait((void *)&cmbuf[count],0,bytes);}
126  count=6*lcbz*linkdbles;
127  if (xmnode!=0){
128  if (Layout::primaryNode())
129  QDPInternal::sendToWait((void *)&cmbuf[count],xmnode,bytes);
130  if (Layout::nodeNumber() == xmnode)
131  QDPInternal::recvFromWait((void *)&cmbuf[count],0,bytes);}
132  count=7*lcbz*linkdbles;
133  if (tmnode!=0){
134  if (Layout::primaryNode())
135  QDPInternal::sendToWait((void *)&cmbuf[count],tmnode,bytes);
136  if (Layout::nodeNumber() == tmnode)
137  QDPInternal::recvFromWait((void *)&cmbuf[count],0,bytes);}}
138  cmtimer.stop();
139 
140  if (Layout::nodeNumber()==cnode){
141 #endif
142  zkeep=coord[2];
143  xshift=lcbz*linkdbles;
144  yshift=2*lcbz*linkdbles;
145  zshift=3*lcbz*linkdbles;
146  shift=0;
147  for (int lz=0;lz<lcbz;lz++,shift+=linkdbles){
148  localPokeCernLink(coord,&cmbuf[xshift+shift],u[0]);
149  localPokeCernLink(coord,&cmbuf[yshift+shift],u[1]);
150  localPokeCernLink(coord,&cmbuf[zshift+shift],u[2]);
151  localPokeCernLink(coord,&cmbuf[shift],u[3]);
152  coord[2]+=2;}
153  coord[2]=zkeep+2*(lcbz-1)-1;
154  shift=4*lcbz*linkdbles;
155  for (int lz=1;lz<lcbz;lz++,shift+=linkdbles){
156  localPokeCernLink(coord,&cmbuf[shift],u[2]);
157  coord[2]-=2;}
158  coord[2]=zkeep;
159 #if defined(ARCH_PARSCALAR) || defined(ARCH_PARSCALARVEC)
160  }
161 
162  if (Layout::nodeNumber()==zmnode){
163 #endif
164  zkeep=coord[2];
165  if (coord[2]==0) coord[2]=NZ-1; else coord[2]--;
166  shift=(5*lcbz-1)*linkdbles;
167  localPokeCernLink(coord,&cmbuf[shift],u[2]);
168  coord[2]=zkeep;
169 #if defined(ARCH_PARSCALAR) || defined(ARCH_PARSCALARVEC)
170  }
171 
172  if (Layout::nodeNumber()==ymnode){
173 #endif
174  zkeep=coord[2];
175  ykeep=coord[1];
176  coord[2]=zkeep+2*(lcbz-1);
177  if (coord[1]==0) coord[1]=NY-1; else coord[1]--;
178  shift=5*lcbz*linkdbles;
179  for (int lz=0;lz<lcbz;lz++,shift+=linkdbles){
180  localPokeCernLink(coord,&cmbuf[shift],u[1]);
181  coord[2]-=2;}
182  coord[1]=ykeep;
183  coord[2]=zkeep;
184 #if defined(ARCH_PARSCALAR) || defined(ARCH_PARSCALARVEC)
185  }
186 
187  if (Layout::nodeNumber()==xmnode){
188 #endif
189  zkeep=coord[2];
190  xkeep=coord[0];
191  coord[2]=zkeep+2*(lcbz-1);
192  if (coord[0]==0) coord[0]=NX-1; else coord[0]--;
193  shift=6*lcbz*linkdbles;
194  for (int lz=0;lz<lcbz;lz++,shift+=linkdbles){
195  localPokeCernLink(coord,&cmbuf[shift],u[0]);
196  coord[2]-=2;}
197  coord[0]=xkeep;
198  coord[2]=zkeep;
199 #if defined(ARCH_PARSCALAR) || defined(ARCH_PARSCALARVEC)
200  }
201 
202  if (Layout::nodeNumber()==tmnode){
203 #endif
204  zkeep=coord[2];
205  tkeep=coord[3];
206  coord[2]=zkeep+2*(lcbz-1);
207  if (coord[3]==0) coord[3]=NT-1; else coord[3]--;
208  shift=7*lcbz*linkdbles;
209  for (int lz=0;lz<lcbz;lz++,shift+=linkdbles){
210  localPokeCernLink(coord,&cmbuf[shift],u[3]);
211  coord[2]-=2;}
212  coord[3]=tkeep;
213  coord[2]=zkeep;
214 #if defined(ARCH_PARSCALAR) || defined(ARCH_PARSCALARVEC)
215  }
216 
217 #endif
218 }
219 
220 
221 
222 
223 void readCERN(multi1d<LatticeColorMatrix>& u, const std::string& cfg_file)
224 {
225  if ((sizeof(int)!=4)||(sizeof(double)!=8)){
226  QDPIO::cout << "CERN files contain 4-byte ints, 8-byte doubles"<<std::endl;
227  QDP_abort(1);}
228  if (QDP::Nc!=3){
229  QDPIO::cout << "readCERN only supports Nc=3"<<std::endl;
230  QDP_abort(1);}
231  if (QDP::Nd!=4){
232  QDPIO::cout << "readCERN only supported for 4 space-time dimensions"<<std::endl;
233  QDP_abort(1);}
234 
235  int NT = QDP::Layout::lattSize()[3];
236  int NZ = QDP::Layout::lattSize()[2];
237  int NY = QDP::Layout::lattSize()[1];
238  int NX = QDP::Layout::lattSize()[0];
239  if (NZ%2){
240  QDPIO::cout << "NZ must be even"<<endl;
241  QDP_abort(1); }
242 
243  QDPIO::cout <<std::endl<< "Beginning read of CERN gauge field ("<<NX<<" x "<<NY
244  <<" x "<<NZ<<" ) x "<<NT<<std::endl;
245  StopWatch rtimer; rtimer.start();
246  StopWatch iotimer;
247  StopWatch cmtimer;
248 
249  // CERN field always stored as little endian
250  iotimer.start();
251  BinaryFileReader fin(cfg_file);
252  int Ncern[4];
253  fin.readArrayLittleEndian((char*)&Ncern[0],sizeof(int),4); // assumes little endian, no checksums
254  iotimer.stop();
255 
256  if ((Ncern[0]!=NT)||(Ncern[1]!=NX)||(Ncern[2]!=NY)||(Ncern[3]!=NZ)){
257  QDPIO::cout << "readCERN: Lattice size mismatch " << std::endl;
258  QDPIO::cout << "read "<<Ncern[1]<<" "<<Ncern[2]<<" "<<Ncern[3]
259  <<" "<<Ncern[0]<<std::endl;
260  QDPIO::cout << "Chroma wants "<<NX<<" "<<NY<<" "<<NZ<<" "<<NT<<std::endl;
261  QDP_abort(1); }
262  double plaq;
263  iotimer.start();
264  fin.readArrayLittleEndian((char*)&plaq,sizeof(double),1);
265  iotimer.stop();
266  plaq/=3.0;
267 
268  multi1d<int> gridsize=Layout::subgridLattSize();
269  int LZ=gridsize[2];
270  size_t dbsize=sizeof(double);
271  size_t linkdbles=2*QDP::Nc*QDP::Nc;
272  size_t rdlinks=4*NZ; // This is 8 directions * (NZ/2) odd sites in the z-direction; NZ must be even
273  size_t rdsize=rdlinks*linkdbles;
274  vector<double> rdbuf(rdsize);
275 
276  multi1d<int> coord(QDP::Nd);
277  for (int it=0;it<NT;it++){
278  coord[3]=it;
279  for (int ix=0;ix<NX;ix++){
280  coord[0]=ix;
281  for (int iy=0;iy<NY;iy++){
282  coord[1]=iy;
283  // read entire z on primary
284  iotimer.start();
285  fin.readArrayPrimaryNodeLittleEndian((char*)&rdbuf[0],dbsize,rdsize);
286  iotimer.stop();
287 
288  int izstart=((it+ix+iy)%2)?0:1;
289  int zshift=0;
290  int zz=0;
291  if (LZ%2){ // careful if LZ is odd
292  zz=(izstart)?1:-1;
293  }
294  for (int iz=izstart;iz<NZ;iz+=LZ+zz){
295  int lcbz=(LZ-zz)/2;
296  coord[2]=iz;
297  pokeCernLinks(&rdbuf[zshift],lcbz,linkdbles,coord,NX,NY,NZ,NT,u,cmtimer);
298  zshift+=8*lcbz*linkdbles;
299  if (zz!=0) zz=-zz;
300  }}}}
301 
302  rtimer.stop();
303  QDPIO::cout << "readCERN: plaq read: " << plaq << std::endl;
304  Double w_plaq,s_plaq,t_plaq,link;
305  MesPlq(u, w_plaq, s_plaq, t_plaq,link);
306  QDPIO::cout << "readCERN: plaq recomputed: " << w_plaq << std::endl;
307  QDPIO::cout << "Read of CERN gauge field done in "<<rtimer.getTimeInSeconds()<<" seconds"<<std::endl;
308  QDPIO::cout << "Time of IO operations = "<<iotimer.getTimeInSeconds()<<" seconds"<<std::endl;
309  QDPIO::cout << "Time of communications = "<<cmtimer.getTimeInSeconds()<<" seconds"<<std::endl<<std::endl;
310 #if defined(ARCH_PARSCALAR) || defined(ARCH_PARSCALARVEC)
311  QMP_barrier();
312 #endif
313 }
314 
315 
316 } // end namespace Chroma
317 
chromabase.h
Primary include file for CHROMA library code.
it
int it
Definition: meslate.cc:33
coord
multi1d< int > coord(Nd)
iy
int iy
Definition: meslate.cc:27
lz
int lz
Definition: meslate.cc:40
iz
int iz
Definition: meslate.cc:32
plaq
Double plaq
Definition: meslate.cc:58
Nd
Nd
Definition: meslate.cc:74
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::u
static multi1d< LatticeColorMatrix > u
Definition: syssolver_linop_qop_mg_w.cc:39
Chroma::pokeCernLinks
void pokeCernLinks(const double *rdbuf, int lcbz, int linkdbles, multi1d< int > &coord, int NX, int NY, int NZ, int NT, multi1d< LatticeColorMatrix > &u, StopWatch &cmtimer)
Definition: readcern.cc:49
Chroma::readCERN
void readCERN(multi1d< LatticeColorMatrix > &, const std::string &)
Definition: readcern.cc:223
Chroma::localPokeCernLink
void localPokeCernLink(const multi1d< int > &coord, const double *buf, LatticeColorMatrix &Umufield)
Definition: readcern.cc:33
Chroma::k
int k
Definition: invbicg.cc:119
Octave::count
int count
Definition: octave.h:14
std
Definition: gtest.h:745
testing::internal::Double
FloatingPoint< double > Double
Definition: gtest.h:7351
testing::internal::string
::std::string string
Definition: gtest.h:1979
link
Double link
Definition: pade_trln_w.cc:146
t_plaq
Double t_plaq
Definition: pade_trln_w.cc:145
w_plaq
Double w_plaq
Definition: pade_trln_w.cc:143
s_plaq
Double s_plaq
Definition: pade_trln_w.cc:144
MesPlq
MesPlq(u, w_plaq, s_plaq, t_plaq, link)
Generated by doxygen 1.9.1