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