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