bug fix for index calculation
[charm.git] / examples / bigsim / sdag / jacobi3d / jacobi3d.C
1 /*****************************************************************************
2  * $Source$
3  * $Author$
4  * $Date$
5  * $Revision$
6  *****************************************************************************/
7
8 /** \file jacobi3d.C
9  *  Author: Abhinav S Bhatele
10  *  Date Created: June 01st, 2009
11  *
12  *  This does a topological placement for a 3d jacobi.
13  *
14  *      
15  *            *****************
16  *         *               *  *
17  *   ^  *****************     *
18  *   |  *               *     *
19  *   |  *               *     *
20  *   |  *               *     *
21  *   Y  *               *     *
22  *   |  *               *     *
23  *   |  *               *     *
24  *   |  *               *  * 
25  *   ~  *****************    Z
26  *      <------ X ------> 
27  *
28  *   X: left, right --> wrap_x
29  *   Y: top, bottom --> wrap_y
30  *   Z: front, back --> wrap_z
31  */
32
33 #include "jacobi3d.decl.h"
34 #include "TopoManager.h"
35
36 // See README for documentation
37
38 /*readonly*/ CProxy_Main mainProxy;
39 /*readonly*/ int arrayDimX;
40 /*readonly*/ int arrayDimY;
41 /*readonly*/ int arrayDimZ;
42 /*readonly*/ int blockDimX;
43 /*readonly*/ int blockDimY;
44 /*readonly*/ int blockDimZ;
45
46 // specify the number of worker chares in each dimension
47 /*readonly*/ int num_chare_x;
48 /*readonly*/ int num_chare_y;
49 /*readonly*/ int num_chare_z;
50
51 /*readonly*/ int globalBarrier;
52
53 static unsigned long next = 1;
54
55 int myrand(int numpes) {
56   next = next * 1103515245 + 12345;
57   return((unsigned)(next/65536) % numpes);
58 }
59
60 // We want to wrap entries around, and because mod operator % 
61 // sometimes misbehaves on negative values. -1 maps to the highest value.
62 #define wrap_x(a)       (((a)+num_chare_x)%num_chare_x)
63 #define wrap_y(a)       (((a)+num_chare_y)%num_chare_y)
64 #define wrap_z(a)       (((a)+num_chare_z)%num_chare_z)
65
66 #define USE_3D_ARRAYS           0
67 #if USE_3D_ARRAYS
68 #define index(a, b, c)  a][b][c 
69 #else
70 #define index(a, b, c)  ( (a)*(blockDimY+2)*(blockDimZ+2) + (b)*(blockDimZ+2) + (c) )
71 #endif
72
73 #define MAX_ITER                26
74 #define WARM_ITER               5
75 #define LEFT                    1
76 #define RIGHT                   2
77 #define TOP                     3
78 #define BOTTOM                  4
79 #define FRONT                   5
80 #define BACK                    6
81 #define DIVIDEBY7               0.14285714285714285714
82
83 double startTime;
84 double endTime;
85
86 /** \class ghostMsg
87  *
88  */
89 class ghostMsg: public CMessage_ghostMsg {
90   public:
91     int dir;
92     int height;
93     int width;
94     double* gh;
95
96     ghostMsg(int _d, int _h, int _w) : dir(_d), height(_h), width(_w) {
97     }
98 };
99
100 /** \class Main
101  *
102  */
103 class Main : public CBase_Main {
104   public:
105     CProxy_Jacobi array;
106     int iterations;
107
108     Main(CkArgMsg* m) {
109       if ( (m->argc != 3) && (m->argc != 7) ) {
110         CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
111         CkPrintf("OR %s [array_size_X] [array_size_Y] [array_size_Z] [block_size_X] [block_size_Y] [block_size_Z]\n", m->argv[0]);
112         CkAbort("Abort");
113       }
114
115       // set iteration counter to zero
116       iterations = 0;
117
118       // store the main proxy
119       mainProxy = thisProxy;
120         
121       if(m->argc == 3) {
122         arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
123         blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]); 
124       }
125       else if (m->argc == 7) {
126         arrayDimX = atoi(m->argv[1]);
127         arrayDimY = atoi(m->argv[2]);
128         arrayDimZ = atoi(m->argv[3]);
129         blockDimX = atoi(m->argv[4]); 
130         blockDimY = atoi(m->argv[5]); 
131         blockDimZ = atoi(m->argv[6]);
132       }
133
134       if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
135         CkAbort("array_size_X % block_size_X != 0!");
136       if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
137         CkAbort("array_size_Y % block_size_Y != 0!");
138       if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
139         CkAbort("array_size_Z % block_size_Z != 0!");
140
141       num_chare_x = arrayDimX / blockDimX;
142       num_chare_y = arrayDimY / blockDimY;
143       num_chare_z = arrayDimZ / blockDimZ;
144
145       // print info
146       CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
147       CkPrintf("Running Jacobi on %d processors with (%d, %d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y, num_chare_z);
148       CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
149       CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
150
151       // Create new array of worker chares
152 #if USE_TOPOMAP
153       CProxy_JacobiMap map = CProxy_JacobiMap::ckNew(num_chare_x, num_chare_y, num_chare_z);
154       CkPrintf("Topology Mapping is being done ... \n");
155       CkArrayOptions opts(num_chare_x, num_chare_y, num_chare_z);
156       opts.setMap(map);
157       array = CProxy_Jacobi::ckNew(opts);
158 #else
159       array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y, num_chare_z);
160 #endif
161
162       TopoManager tmgr;
163       CkArray *jarr = array.ckLocalBranch();
164       int jmap[num_chare_x][num_chare_y][num_chare_z];
165
166       int hops=0, p;
167       for(int i=0; i<num_chare_x; i++)
168         for(int j=0; j<num_chare_y; j++)
169           for(int k=0; k<num_chare_z; k++) {
170             jmap[i][j][k] = jarr->procNum(CkArrayIndex3D(i, j, k));
171           }
172
173       for(int i=0; i<num_chare_x; i++)
174         for(int j=0; j<num_chare_y; j++)
175           for(int k=0; k<num_chare_z; k++) {
176             p = jmap[i][j][k];
177             hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i+1)][j][k]);
178             hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i-1)][j][k]);
179             hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j+1)][k]);
180             hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j-1)][k]);
181             hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k+1)]);
182             hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k-1)]);
183           }
184       CkPrintf("Total Hops: %d\n", hops);
185
186       //Start the computation
187       array.doStep();
188     }
189
190     // Each worker reports back to here when it completes an iteration
191     void report() {
192       iterations++;
193       if (iterations <= WARM_ITER) {
194         if (iterations == WARM_ITER)
195           startTime = CmiWallTimer();
196         array.doStep();
197       }
198       else {
199         CkPrintf("Completed %d iterations\n", MAX_ITER-1);
200         endTime = CmiWallTimer();
201         CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(MAX_ITER-1-WARM_ITER));
202         CkExit();
203       }
204     }
205 };
206
207 /** \class Jacobi
208  *
209  */
210
211 class Jacobi: public CBase_Jacobi {
212   Jacobi_SDAG_CODE
213
214   public:
215     int iterations;
216     int imsg;
217
218 #if USE_3D_ARRAYS
219     double ***temperature;
220     double ***new_temperature;
221 #else
222     double *temperature;
223     double *new_temperature;
224 #endif
225
226     // Constructor, initialize values
227     Jacobi() {
228       __sdag_init();
229
230       int i, j, k;
231       // allocate a three dimensional array
232 #if USE_3D_ARRAYS
233       temperature = new double**[blockDimX+2];
234       new_temperature = new double**[blockDimX+2];
235       for (i=0; i<blockDimX+2; i++) {
236         temperature[i] = new double*[blockDimY+2];
237         new_temperature[i] = new double*[blockDimY+2];
238         for(j=0; j<blockDimY+2; j++) {
239           temperature[i][j] = new double[blockDimZ+2];
240           new_temperature[i][j] = new double[blockDimZ+2];
241         }
242       }
243 #else
244       temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
245       new_temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
246 #endif
247
248       for(i=0; i<blockDimX+2; ++i) {
249         for(j=0; j<blockDimY+2; ++j) {
250           for(k=0; k<blockDimZ+2; ++k) {
251             temperature[index(i, j, k)] = 0.0;
252           }
253         } 
254       }
255
256       iterations = 0;
257       imsg = 0;
258       constrainBC();
259     }
260
261     Jacobi(CkMigrateMessage* m) {}
262
263     ~Jacobi() { 
264 #if USE_3D_ARRAYS
265       for (int i=0; i<blockDimX+2; i++) {
266         for(int j=0; j<blockDimY+2; j++) {
267           delete [] temperature[i][j];
268           delete [] new_temperature[i][j];
269         }
270         delete [] temperature[i];
271         delete [] new_temperature[i];
272       }
273       delete [] temperature; 
274       delete [] new_temperature; 
275 #else
276       delete [] temperature; 
277       delete [] new_temperature; 
278 #endif
279     }
280
281     // Send ghost faces to the six neighbors
282     void begin_iteration(void) {
283       if (thisIndex.x == 0 && thisIndex.y == 0 && thisIndex.z == 0) {
284           CkPrintf("Start of iteration %d\n", iterations);
285           BgPrintf("BgPrint> Start of iteration at %f\n");
286       }
287       iterations++;
288
289       // Copy different faces into messages
290       ghostMsg *leftMsg = new (blockDimY*blockDimZ) ghostMsg(RIGHT, blockDimY, blockDimZ);
291       ghostMsg *rightMsg = new (blockDimY*blockDimZ) ghostMsg(LEFT, blockDimY, blockDimZ);
292       ghostMsg *topMsg = new (blockDimX*blockDimZ) ghostMsg(BOTTOM, blockDimX, blockDimZ);
293       ghostMsg *bottomMsg = new (blockDimX*blockDimZ) ghostMsg(TOP, blockDimX, blockDimZ);
294       ghostMsg *frontMsg = new (blockDimX*blockDimY) ghostMsg(BACK, blockDimX, blockDimY);
295       ghostMsg *backMsg = new (blockDimX*blockDimY) ghostMsg(FRONT, blockDimX, blockDimY);
296
297       CkSetRefNum(leftMsg, iterations);
298       CkSetRefNum(rightMsg, iterations);
299       CkSetRefNum(topMsg, iterations);
300       CkSetRefNum(bottomMsg, iterations);
301       CkSetRefNum(frontMsg, iterations);
302       CkSetRefNum(backMsg, iterations);
303
304       for(int j=0; j<blockDimY; ++j) 
305         for(int k=0; k<blockDimZ; ++k) {
306           leftMsg->gh[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
307           rightMsg->gh[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
308       }
309
310       for(int i=0; i<blockDimX; ++i) 
311         for(int k=0; k<blockDimZ; ++k) {
312           topMsg->gh[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
313           bottomMsg->gh[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
314       }
315
316       for(int i=0; i<blockDimX; ++i) 
317         for(int j=0; j<blockDimY; ++j) {
318           frontMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
319           backMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
320       }
321
322       // Send my left face
323       thisProxy(wrap_x(thisIndex.x-1), thisIndex.y, thisIndex.z).receiveGhosts(leftMsg);
324       // Send my right face
325       thisProxy(wrap_x(thisIndex.x+1), thisIndex.y, thisIndex.z).receiveGhosts(rightMsg);
326       // Send my top face
327       thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z).receiveGhosts(topMsg);
328       // Send my bottom face
329       thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z).receiveGhosts(bottomMsg);
330       // Send my front face
331       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z-1)).receiveGhosts(frontMsg);
332       // Send my back face
333       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z+1)).receiveGhosts(backMsg);
334     }
335
336     void processGhosts(ghostMsg *gmsg) {
337       int height = gmsg->height;
338       int width = gmsg->width;
339
340       switch(gmsg->dir) {
341         case LEFT:
342           for(int j=0; j<height; ++j) 
343             for(int k=0; k<width; ++k) {
344               temperature[index(0, j+1, k+1)] = gmsg->gh[k*height+j];
345           }
346           break;
347         case RIGHT:
348           for(int j=0; j<height; ++j) 
349             for(int k=0; k<width; ++k) {
350               temperature[index(blockDimX+1, j+1, k+1)] = gmsg->gh[k*height+j];
351           }
352           break;
353         case TOP:
354           for(int i=0; i<height; ++i) 
355             for(int k=0; k<width; ++k) {
356               temperature[index(i+1, 0, k+1)] = gmsg->gh[k*height+i];
357           }
358           break;
359         case BOTTOM:
360           for(int i=0; i<height; ++i) 
361             for(int k=0; k<width; ++k) {
362               temperature[index(i+1, blockDimY+1, k+1)] = gmsg->gh[k*height+i];
363           }
364           break;
365         case FRONT:
366           for(int i=0; i<height; ++i) 
367             for(int j=0; j<width; ++j) {
368               temperature[index(i+1, j+1, blockDimZ+1)] = gmsg->gh[j*height+i];
369           }
370           break;
371         case BACK:
372           for(int i=0; i<height; ++i) 
373             for(int j=0; j<width; ++j) {
374               temperature[index(i+1, j+1, 0)] = gmsg->gh[j*height+i];
375           }
376           break;
377         default:
378           CkAbort("ERROR\n");
379       }
380
381       delete gmsg;
382     }
383
384
385     void check_and_compute() {
386       compute_kernel();
387
388       // calculate error
389       // not being done right now since we are doing a fixed no. of iterations
390
391 #if USE_3D_ARRAYS
392       double ***tmp;
393 #else
394       double *tmp;
395 #endif
396       tmp = temperature;
397       temperature = new_temperature;
398       new_temperature = tmp;
399
400       constrainBC();
401
402       if (iterations <= WARM_ITER || iterations >= MAX_ITER)
403         contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
404       else
405         doStep();
406     }
407
408     // Check to see if we have received all neighbor values yet
409     // If all neighbor values have been received, we update our values and proceed
410     void compute_kernel() {
411 #pragma unroll    
412       for(int i=1; i<blockDimX+1; ++i) {
413         for(int j=1; j<blockDimY+1; ++j) {
414           for(int k=1; k<blockDimZ+1; ++k) {
415             // update my value based on the surrounding values
416             new_temperature[index(i, j, k)] = (temperature[index(i-1, j, k)] 
417                                             +  temperature[index(i+1, j, k)]
418                                             +  temperature[index(i, j-1, k)]
419                                             +  temperature[index(i, j+1, k)]
420                                             +  temperature[index(i, j, k-1)]
421                                             +  temperature[index(i, j, k+1)]
422                                             +  temperature[index(i, j, k)] ) * DIVIDEBY7;
423           }
424         }
425       }
426     }
427
428     // Enforce some boundary conditions
429     void constrainBC() {
430       // Heat left, top and front faces of each chare's block
431       for(int i=1; i<blockDimX+1; ++i)
432         for(int k=1; k<blockDimZ+1; ++k)
433           temperature[index(i, 1, k)] = 255.0;
434       for(int j=1; j<blockDimY+1; ++j)
435         for(int k=1; k<blockDimZ+1; ++k)
436           temperature[index(1, j, k)] = 255.0;
437       for(int i=1; i<blockDimX+1; ++i)
438         for(int j=1; j<blockDimY+1; ++j)
439           temperature[index(i, j, 1)] = 255.0;
440     }
441
442 };
443
444 /** \class JacobiMap
445  *
446  */
447
448 class JacobiMap : public CkArrayMap {
449   public:
450     int X, Y, Z;
451     int *mapping;
452
453     JacobiMap(int x, int y, int z) {
454       X = x; Y = y; Z = z;
455       mapping = new int[X*Y*Z];
456
457       // we are assuming that the no. of chares in each dimension is a 
458       // multiple of the torus dimension
459
460       TopoManager tmgr;
461       int dimNX, dimNY, dimNZ, dimNT;
462
463       dimNX = tmgr.getDimNX();
464       dimNY = tmgr.getDimNY();
465       dimNZ = tmgr.getDimNZ();
466       dimNT = tmgr.getDimNT();
467
468       // we are assuming that the no. of chares in each dimension is a 
469       // multiple of the torus dimension
470       int numCharesPerPe = X*Y*Z/CkNumPes();
471
472       int numCharesPerPeX = X / dimNX;
473       int numCharesPerPeY = Y / dimNY;
474       int numCharesPerPeZ = Z / dimNZ;
475       int pe = 0, pes = CkNumPes();
476
477 #if USE_BLOCK_RNDMAP
478       int used[pes];
479       for(int i=0; i<pes; i++)
480         used[i] = 0;
481 #endif
482
483       if(dimNT < 2) {   // one core per node
484         if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ); 
485         for(int i=0; i<dimNX; i++)
486           for(int j=0; j<dimNY; j++)
487             for(int k=0; k<dimNZ; k++)
488             {
489 #if USE_BLOCK_RNDMAP
490               pe = myrand(pes); 
491               while(used[pe]!=0) {
492                 pe = myrand(pes); 
493               }
494               used[pe] = 1;
495 #endif
496
497               for(int ci=i*numCharesPerPeX; ci<(i+1)*numCharesPerPeX; ci++)
498                 for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
499                   for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
500 #if USE_TOPOMAP
501                     mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k);
502 #elif USE_BLOCK_RNDMAP
503                     mapping[ci*Y*Z + cj*Z + ck] = pe;
504 #endif
505                   }
506             }
507       } else {          // multiple cores per node
508         // In this case, we split the chares in the X dimension among the
509         // cores on the same node. The strange thing I figured out is that
510         // doing this in the Z dimension is not as good.
511         numCharesPerPeX /= dimNT;
512         if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ);
513
514         for(int i=0; i<dimNX; i++)
515           for(int j=0; j<dimNY; j++)
516             for(int k=0; k<dimNZ; k++)
517               for(int l=0; l<dimNT; l++)
518                 for(int ci=(dimNT*i+l)*numCharesPerPeX; ci<(dimNT*i+l+1)*numCharesPerPeX; ci++)
519                   for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
520                     for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
521                       mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k, l);
522                     }
523       } // end of if
524
525       if(CkMyPe() == 0) CkPrintf("Map generated ... \n");
526     }
527
528     ~JacobiMap() { 
529       delete [] mapping;
530     }
531
532     int procNum(int, const CkArrayIndex &idx) {
533       int *index = (int *)idx.data();
534       return mapping[index[0]*Y*Z + index[1]*Z + index[2]]; 
535     }
536 };
537
538 #include "jacobi3d.def.h"