su2hmc
Loading...
Searching...
No Matches
Site Indexing
Collaboration diagram for Site Indexing:

Functions

int Addrc (unsigned int *iu, unsigned int *id)
 Loads the addresses required during the update.
int ia (int x, int y, int z, int t)
 Described as a 21st Century address calculator, it gets the memory address of an array entry.
int Check_addr (unsigned int *table, int lns, int lnt, int imin, int imax)
 Checks that the addresses are within bounds before an update.
int Index2lcoord (int index, int *coord)
 Converts the index of a point in memory to the equivalent point in the 4 dimensional array, where the time index is the last coordinate in the array.
int Index2gcoord (int index, int *coord)
 Converts the index of a point in memory to the equivalent point in the 4 dimensional array, where the time index is the last coordinate in the array.
int Coord2lindex (int ix, int iy, int iz, int it)
 Converts the coordinates of a local lattice point to its index in the computer memory.
int Coord2gindex (int ix, int iy, int iz, int it)
 Converts the coordinates of a global lattice point to its index in the computer memory.
int Testlcoord (int cap)
 Tests if the local coordinate transformation functions are working.
int Testgcoord (int cap)
 This is completely new and missing from the original code.

Detailed Description

Function Documentation

◆ Addrc()

int Addrc ( unsigned int * iu,
unsigned int * id )

Loads the addresses required during the update.

Parameters
[out]iuUpper halo indices
[out]idLower halo indices
See also
hu, hd, h1u, h1d, h2u, h2d, halosize
Returns
Zero on success, integer error code otherwise

In original AOS code halo was stored after the sublattice

Definition at line 13 of file coord.c.

13 {
14 const char funcname[] = "Addrc";
15 //Rather than having 8 ih variables I'm going to use a 2x4 array
16 //down is 0, up is 1
17 int ih[2][4] = {{-1,-1,-1,-1},{-1,-1,-1,-1}};
18 if(nproc>1){
19 hd = (unsigned int *)aligned_alloc(AVX,ndim*halo*sizeof(int));
20 hu = (unsigned int *)aligned_alloc(AVX,ndim*halo*sizeof(int));
21 }
22
23 //Do the lookups appropriate for over indexing into halos
24 //order is down, up for each x y z t
25 //
26 // Since the h2? terms are related to the h1? terms I've dropped them in the C Version
27 // saving about 4 billionths of a second in the process
28 //
29 //Need to watch these +/- 1 at the end. Is that a FORTRAN thing or a program thing?
30 //The only time I see h1d called that +1 term gets cancelled by a -1 so I'm going
31 //to omit it here at my own peril. (Turned out I was right)
33 h1d[0]=kvol;
34 //For SOA each direction gets its own halo. Instead we have a j*kvol+h1?[i] term in the halo exchange
35 //h1d[0]=0;
36 h1u[0]=h1d[0]+halox;
37 halosize[0]=halox;
38
39 h1d[1]=h1u[0]+halox; h1u[1]=h1d[1]+haloy;
40 halosize[1]=haloy;
41
42 h1d[2]=h1u[1]+haloy; h1u[2]=h1d[2]+haloz;
43 halosize[2]=haloz;
44
45 h1d[3]=h1u[2]+haloz; h1u[3]=h1d[3]+halot;
46 halosize[3]=halot;
47
48 //Time for the nitty-gritty
49 /*
50 * Variables are:
51 *
52 * h1d(mu) = starting point in tail of down halo in direction mu
53 * h2d(mu) = finishing point in tail of down halo in direction mu
54 *
55 * h1u(mu) = starting point in tail of up halo in direction mu
56 * h2u(mu) = finishing point in tail of up halo in direction mu
57 *
58 * hd(i,mu) = index in core of point that should be packed into the
59 * ith location of the down halo in direction mu
60 *
61 * hu(i,mu) = index in core of point that should be packed into the
62 * ith location of the up halo in direction mu
63 *
64 * Note that hd and hu should be used for PACKING before SENDING
65 *
66 * Unpacking would be done with a loop over ALL the core sites with
67 * reference to normal dn/up lookups, ie we DO NOT have a list of
68 * where in the halo the core point i should go
69 *
70 * Halo points are ordered "as they come" in the linear loop over
71 * core sites
72 */
73 int iaddr, ic;
74 //if using ic++ inside the loop instead
75 // ic=-1;
76#ifdef _DEBUG
77 printf("ksizex = %i, ksizet=%i\n", ksizex, ksizet);
78#endif
79 //The loop order here matters as it affects the value of ic corresponding to each entry
80 for(int jt=0;jt<ksizet;jt++)
81 for(int jz=0;jz<ksizez;jz++)
82 for(int jy=0;jy<ksizey;jy++)
83 for(int jx=0;jx<ksizex;jx++){
84 //First value of ic is zero as planned.
85 //ic++;
86 ic=((jt*ksizez+jz)*ksizey+jy)*ksizex+jx;
87 //jx!=0 is logically equivalent to if(jx)
88 //If we're inside the sublattice, take the down nearest neightbour from inside the sublattice
89 if(jx)
90 iaddr = ia(jx-1,jy,jz,jt);
91 //Else if we're at the "down" edge, the down nearest neighbour is in the halo
92 else{
93 ih[0][0]++;
94#if npx>1
95 if(ih[0][0]>= halo){
96 fprintf(stderr, "Error %i in %s: Index ih[%i][%i]=%i is larger than the halo size %i."\
97 "\nExiting...\n\n", HALOLIM, funcname, 0, 0, ih[0][0], halo);
98#if(nproc>1)
99 MPI_Abort(comm,HALOLIM);
100#else
101 exit(HALOLIM);
102#endif
103 }
104 hd[0+ndim*ih[0][0]]=ic;
105 iaddr=h1d[0]+ih[0][0];
106#elif npx==1
107 iaddr = ia(jx-1,jy,jz,jt);
108#endif
109 }
110 id[0*kvol+ic]=iaddr;
111
112 if(jx<ksize-1)
113 iaddr = ia(jx+1,jy,jz,jt);
114 else{
115 ih[1][0]++;
116#if npx>1
117 if(ih[1][0]>= halo){
118 fprintf(stderr, "Error %i in %s: Index ih[%i][%i]=%i is larger than the halo size %i."
119 "\nExiting...\n\n", HALOLIM, funcname, 1, 0, ih[1][0], halo);
120#if(nproc>1)
121 MPI_Abort(comm,HALOLIM);
122#else
123 exit(HALOLIM);
124#endif
125 }
126 hu[0+ndim*ih[1][0]]=ic;
127 iaddr=ih[1][0]+h1u[0];
128#elif npx==1
129 iaddr = ia(jx+1,jy,jz,jt);
130#endif
131 }
132 iu[0*kvol+ic]=iaddr;
133
134 if(jy)
135 iaddr = ia(jx,jy-1,jz,jt);
136 else{
137 ih[0][1]++;
138#if npy>1
139 if(ih[0][1]>= halo){
140 fprintf(stderr, "Error %i in %s: Index ih[%i][%i]=%i is larger than the halo size %i."\
141 "\nExiting...\n\n", HALOLIM, funcname, 0, 1, ih[0][1], halo);
142#if(nproc>1)
143 MPI_Abort(comm,HALOLIM);
144#else
145 exit(HALOLIM);
146#endif
147 }
148 hd[1+ndim*ih[0][1]]=ic;
149 iaddr=h1d[1]+ih[0][1];
150#elif npy==1
151 iaddr = ia(jx,jy-1,jz,jt);
152#endif
153 }
154 id[1*kvol+ic]=iaddr;
155
156 if(jy<ksize-1)
157 iaddr = ia(jx,jy+1,jz,jt);
158 else{
159 ih[1][1]++;
160#if npy>1
161 if(ih[1][1]>= halo){
162 fprintf(stderr, "Error %i in %s: Index ih[%i][%i]=%i is larger than the halo size %i."
163 "\nExiting...\n\n", HALOLIM, funcname, 1, 1, ih[1][1], halo);
164#if(nproc>1)
165 MPI_Abort(comm,HALOLIM);
166#else
167 exit(HALOLIM);
168#endif
169 }
170 hu[1+ndim*ih[1][1]]=ic;
171 iaddr=ih[1][1]+h1u[1];
172#elif npy==1
173 iaddr = ia(jx,jy+1,jz,jt);
174#endif
175 }
176 iu[1*kvol+ic]=iaddr;
177
178 if(jz)
179 iaddr = ia(jx,jy,jz-1,jt);
180 else{
181 ih[0][2]++;
182#if npz>1
183 if(ih[0][2]>= halo){
184 fprintf(stderr, "Error %i in %s: Index ih[%i][%i]=%i is larger than the halo size %i."\
185 "\nExiting...\n\n", HALOLIM, funcname, 0, 2, ih[0][2], halo);
186#if(nproc>1)
187 MPI_Abort(comm,HALOLIM);
188#else
189 exit(HALOLIM);
190#endif
191 }
192 hd[2+ndim*ih[0][2]]=ic;
193 iaddr=h1d[2]+ih[0][2];
194#elif npz==1
195 iaddr = ia(jx,jy,jz-1,jt);
196#endif
197 }
198 id[2*kvol+ic]=iaddr;
199
200 if(jz<ksize-1)
201 iaddr = ia(jx,jy,jz+1,jt);
202 else{
203 ih[1][2]++;
204#if npz>1
205 if(ih[1][2]>= halo){
206 fprintf(stderr, "Error %i in %s: Index ih[%i][%i]=%i is larger than the halo size %i."
207 "\nExiting...\n\n", HALOLIM, funcname, 1, 2, ih[1][2], halo);
208#if(nproc>1)
209 MPI_Abort(comm,HALOLIM);
210#else
211 exit(HALOLIM);
212#endif
213 }
214 hu[2+ndim*ih[1][2]]=ic;
215 iaddr=ih[1][2]+h1u[2];
216#elif npz==1
217 iaddr = ia(jx,jy,jz+1,jt);
218#endif
219 }
220 iu[2*kvol+ic]=iaddr;
221
222 if(jt)
223 iaddr = ia(jx,jy,jz,jt-1);
224 else{
225 ih[0][3]++;
226#if npt>1
227 if(ih[0][3]>= halo){
228 fprintf(stderr, "Error %i in %s: Index ih[%i][%i]=%i is larger than the halo size %i."\
229 "\nExiting...\n\n", HALOLIM, funcname, 0, 3, ih[0][3], halo);
230#if(nproc>1)
231 MPI_Abort(comm,HALOLIM);
232#else
233 exit(HALOLIM);
234#endif
235 }
236 hd[3+ndim*ih[0][3]]=ic;
237 iaddr=h1d[3]+ih[0][3];
238#elif npt==1
239 iaddr = ia(jx,jy,jz,jt-1);
240#endif
241 }
242 id[3*kvol+ic]=iaddr;
243
244 if(jt<ksizet-1)
245 iaddr = ia(jx,jy,jz,jt+1);
246 else{
247 ih[1][3]++;
248#if npt>1
249 if(ih[1][3]>= halo){
250 fprintf(stderr, "Error %i in %s: Index ih[%i][%i]=%i is larger than the halo size %i."
251 "\nExiting...\n\n", HALOLIM, funcname, 1, 3, ih[1][3], halo);
252#if(nproc>1)
253 MPI_Abort(comm,HALOLIM);
254#else
255 exit(HALOLIM);
256#endif
257 }
258 hu[3+ndim*ih[1][3]]=ic;
259 iaddr=ih[1][3]+h1u[3];
260#elif npt==1
261 iaddr = ia(jx,jy,jz,jt+1);
262#endif
263 }
264 iu[3*kvol+ic]=iaddr;
265 }
266 //Print iu and id for diagnostics
267#ifdef _DEBUG
268#pragma omp parallel sections
269 {
270#pragma omp section
271 {
272 FILE *id_out = fopen("id_out", "w");
273 fprintf(id_out,"x\ty\tz\ti\n");
274 for(int i=0;i<kvol;i++)
275 fprintf(id_out,"%03i\t%03i\t%03i\t%03i\n",id[i],id[i+kvol*1],id[i+kvol*2],id[i+kvol*3]);
276 fclose(id_out);
277 }
278#pragma omp section
279 {
280 FILE *iu_out = fopen("iu_out", "w");
281 fprintf(iu_out,"x\ty\tz\ti\n");
282 for(int i=0;i<kvol;i++)
283 fprintf(iu_out,"%03i\t%03i\t%03i\t%03i\n",iu[i],iu[i+kvol*1],iu[i+kvol*2],iu[i+kvol*3]);
284 fclose(iu_out);
285
286 }
287
288 }
289#endif
290 return 0;
291}
unsigned int halosize[ndim]
Array containing the size of the halo in each direction.
Definition coord.c:12
unsigned int * hd
Down halo indices.
Definition coord.c:11
unsigned int * hu
Up halo indices.
Definition coord.c:11
unsigned int h1d[ndim]
Down halo starting element.
Definition coord.c:12
unsigned int h1u[ndim]
Up halo starting element.
Definition coord.c:12
#define HALOLIM
Index goes beyond the halo.
Definition errorcodes.h:117
int ia(int x, int y, int z, int t)
Described as a 21st Century address calculator, it gets the memory address of an array entry.
Definition coord.c:293
#define ksizex
Sublattice x extent.
Definition sizes.h:148
#define AVX
Alignment of arrays. 64 for AVX-512, 32 for AVX/AVX2. 16 for SSE. Since AVX is standard on modern x86...
Definition sizes.h:279
#define nproc
Number of processors for MPI.
Definition sizes.h:138
#define halot
t Halo size
Definition sizes.h:228
#define ksizet
Sublattice t extent.
Definition sizes.h:158
#define halox
x Halo size
Definition sizes.h:210
#define kvol
Sublattice volume.
Definition sizes.h:163
#define haloz
z Halo size
Definition sizes.h:222
#define haloy
y Halo size
Definition sizes.h:216
#define ksize
Sublattice spatial extent for a cubic lattice.
Definition sizes.h:155
#define halo
Total Halo size.
Definition sizes.h:231
#define ksizez
Sublattice z extent.
Definition sizes.h:152
#define ndim
Dimensions.
Definition sizes.h:188
#define ksizey
Sublattice y extent.
Definition sizes.h:150

References AVX, h1d, h1u, halo, HALOLIM, halosize, halot, halox, haloy, haloz, hd, hu, ia(), ksize, ksizet, ksizex, ksizey, ksizez, kvol, ndim, and nproc.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Check_addr()

int Check_addr ( unsigned int * table,
int lns,
int lnt,
int imin,
int imax )

Checks that the addresses are within bounds before an update.

Parameters
[in]tablePointer to the table in question
[in]lnsSize of each spacial dimension
[in]lntSize of the time dimension
[in]iminLower bound for element of the table
[in]imaxUpper bound for an element of the table
Returns
Zero on success, integer error code otherwise.

Definition at line 305 of file coord.c.

305 {
306 const char funcname[] = "Check_addr";
307 //Get the total number of elements in each dimension of the table
308 int ntable = lns*lns*lns*lnt;
309 int iaddr;
310 //Collapsing two for loops together
311 for(int j=0; j<ntable*ndim; j++){
312 iaddr = table[j];
313 if((iaddr<imin) || (iaddr>= imax)){
314 fprintf(stderr, "Error %i in %s: %i is out of the bounds of (%i,%i)\n"\
315 "for a table of size %i^3 *%i.\nExiting...\n\n",\
316 BOUNDERROR,funcname,iaddr,imin,imax,lns,lnt);
317#if(nproc>1)
318 MPI_Abort(comm,BOUNDERROR);
319#else
320 exit(BOUNDERROR);
321#endif
322 }
323 }
324 return 0;
325}
#define BOUNDERROR
Accessing out of bounds element.
Definition errorcodes.h:62

References BOUNDERROR, and ndim.

Here is the caller graph for this function:

◆ Coord2gindex()

int Coord2gindex ( int ix,
int iy,
int iz,
int it )
inline

Converts the coordinates of a global lattice point to its index in the computer memory.

This is a rather nuanced function, as C and Fortran are rather different in how they store arrays. C starts with index 0 and Fortran (by default) starts with index 1

Also C and Fortran store data in the opposite memory order so be careful when calling this function!

Parameters
[in]ix,iy,iz,itIndex in each direction
Returns
The position of the point in the flattened array

Definition at line 365 of file coord.c.

365 {
366 const char funcname[] = "Coord2gindex";
367
368 //I've factorised this function compared to its original
369 //implementation to reduce the number of multiplications
370 //and hopefully improve performance
371 // return it+nt*(iz+nz*(iy+ny*ix));
372 return ix+nx*(iy+ny*(iz+nz*it));
373}
#define nx
Lattice x extent.
Definition sizes.h:72
#define nz
Lattice z extent. We normally use cubic lattices so this is the same as nx.
Definition sizes.h:86
#define ny
Lattice y extent. We normally use cubic lattices so this is the same as nx.
Definition sizes.h:80

References nx, ny, and nz.

Here is the caller graph for this function:

◆ Coord2lindex()

int Coord2lindex ( int ix,
int iy,
int iz,
int it )
inline

Converts the coordinates of a local lattice point to its index in the computer memory.

This is a rather nuanced function, as C and Fortran are rather different in how they store arrays. C starts with index 0 and Fortran (by default) starts with index 1

Also C and Fortran store data in the opposite memory order so be careful when calling this function!

Parameters
[in]ix,iy,iz,itIndex in each direction
Returns
The position of the point in the flattened array

Definition at line 356 of file coord.c.

356 {
357 const char funcname[] = "Coord2gindex";
358
359 //I've factorised this function compared to its original
360 //implementation to reduce the number of multiplications
361 //and hopefully improve performance
362 //int index = coord[3]+ksizez*(coord[2]+ksizey*(coord[1]+ksizex*coord[0]));
363 return it+ksizet*(iz+ksizez*(iy+ksizey*ix));
364}

References ksizet, ksizey, and ksizez.

◆ ia()

int ia ( int x,
int y,
int z,
int t )
inline

Described as a 21st Century address calculator, it gets the memory address of an array entry.

Parameters
[in]x,y,z,tThe coordinates
Returns
An integer corresponding to the position of the entry in a flattened row-major array.
Todo
Future... Switch for Row and column major, and zero or one indexing

Definition at line 293 of file coord.c.

293 {
294 const char funcname[] = "ia";
295 //We need to ensure that the indices aren't out of bounds using while loops
296 while(x<0) x+=ksizex; while(x>=ksizex) x-= ksizex;
297 while(y<0) y+=ksizey; while(y>=ksizey) y-= ksizey;
298 while(z<0) z+=ksizez; while(z>=ksizez) z-= ksizez;
299 while(t<0) t+=ksizet; while(t>=ksizet) t-= ksizet;
300
301 //And flattening.
302 //return t+ksizet*(z+ksizez*(y+ksizey*x));
303 return ((t*ksizez+z)*ksizey+y)*ksizex+x;
304}

References ksizet, ksizex, ksizey, and ksizez.

Here is the caller graph for this function:

◆ Index2gcoord()

int Index2gcoord ( int index,
int * coord )
inline

Converts the index of a point in memory to the equivalent point in the 4 dimensional array, where the time index is the last coordinate in the array.

This is a rather nuanced function, as C and Fortran are rather different in how they store arrays. C starts with index 0 and Fortran (by default) starts with index 1

Also C and Fortran store data in the opposite memory order so be careful when calling this function!

Parameters
[in]indexThe index of the point as stored linearly in computer memory
[out]coordThe 4-array for the coordinates. The first three spots are for the time index.
Returns
Zero on success. Integer Error code otherwise

Definition at line 341 of file coord.c.

341 {
342
343 const char funcname[] = "Index2gcoord";
344 //A divide and conquer approach. Going from the deepest coordinate
345 //to the least deep coordinate, we take the modulo of the index by
346 //the length of that axis to get the coordinate, and then divide
347 //the index by the length of that coordinate to set up for the
348 //next coordinate. This works since int/int gives an int.
349 coord[3] = index%nt; index/=nt;
350 coord[2] = index%nz; index/=nz;
351 coord[1] = index%ny; index/=ny;
352 coord[0] = index; //No need to divide by nt since were done.
353
354 return 0;
355}
#define nt
Lattice temporal extent. This also corresponds to the inverse temperature.
Definition sizes.h:92

References nt, ny, and nz.

Here is the caller graph for this function:

◆ Index2lcoord()

int Index2lcoord ( int index,
int * coord )
inline

Converts the index of a point in memory to the equivalent point in the 4 dimensional array, where the time index is the last coordinate in the array.

This is a rather nuanced function, as C and Fortran are rather different in how they store arrays. C starts with index 0 and Fortran (by default) starts with index 1

Also C and Fortran store data in the opposite memory order so be careful when calling this function!

Parameters
[in]indexThe index of the point as stored linearly in computer memory
[out]coordThe 4-array for the coordinates. The first three spots are for the time index.
Returns
Zero on success. Integer Error code otherwise

Definition at line 326 of file coord.c.

326 {
327
328 const char funcname[] = "Index2lcoord";
329 //A divide and conquer approach. Going from the deepest coordinate
330 //to the least deep coordinate, we take the modulo of the index by
331 //the length of that axis to get the coordinate, and then divide
332 //the index by the length of that coordinate to set up for the
333 //next coordinate. This works since int/int gives an int.
334 coord[3] = index%ksizet; index/=ksizet;
335 coord[2] = index%ksizez; index/=ksizez;
336 coord[1] = index%ksizey; index/=ksizey;
337 coord[0] = index; //No need to divide by nt since were done.
338
339 return 0;
340}

References ksizet, ksizey, and ksizez.

Here is the caller graph for this function:

◆ Testgcoord()

int Testgcoord ( int cap)

This is completely new and missing from the original code.

We test the coordinate conversion functions by doing the following

  1. Convert from int to gcoord (new)
  2. Convert from gcoord to int (also new) and compare to input. If we get the same value we started with then we're probably doing something right

The code is basically the same as the previous function with different magic numbers.

Parameters
[in]capThe max value the index can take on. Should be the size of our array
Returns
Zero on success, integer error code otherwise

Definition at line 396 of file coord.c.

396 {
397 const char funcname[] = "Testgcoord";
398 int coord[4], index, index2;
399#pragma omp parallel for private(coord, index, index2)
400 for(index=0; index<cap; index++){
401 Index2gcoord(index, coord);
402#pragma omp critical
403 printf("Coordinates for %i are (x,y,z,t):[%i,%i,%i,%i].\n", index,\
404 coord[0], coord[1], coord[2], coord[3]);
405 //index2 = Coord2gindex(coord);
406 if(!(index==index2)){
407 fprintf(stderr, "Error %i in %s: Converted index %i does not match "\
408 "original index %i.\nExiting...\n\n",\
409 INDTOCOORD, funcname, index2, index);
410#if(nproc>1)
411 MPI_Abort(comm,INDTOCOORD);
412#else
413 exit(INDTOCOORD);
414#endif
415 }
416 }
417 return 0;
418}
#define INDTOCOORD
Issues converting index to coordinates.
Definition errorcodes.h:58
int Index2gcoord(int index, int *coord)
Converts the index of a point in memory to the equivalent point in the 4 dimensional array,...
Definition coord.c:341

References Index2gcoord(), and INDTOCOORD.

Here is the call graph for this function:

◆ Testlcoord()

int Testlcoord ( int cap)

Tests if the local coordinate transformation functions are working.

Going to expand a little on the original here and do the following

  1. Convert from int to lcoord (the original code) And the planned additional features
  2. Convert from lcoord to int (new, function doesn't exist in the original If we get the same value we started with then we're probably doing something right.
Parameters
[in]capThe max value the index can take on. Should be the size of the array
Returns
Zero on success, integer error code otherwise.

Definition at line 374 of file coord.c.

374 {
375 const char funcname[] = "Testlcoord";
376 //The storage array for the coordinates, and the index and its test value.
377 int coord[4], index, index2;
378 for(index =0; index<cap; index++){
379 Index2lcoord(index, coord);
380 printf("Coordinates for %i are (x,y,z,t):[%i,%i,%i,%i].\n", index,\
381 coord[0], coord[1], coord[2], coord[3]);
382 //index2 = Coord2lindex(coord);
383 if(!(index==index2)){
384 fprintf(stderr, "Error %i in %s: Converted index %i does not match "
385 "original index %i.\nExiting...\n\n",\
386 INDTOCOORD, funcname, index2, index);
387#if(nproc>1)
388 MPI_Abort(comm,INDTOCOORD);
389#else
390 exit(INDTOCOORD);
391#endif
392 }
393 }
394 return 0;
395}
int Index2lcoord(int index, int *coord)
Converts the index of a point in memory to the equivalent point in the 4 dimensional array,...
Definition coord.c:326

References Index2lcoord(), and INDTOCOORD.

Here is the call graph for this function: