add Jacobi3d and gauss-siedal combination
[charm.git] / tests / charm++ / jacobi3d-gausssiedal / 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 #ifdef JACOBI_OPENMP
36 #include <omp.h>
37 #endif
38
39 // See README for documentation
40
41 /*readonly*/ CProxy_Main mainProxy;
42 /*readonly*/ int arrayDimX;
43 /*readonly*/ int arrayDimY;
44 /*readonly*/ int arrayDimZ;
45 /*readonly*/ int blockDimX;
46 /*readonly*/ int blockDimY;
47 /*readonly*/ int blockDimZ;
48
49 // specify the number of worker chares in each dimension
50 /*readonly*/ int num_chare_x;
51 /*readonly*/ int num_chare_y;
52 /*readonly*/ int num_chare_z;
53
54 /*readonly*/ int globalBarrier;
55
56 static unsigned long next = 1;
57
58 int myrand(int numpes) {
59   next = next * 1103515245 + 12345;
60   return((unsigned)(next/65536) % numpes);
61 }
62
63 // We want to wrap entries around, and because mod operator % 
64 // sometimes misbehaves on negative values. -1 maps to the highest value.
65 #define wrap_x(a)       (((a)+num_chare_x)%num_chare_x)
66 #define wrap_y(a)       (((a)+num_chare_y)%num_chare_y)
67 #define wrap_z(a)       (((a)+num_chare_z)%num_chare_z)
68
69 #define index(a, b, c)  ( (a)*(blockDimY+2)*(blockDimZ+2) + (b)*(blockDimZ+2) + (c) )
70
71 #define PRINT_FREQ      100
72 #define CKP_FREQ                100
73 #define MAX_ITER         400    
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 #ifdef JACOBI_OPENMP 
187       CProxy_OmpInitializer ompInit = 
188         CProxy_OmpInitializer::ckNew(4);
189 #else    
190       //Start the computation
191       start();
192 #endif
193     }
194
195     void start() {
196       startTime = CmiWallTimer();                
197       array.doStep();
198     }
199
200     // Each worker reports back to here when it completes an iteration
201         void report() {
202         iterations += CKP_FREQ;
203         if (iterations < MAX_ITER) {
204 #ifdef CMK_MEM_CHECKPOINT
205             CkCallback cb (CkIndex_Jacobi::doStep(), array);
206             CkStartMemCheckpoint(cb);
207 #else
208                         array.doStep();
209 #endif
210         } else {
211             CkPrintf("Completed %d iterations\n", MAX_ITER-1);
212             endTime = CmiWallTimer();
213             CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(MAX_ITER-1-WARM_ITER));
214             CkExit();
215         }
216     }
217
218 };
219
220 /** \class Jacobi
221  *
222  */
223
224 class Jacobi: public CBase_Jacobi {
225   Jacobi_SDAG_CODE
226
227   public:
228     int iterations;
229     int imsg;
230
231     double *temperature;
232         double timing,average;
233
234     // Constructor, initialize values
235     Jacobi() {
236       // This call is an anachronism - the need to call __sdag_init() has been
237       // removed. We still call it here to test backward compatibility.
238       __sdag_init();
239
240       int i, j, k;
241       // allocate a three dimensional array
242       temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
243
244       for(i=0; i<blockDimX+2; ++i) {
245           for(j=0; j<blockDimY+2; ++j) {
246               for(k=0; k<blockDimZ+2; ++k) {
247                   temperature[index(i, j, k)] = 0.0;
248               }
249           } 
250       }
251
252       iterations = 0;
253       imsg = 0;
254       constrainBC();
255
256       usesAtSync = CmiTrue;
257
258     }
259
260     Jacobi(CkMigrateMessage* m): CBase_Jacobi(m) { }
261
262     ~Jacobi() { 
263         delete [] temperature; 
264     }
265
266         // Pupping function for migration and fault tolerance
267         // Condition: assuming the 3D Chare Arrays are NOT used
268         void pup(PUP::er &p){
269                 
270                 // calling parent's pup
271                 CBase_Jacobi::pup(p);
272         
273                 __sdag_pup(p);
274         
275                 // pupping properties of this class
276                 p | iterations;
277                 p | imsg;
278
279                 // if unpacking, allocate the memory space
280                 if(p.isUnpacking()){
281                         temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
282         
283                 }
284
285                 // pupping the arrays
286                 p((char *)temperature, (blockDimX+2) * (blockDimY+2) * (blockDimZ+2) * sizeof(double));
287
288         }
289
290
291     // Send ghost faces to the six neighbors
292     void begin_iteration(void) {
293       if (thisIndex.x == 0 && thisIndex.y == 0 && thisIndex.z == 0) {
294 //          CkPrintf("Start of iteration %d\n", iterations);
295                         if(iterations % PRINT_FREQ == 0){
296                                 average = timing;
297                                 timing = CmiWallTimer();
298                                 average = (timing - average)/(double)PRINT_FREQ;
299                                 CkPrintf("time=%.2f it=%d avg=%.4f\n",timing,iterations,average);
300                         }
301       }
302       iterations++;
303
304       // Copy different faces into messages
305       ghostMsg *leftMsg = new (blockDimY*blockDimZ) ghostMsg(RIGHT, blockDimY, blockDimZ);
306       ghostMsg *rightMsg = new (blockDimY*blockDimZ) ghostMsg(LEFT, blockDimY, blockDimZ);
307       ghostMsg *topMsg = new (blockDimX*blockDimZ) ghostMsg(BOTTOM, blockDimX, blockDimZ);
308       ghostMsg *bottomMsg = new (blockDimX*blockDimZ) ghostMsg(TOP, blockDimX, blockDimZ);
309       ghostMsg *frontMsg = new (blockDimX*blockDimY) ghostMsg(BACK, blockDimX, blockDimY);
310       ghostMsg *backMsg = new (blockDimX*blockDimY) ghostMsg(FRONT, blockDimX, blockDimY);
311
312       CkSetRefNum(leftMsg, iterations);
313       CkSetRefNum(rightMsg, iterations);
314       CkSetRefNum(topMsg, iterations);
315       CkSetRefNum(bottomMsg, iterations);
316       CkSetRefNum(frontMsg, iterations);
317       CkSetRefNum(backMsg, iterations);
318
319       for(int j=0; j<blockDimY; ++j) 
320           for(int k=0; k<blockDimZ; ++k) {
321               leftMsg->gh[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
322               rightMsg->gh[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
323           }
324
325       for(int i=0; i<blockDimX; ++i) 
326           for(int k=0; k<blockDimZ; ++k) {
327               topMsg->gh[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
328               bottomMsg->gh[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
329           }
330
331       for(int i=0; i<blockDimX; ++i) 
332           for(int j=0; j<blockDimY; ++j) {
333               frontMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
334               backMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
335           }
336
337       // Send my left face
338        thisProxy(wrap_x(thisIndex.x-1), thisIndex.y, thisIndex.z).receiveGhosts(leftMsg);
339       // Send my right face
340       thisProxy(wrap_x(thisIndex.x+1), thisIndex.y, thisIndex.z).receiveGhosts(rightMsg);
341       // Send my top face
342       thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z).receiveGhosts(topMsg);
343       // Send my bottom face
344       thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z).receiveGhosts(bottomMsg);
345       // Send my front face
346       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z-1)).receiveGhosts(frontMsg);
347       // Send my back face
348       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z+1)).receiveGhosts(backMsg);
349     }
350
351     void processGhosts(ghostMsg *gmsg) {
352       int height = gmsg->height;
353       int width = gmsg->width;
354
355       switch(gmsg->dir) {
356       case LEFT:
357           for(int j=0; j<height; ++j) 
358               for(int k=0; k<width; ++k) {
359                   temperature[index(0, j+1, k+1)] = gmsg->gh[k*height+j];
360               }
361           break;
362       case RIGHT:
363           for(int j=0; j<height; ++j) 
364               for(int k=0; k<width; ++k) {
365                   temperature[index(blockDimX+1, j+1, k+1)] = gmsg->gh[k*height+j];
366               }
367           break;
368       case TOP:
369           for(int i=0; i<height; ++i) 
370               for(int k=0; k<width; ++k) {
371                   temperature[index(i+1, 0, k+1)] = gmsg->gh[k*height+i];
372               }
373           break;
374       case BOTTOM:
375           for(int i=0; i<height; ++i) 
376               for(int k=0; k<width; ++k) {
377                   temperature[index(i+1, blockDimY+1, k+1)] = gmsg->gh[k*height+i];
378               }
379           break;
380       case FRONT:
381           for(int i=0; i<height; ++i) 
382               for(int j=0; j<width; ++j) {
383                   temperature[index(i+1, j+1, blockDimZ+1)] = gmsg->gh[j*height+i];
384               }
385           break;
386       case BACK:
387           for(int i=0; i<height; ++i) 
388               for(int j=0; j<width; ++j) {
389                   temperature[index(i+1, j+1, 0)] = gmsg->gh[j*height+i];
390               }
391           break;
392       default:
393           CkAbort("ERROR\n");
394       }
395
396       delete gmsg;
397     }
398
399
400         void check_and_compute() {
401                 compute_kernel();
402
403                 // calculate error
404                 // not being done right now since we are doing a fixed no. of iterations
405
406                 constrainBC();
407
408                 if (iterations % CKP_FREQ == 0 || iterations > MAX_ITER){
409 #ifdef CMK_MEM_CHECKPOINT
410                         contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
411 #elif CMK_MESSAGE_LOGGING
412                         if(iterations > MAX_ITER)
413                                 contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
414                         else
415                                 AtSync();
416 #else
417             contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
418 #endif
419         } else {
420                         doStep();
421         }
422     }
423
424         void ResumeFromSync(){
425                 doStep();
426         }
427
428
429     // Check to see if we have received all neighbor values yet
430     // If all neighbor values have been received, we update our values and proceed
431     void compute_kernel() {     //Gauss-Siedal compute
432         int i;
433   #pragma omp parallel for schedule(dynamic) 
434         for(i=1; i<blockDimX+1; ++i) {
435           //printf("[%d] did  %d iteration out of %d \n", omp_get_thread_num(), i, blockDimX+1); 
436           for(int j=1; j<blockDimY+1; ++j) {
437               for(int k=1; k<blockDimZ+1; ++k) {
438                   // update my value based on the surrounding values
439                    temperature[index(i, j, k)] = (temperature[index(i-1, j, k)] 
440                                             +  temperature[index(i+1, j, k)]
441                                             +  temperature[index(i, j-1, k)]
442                                             +  temperature[index(i, j+1, k)]
443                                             +  temperature[index(i, j, k-1)]
444                                             +  temperature[index(i, j, k+1)]
445                                             +  temperature[index(i, j, k)] ) * DIVIDEBY7;
446               }
447           }
448       }
449     }
450
451     // Enforce some boundary conditions
452     void constrainBC() {
453       // Heat left, top and front faces of each chare's block
454       for(int i=1; i<blockDimX+1; ++i)
455       for(int k=1; k<blockDimZ+1; ++k)
456               temperature[index(i, 1, k)] = 255.0;
457       for(int j=1; j<blockDimY+1; ++j)
458           for(int k=1; k<blockDimZ+1; ++k)
459               temperature[index(1, j, k)] = 255.0;
460       for(int i=1; i<blockDimX+1; ++i)
461           for(int j=1; j<blockDimY+1; ++j)
462               temperature[index(i, j, 1)] = 255.0;
463     }
464
465 };
466
467 /** \class JacobiMap
468  *
469  */
470
471 class JacobiMap : public CkArrayMap {
472   public:
473     int X, Y, Z;
474     int *mapping;
475
476     JacobiMap(int x, int y, int z) {
477       X = x; Y = y; Z = z;
478       mapping = new int[X*Y*Z];
479
480       // we are assuming that the no. of chares in each dimension is a 
481       // multiple of the torus dimension
482
483       TopoManager tmgr;
484       int dimNX, dimNY, dimNZ, dimNT;
485
486       dimNX = tmgr.getDimNX();
487       dimNY = tmgr.getDimNY();
488       dimNZ = tmgr.getDimNZ();
489       dimNT = tmgr.getDimNT();
490
491       // we are assuming that the no. of chares in each dimension is a 
492       // multiple of the torus dimension
493       int numCharesPerPe = X*Y*Z/CkNumPes();
494
495       int numCharesPerPeX = X / dimNX;
496       int numCharesPerPeY = Y / dimNY;
497       int numCharesPerPeZ = Z / dimNZ;
498       int pe = 0, pes = CkNumPes();
499
500 #if USE_BLOCK_RNDMAP
501       int used[pes];
502       for(int i=0; i<pes; i++)
503         used[i] = 0;
504 #endif
505
506       if(dimNT < 2) {   // one core per node
507         if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ); 
508         for(int i=0; i<dimNX; i++)
509           for(int j=0; j<dimNY; j++)
510             for(int k=0; k<dimNZ; k++)
511             {
512 #if USE_BLOCK_RNDMAP
513               pe = myrand(pes); 
514               while(used[pe]!=0) {
515                 pe = myrand(pes); 
516               }
517               used[pe] = 1;
518 #endif
519
520               for(int ci=i*numCharesPerPeX; ci<(i+1)*numCharesPerPeX; ci++)
521                 for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
522                   for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
523 #if USE_TOPOMAP
524                     mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k);
525 #elif USE_BLOCK_RNDMAP
526                     mapping[ci*Y*Z + cj*Z + ck] = pe;
527 #endif
528                   }
529             }
530       } else {          // multiple cores per node
531         // In this case, we split the chares in the X dimension among the
532         // cores on the same node. The strange thing I figured out is that
533         // doing this in the Z dimension is not as good.
534         numCharesPerPeX /= dimNT;
535         if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ);
536
537         for(int i=0; i<dimNX; i++)
538           for(int j=0; j<dimNY; j++)
539             for(int k=0; k<dimNZ; k++)
540               for(int l=0; l<dimNT; l++)
541                 for(int ci=(dimNT*i+l)*numCharesPerPeX; ci<(dimNT*i+l+1)*numCharesPerPeX; ci++)
542                   for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
543                     for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
544                       mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k, l);
545                     }
546       } // end of if
547
548       if(CkMyPe() == 0) CkPrintf("Map generated ... \n");
549     }
550
551     ~JacobiMap() { 
552       delete [] mapping;
553     }
554
555     int procNum(int, const CkArrayIndex &idx) {
556       int *index = (int *)idx.data();
557       return mapping[index[0]*Y*Z + index[1]*Z + index[2]]; 
558     }
559 };
560
561 class OmpInitializer : public CBase_OmpInitializer{
562 public:
563         OmpInitializer(int numThreads) {
564 #ifdef JACOBI_OPENMP
565           if (numThreads < 1) {
566             numThreads = 1; 
567           }
568           //omp_set_num_threads(numThreads);
569           if(CkMyPe() == 0)
570           {
571 #pragma omp parallel
572               { if(omp_get_thread_num() == 0 )
573                   CkPrintf("Computation loop will be parallelized"
574                " with %d OpenMP threads\n", omp_get_num_threads());
575               }
576           }
577 #endif
578           mainProxy.start();
579         }
580 };
581
582
583 #include "jacobi3d.def.h"