d75b7586d296e0f041f778c7956bf011f2304a0e
[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           BgPrintf("BgPrint> Start of iteration at %f\n");
277       }
278       iterations++;
279
280       // Copy different faces into messages
281       ghostMsg *leftMsg = new (blockDimY*blockDimZ) ghostMsg(RIGHT, blockDimY, blockDimZ);
282       ghostMsg *rightMsg = new (blockDimY*blockDimZ) ghostMsg(LEFT, blockDimY, blockDimZ);
283       ghostMsg *topMsg = new (blockDimX*blockDimZ) ghostMsg(BOTTOM, blockDimX, blockDimZ);
284       ghostMsg *bottomMsg = new (blockDimX*blockDimZ) ghostMsg(TOP, blockDimX, blockDimZ);
285       ghostMsg *frontMsg = new (blockDimX*blockDimY) ghostMsg(BACK, blockDimX, blockDimY);
286       ghostMsg *backMsg = new (blockDimX*blockDimY) ghostMsg(FRONT, blockDimX, blockDimY);
287
288       CkSetRefNum(leftMsg, iterations);
289       CkSetRefNum(rightMsg, iterations);
290       CkSetRefNum(topMsg, iterations);
291       CkSetRefNum(bottomMsg, iterations);
292       CkSetRefNum(frontMsg, iterations);
293       CkSetRefNum(backMsg, iterations);
294
295       for(int j=0; j<blockDimY; ++j) 
296         for(int k=0; k<blockDimZ; ++k) {
297           leftMsg->gh[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
298           rightMsg->gh[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
299       }
300
301       for(int i=0; i<blockDimX; ++i) 
302         for(int k=0; k<blockDimZ; ++k) {
303           topMsg->gh[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
304           bottomMsg->gh[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
305       }
306
307       for(int i=0; i<blockDimX; ++i) 
308         for(int j=0; j<blockDimY; ++j) {
309           frontMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
310           backMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
311       }
312
313       // Send my left face
314       thisProxy(wrap_x(thisIndex.x-1), thisIndex.y, thisIndex.z).receiveGhosts(leftMsg);
315       // Send my right face
316       thisProxy(wrap_x(thisIndex.x+1), thisIndex.y, thisIndex.z).receiveGhosts(rightMsg);
317       // Send my top face
318       thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z).receiveGhosts(topMsg);
319       // Send my bottom face
320       thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z).receiveGhosts(bottomMsg);
321       // Send my front face
322       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z-1)).receiveGhosts(frontMsg);
323       // Send my back face
324       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z+1)).receiveGhosts(backMsg);
325     }
326
327     void processGhosts(ghostMsg *gmsg) {
328       int height = gmsg->height;
329       int width = gmsg->width;
330
331       switch(gmsg->dir) {
332         case LEFT:
333           for(int j=0; j<height; ++j) 
334             for(int k=0; k<width; ++k) {
335               temperature[index(0, j+1, k+1)] = gmsg->gh[k*height+j];
336           }
337           break;
338         case RIGHT:
339           for(int j=0; j<height; ++j) 
340             for(int k=0; k<width; ++k) {
341               temperature[index(blockDimX+1, j+1, k+1)] = gmsg->gh[k*height+j];
342           }
343           break;
344         case TOP:
345           for(int i=0; i<height; ++i) 
346             for(int k=0; k<width; ++k) {
347               temperature[index(i+1, 0, k+1)] = gmsg->gh[k*height+i];
348           }
349           break;
350         case BOTTOM:
351           for(int i=0; i<height; ++i) 
352             for(int k=0; k<width; ++k) {
353               temperature[index(i+1, blockDimY+1, k+1)] = gmsg->gh[k*height+i];
354           }
355           break;
356         case FRONT:
357           for(int i=0; i<height; ++i) 
358             for(int j=0; j<width; ++j) {
359               temperature[index(i+1, j+1, blockDimZ+1)] = gmsg->gh[j*height+i];
360           }
361           break;
362         case BACK:
363           for(int i=0; i<height; ++i) 
364             for(int j=0; j<width; ++j) {
365               temperature[index(i+1, j+1, 0)] = gmsg->gh[j*height+i];
366           }
367           break;
368         default:
369           CkAbort("ERROR\n");
370       }
371
372       delete gmsg;
373     }
374
375
376     void check_and_compute() {
377       compute_kernel();
378
379       // calculate error
380       // not being done right now since we are doing a fixed no. of iterations
381
382 #if USE_3D_ARRAYS
383       double ***tmp;
384 #else
385       double *tmp;
386 #endif
387       tmp = temperature;
388       temperature = new_temperature;
389       new_temperature = tmp;
390
391       constrainBC();
392
393       if (iterations <= WARM_ITER || iterations >= MAX_ITER)
394         contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
395       else
396         doStep();
397     }
398
399     // Check to see if we have received all neighbor values yet
400     // If all neighbor values have been received, we update our values and proceed
401     void compute_kernel() {
402 #pragma unroll    
403       for(int i=1; i<blockDimX+1; ++i) {
404         for(int j=1; j<blockDimY+1; ++j) {
405           for(int k=1; k<blockDimZ+1; ++k) {
406             // update my value based on the surrounding values
407             new_temperature[index(i, j, k)] = (temperature[index(i-1, j, k)] 
408                                             +  temperature[index(i+1, j, k)]
409                                             +  temperature[index(i, j-1, k)]
410                                             +  temperature[index(i, j+1, k)]
411                                             +  temperature[index(i, j, k-1)]
412                                             +  temperature[index(i, j, k+1)]
413                                             +  temperature[index(i, j, k)] ) * DIVIDEBY7;
414           }
415         }
416       }
417     }
418
419     // Enforce some boundary conditions
420     void constrainBC() {
421       // Heat left, top and front faces of each chare's block
422       for(int i=1; i<blockDimX+1; ++i)
423         for(int k=1; k<blockDimZ+1; ++k)
424           temperature[index(i, 1, k)] = 255.0;
425       for(int j=1; j<blockDimY+1; ++j)
426         for(int k=1; k<blockDimZ+1; ++k)
427           temperature[index(1, j, k)] = 255.0;
428       for(int i=1; i<blockDimX+1; ++i)
429         for(int j=1; j<blockDimY+1; ++j)
430           temperature[index(i, j, 1)] = 255.0;
431     }
432
433 };
434
435 /** \class JacobiMap
436  *
437  */
438
439 class JacobiMap : public CkArrayMap {
440   public:
441     int X, Y, Z;
442     int *mapping;
443
444     JacobiMap(int x, int y, int z) {
445       X = x; Y = y; Z = z;
446       mapping = new int[X*Y*Z];
447
448       // we are assuming that the no. of chares in each dimension is a 
449       // multiple of the torus dimension
450
451       TopoManager tmgr;
452       int dimNX, dimNY, dimNZ, dimNT;
453
454       dimNX = tmgr.getDimNX();
455       dimNY = tmgr.getDimNY();
456       dimNZ = tmgr.getDimNZ();
457       dimNT = tmgr.getDimNT();
458
459       // we are assuming that the no. of chares in each dimension is a 
460       // multiple of the torus dimension
461       int numCharesPerPe = X*Y*Z/CkNumPes();
462
463       int numCharesPerPeX = X / dimNX;
464       int numCharesPerPeY = Y / dimNY;
465       int numCharesPerPeZ = Z / dimNZ;
466       int pe = 0, pes = CkNumPes();
467
468 #if USE_BLOCK_RNDMAP
469       int used[pes];
470       for(int i=0; i<pes; i++)
471         used[i] = 0;
472 #endif
473
474       if(dimNT < 2) {   // one core per node
475         if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ); 
476         for(int i=0; i<dimNX; i++)
477           for(int j=0; j<dimNY; j++)
478             for(int k=0; k<dimNZ; k++)
479             {
480 #if USE_BLOCK_RNDMAP
481               pe = myrand(pes); 
482               while(used[pe]!=0) {
483                 pe = myrand(pes); 
484               }
485               used[pe] = 1;
486 #endif
487
488               for(int ci=i*numCharesPerPeX; ci<(i+1)*numCharesPerPeX; ci++)
489                 for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
490                   for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
491 #if USE_TOPOMAP
492                     mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k);
493 #elif USE_BLOCK_RNDMAP
494                     mapping[ci*Y*Z + cj*Z + ck] = pe;
495 #endif
496                   }
497             }
498       } else {          // multiple cores per node
499         // In this case, we split the chares in the X dimension among the
500         // cores on the same node. The strange thing I figured out is that
501         // doing this in the Z dimension is not as good.
502         numCharesPerPeX /= dimNT;
503         if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ);
504
505         for(int i=0; i<dimNX; i++)
506           for(int j=0; j<dimNY; j++)
507             for(int k=0; k<dimNZ; k++)
508               for(int l=0; l<dimNT; l++)
509                 for(int ci=(dimNT*i+l)*numCharesPerPeX; ci<(dimNT*i+l+1)*numCharesPerPeX; ci++)
510                   for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
511                     for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
512                       mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k, l);
513                     }
514       } // end of if
515
516       if(CkMyPe() == 0) CkPrintf("Map generated ... \n");
517     }
518
519     ~JacobiMap() { 
520       delete [] mapping;
521     }
522
523     int procNum(int, const CkArrayIndex &idx) {
524       int *index = (int *)idx.data();
525       return mapping[index[0]*Y*Z + index[1]*Z + index[2]]; 
526     }
527 };
528
529 #include "jacobi3d.def.h"