Jacobi example: Remove old CVS headers
[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 /*readonly*/ int globalBarrier;
45
46 static unsigned long next = 1;
47
48 int myrand(int numpes) {
49   next = next * 1103515245 + 12345;
50   return((unsigned)(next/65536) % numpes);
51 }
52
53 // We want to wrap entries around, and because mod operator % 
54 // sometimes misbehaves on negative values. -1 maps to the highest value.
55 #define wrap_x(a)       (((a)+num_chare_x)%num_chare_x)
56 #define wrap_y(a)       (((a)+num_chare_y)%num_chare_y)
57 #define wrap_z(a)       (((a)+num_chare_z)%num_chare_z)
58
59 #define USE_3D_ARRAYS           0
60 #if USE_3D_ARRAYS
61 #define index(a, b, c)  a][b][c 
62 #else
63 #define index(a, b, c)  ( (a)*(blockDimY+2)*(blockDimZ+2) + (b)*(blockDimZ+2) + (c) )
64 #endif
65
66 #define MAX_ITER                26
67 #define WARM_ITER               5
68 #define LEFT                    1
69 #define RIGHT                   2
70 #define TOP                     3
71 #define BOTTOM                  4
72 #define FRONT                   5
73 #define BACK                    6
74 #define DIVIDEBY7               0.14285714285714285714
75
76 double startTime;
77 double endTime;
78
79 /** \class ghostMsg
80  *
81  */
82 class ghostMsg: public CMessage_ghostMsg {
83   public:
84     int dir;
85     int height;
86     int width;
87     double* gh;
88
89     ghostMsg(int _d, int _h, int _w) : dir(_d), height(_h), width(_w) {
90     }
91 };
92
93 /** \class Main
94  *
95  */
96 class Main : public CBase_Main {
97   public:
98     CProxy_Jacobi array;
99     int iterations;
100
101     Main(CkArgMsg* m) {
102       if ( (m->argc != 3) && (m->argc != 7) ) {
103         CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
104         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]);
105         CkAbort("Abort");
106       }
107
108       // set iteration counter to zero
109       iterations = 0;
110
111       // store the main proxy
112       mainProxy = thisProxy;
113         
114       if(m->argc == 3) {
115         arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
116         blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]); 
117       }
118       else if (m->argc == 7) {
119         arrayDimX = atoi(m->argv[1]);
120         arrayDimY = atoi(m->argv[2]);
121         arrayDimZ = atoi(m->argv[3]);
122         blockDimX = atoi(m->argv[4]); 
123         blockDimY = atoi(m->argv[5]); 
124         blockDimZ = atoi(m->argv[6]);
125       }
126
127       if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
128         CkAbort("array_size_X % block_size_X != 0!");
129       if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
130         CkAbort("array_size_Y % block_size_Y != 0!");
131       if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
132         CkAbort("array_size_Z % block_size_Z != 0!");
133
134       num_chare_x = arrayDimX / blockDimX;
135       num_chare_y = arrayDimY / blockDimY;
136       num_chare_z = arrayDimZ / blockDimZ;
137
138       // print info
139       CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
140       CkPrintf("Running Jacobi on %d processors with (%d, %d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y, num_chare_z);
141       CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
142       CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
143
144       // Create new array of worker chares
145 #if USE_TOPOMAP
146       CProxy_JacobiMap map = CProxy_JacobiMap::ckNew(num_chare_x, num_chare_y, num_chare_z);
147       CkPrintf("Topology Mapping is being done ... \n");
148       CkArrayOptions opts(num_chare_x, num_chare_y, num_chare_z);
149       opts.setMap(map);
150       array = CProxy_Jacobi::ckNew(opts);
151 #else
152       array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y, num_chare_z);
153 #endif
154
155       TopoManager tmgr;
156       CkArray *jarr = array.ckLocalBranch();
157       int jmap[num_chare_x][num_chare_y][num_chare_z];
158
159       int hops=0, p;
160       for(int i=0; i<num_chare_x; i++)
161         for(int j=0; j<num_chare_y; j++)
162           for(int k=0; k<num_chare_z; k++) {
163             jmap[i][j][k] = jarr->procNum(CkArrayIndex3D(i, j, k));
164           }
165
166       for(int i=0; i<num_chare_x; i++)
167         for(int j=0; j<num_chare_y; j++)
168           for(int k=0; k<num_chare_z; k++) {
169             p = jmap[i][j][k];
170             hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i+1)][j][k]);
171             hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i-1)][j][k]);
172             hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j+1)][k]);
173             hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j-1)][k]);
174             hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k+1)]);
175             hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k-1)]);
176           }
177       CkPrintf("Total Hops: %d\n", hops);
178
179       //Start the computation
180       array.doStep();
181     }
182
183     // Each worker reports back to here when it completes an iteration
184     void report() {
185       iterations++;
186       if (iterations <= WARM_ITER) {
187         if (iterations == WARM_ITER)
188           startTime = CmiWallTimer();
189         array.doStep();
190       }
191       else {
192         CkPrintf("Completed %d iterations\n", MAX_ITER-1);
193         endTime = CmiWallTimer();
194         CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(MAX_ITER-1-WARM_ITER));
195         CkExit();
196       }
197     }
198 };
199
200 /** \class Jacobi
201  *
202  */
203
204 class Jacobi: public CBase_Jacobi {
205   Jacobi_SDAG_CODE
206
207   public:
208     int iterations;
209     int imsg;
210
211 #if USE_3D_ARRAYS
212     double ***temperature;
213     double ***new_temperature;
214 #else
215     double *temperature;
216     double *new_temperature;
217 #endif
218
219     // Constructor, initialize values
220     Jacobi() {
221       __sdag_init();
222
223       int i, j, k;
224       // allocate a three dimensional array
225 #if USE_3D_ARRAYS
226       temperature = new double**[blockDimX+2];
227       new_temperature = new double**[blockDimX+2];
228       for (i=0; i<blockDimX+2; i++) {
229         temperature[i] = new double*[blockDimY+2];
230         new_temperature[i] = new double*[blockDimY+2];
231         for(j=0; j<blockDimY+2; j++) {
232           temperature[i][j] = new double[blockDimZ+2];
233           new_temperature[i][j] = new double[blockDimZ+2];
234         }
235       }
236 #else
237       temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
238       new_temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
239 #endif
240
241       for(i=0; i<blockDimX+2; ++i) {
242         for(j=0; j<blockDimY+2; ++j) {
243           for(k=0; k<blockDimZ+2; ++k) {
244             temperature[index(i, j, k)] = 0.0;
245           }
246         } 
247       }
248
249       iterations = 0;
250       imsg = 0;
251       constrainBC();
252     }
253
254     Jacobi(CkMigrateMessage* m) {}
255
256     ~Jacobi() { 
257 #if USE_3D_ARRAYS
258       for (int i=0; i<blockDimX+2; i++) {
259         for(int j=0; j<blockDimY+2; j++) {
260           delete [] temperature[i][j];
261           delete [] new_temperature[i][j];
262         }
263         delete [] temperature[i];
264         delete [] new_temperature[i];
265       }
266       delete [] temperature; 
267       delete [] new_temperature; 
268 #else
269       delete [] temperature; 
270       delete [] new_temperature; 
271 #endif
272     }
273
274     // Send ghost faces to the six neighbors
275     void begin_iteration(void) {
276       if (thisIndex.x == 0 && thisIndex.y == 0 && thisIndex.z == 0) {
277           CkPrintf("Start of iteration %d\n", iterations);
278           BgPrintf("BgPrint> Start of iteration at %f\n");
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 top face
320       thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z).receiveGhosts(topMsg);
321       // Send my bottom face
322       thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z).receiveGhosts(bottomMsg);
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 TOP:
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 BOTTOM:
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, blockDimZ+1)] = 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, 0)] = 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"