bug fixes
[charm.git] / examples / bigsim / charm / jacobi2d / jacobi2d.C
1 /*****************************************************************************
2  * $Source$
3  * $Author$
4  * $Date$
5  * $Revision$
6  *****************************************************************************/
7
8 /** \file jacobi2d.C
9  *  Author: Abhinav S Bhatele
10  *  Date Created: October 24th, 2007
11  *
12  *
13  *    ***********  ^
14  *    *         *  |
15  *    *         *  |
16  *    *         *  X
17  *    *         *  |
18  *    *         *  |
19  *    ***********  ~
20  *    <--- Y --->
21  *
22  *    X: blockDimX, arrayDimX --> wrap_x
23  *    Y: blockDimY, arrayDimY --> wrap_y
24  */
25
26 #include "jacobi2d.decl.h"
27 #include "TopoManager.h"
28
29 /*readonly*/ CProxy_Main mainProxy;
30 /*readonly*/ int blockDimX;
31 /*readonly*/ int blockDimY;
32 /*readonly*/ int arrayDimX;
33 /*readonly*/ int arrayDimY;
34
35 // specify the number of worker chares in each dimension
36 /*readonly*/ int num_chare_x;
37 /*readonly*/ int num_chare_y;
38
39 /*readonly*/ int globalBarrier;
40 /*readonly*/ int localBarrier;
41
42 // We want to wrap entries around, and because mod operator % 
43 // sometimes misbehaves on negative values. -1 maps to the highest value.
44 #define wrap_x(a)  (((a)+num_chare_x)%num_chare_x)
45 #define wrap_y(a)  (((a)+num_chare_y)%num_chare_y)
46
47 #define MAX_ITER        25
48 #define WARM_ITER       5
49 #define LEFT            1
50 #define RIGHT           2
51 #define TOP             3
52 #define BOTTOM          4
53
54 double startTime;
55 double endTime;
56
57 class Main : public CBase_Main
58 {
59   public:
60     CProxy_Jacobi array;
61     int iterations;
62
63     Main(CkArgMsg* m) {
64       if ( (m->argc < 4) || (m->argc > 7) ) {
65         CkPrintf("%s [array_size] [block_size] +[no]barrier [+[no]local]\n", m->argv[0]);
66         CkPrintf("%s [array_size_X] [array_size_Y] [block_size_X] [block_size_Y] +[no]barrier [+[no]local]\n", m->argv[0]);
67         CkAbort("Abort");
68       }
69
70       if(m->argc < 6) {
71         arrayDimY = arrayDimX = atoi(m->argv[1]);
72         blockDimY = blockDimX = atoi(m->argv[2]);
73         if(strcasecmp(m->argv[3], "+nobarrier") == 0) {
74           globalBarrier = 0;
75           if(strcasecmp(m->argv[4], "+nolocal") == 0)
76             localBarrier = 0;
77           else
78             localBarrier = 1;
79         }
80         else
81           globalBarrier = 1;
82         if(globalBarrier == 0) CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
83       }
84       else {
85         arrayDimX = atoi(m->argv[1]);
86         arrayDimY = atoi(m->argv[2]);
87         blockDimX = atoi(m->argv[3]);
88         blockDimY = atoi(m->argv[4]);
89         if(strcasecmp(m->argv[5], "+nobarrier") == 0) {
90           globalBarrier = 0;
91           if(strcasecmp(m->argv[6], "+nolocal") == 0)
92             localBarrier = 0;
93           else
94             localBarrier = 1;
95         }
96         else
97           globalBarrier = 1;
98         if(globalBarrier == 0 && localBarrier == 0)
99           CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
100         else
101           CkPrintf("\nSTENCIL COMPUTATION WITH LOCAL BARRIERS\n");
102       }
103
104       if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
105         CkAbort("array_size_X % block_size_X != 0!");
106       if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
107         CkAbort("array_size_Y % block_size_Y != 0!");
108
109       // store the main proxy
110       mainProxy = thisProxy;
111
112       num_chare_x = arrayDimX / blockDimX;
113       num_chare_y = arrayDimY / blockDimY;
114
115       // print info
116       CkPrintf("Running Jacobi on %d processors with (%d, %d) elements\n", CkNumPes(), num_chare_x, num_chare_y);
117       CkPrintf("Array Dimensions: %d %d\n", arrayDimX, arrayDimY);
118       CkPrintf("Block Dimensions: %d %d\n", blockDimX, blockDimY);
119
120       // Create new array of worker chares
121 #if USE_TOPOMAP
122       CProxy_JacobiMap map = CProxy_JacobiMap::ckNew(num_chare_x, num_chare_y);
123       CkPrintf("Topology Mapping is being done ... \n");
124       CkArrayOptions opts(num_chare_x, num_chare_y);
125       opts.setMap(map);
126       array = CProxy_Jacobi::ckNew(opts);
127 #else
128       array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y);
129 #endif
130
131       //Start the computation
132       iterations = 0;
133       array.begin_iteration();
134     }
135
136     // Each worker reports back to here when it completes an iteration
137     void report(CkReductionMsg *msg) {
138       iterations++;
139       if(iterations == WARM_ITER)
140         startTime = CmiWallTimer();
141       double error = *((double *)msg->getData());
142
143       if((globalBarrier == 1 && iterations < MAX_ITER) || (globalBarrier == 0 && iterations <= WARM_ITER)) {
144         if(iterations > WARM_ITER) CkPrintf("Start of iteration %d\n", iterations);
145         BgPrintf("BgPrint> Start of iteration at %f\n");
146         array.begin_iteration();
147       } else {
148         CkPrintf("Completed %d iterations\n", MAX_ITER);
149         endTime = CmiWallTimer();
150         CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(MAX_ITER-WARM_ITER));
151         CkExit();
152       }
153     }
154
155 };
156
157 class Jacobi: public CBase_Jacobi {
158   public:
159     int arrived_left;
160     int arrived_right;
161     int arrived_top;
162     int arrived_bottom;
163     int readyToSend;
164
165     double **temperature;
166     double **new_temperature;
167     void *sendLogs[4];
168     void *ackLogs[5];
169     int iterations;
170
171     // Constructor, initialize values
172     Jacobi() {
173       int i,j;
174       // allocate two dimensional arrays
175       temperature = new double*[blockDimX+2];
176       new_temperature = new double*[blockDimX+2];
177       for (i=0; i<blockDimX+2; i++) {
178         temperature[i] = new double[blockDimY+2];
179         new_temperature[i] = new double[blockDimY+2];
180       }
181       for(i=0;i<blockDimX+2; i++) {
182         for(j=0;j<blockDimY+2; j++) {
183           temperature[i][j] = 0.5;
184           new_temperature[i][j] = 0.5;
185         }
186       }
187
188       arrived_left = 0;
189       arrived_right = 0;
190       arrived_top = 0;
191       arrived_bottom = 0;
192       readyToSend = 5;
193       iterations = 0;
194       constrainBC();
195     }
196
197     Jacobi(CkMigrateMessage* m) {}
198
199     ~Jacobi() { 
200       for (int i=0; i<blockDimX+2; i++) {
201         delete [] temperature[i];
202         delete [] new_temperature[i];
203       }
204       delete [] temperature; 
205       delete [] new_temperature; 
206     }
207
208     // Perform one iteration of work
209     void begin_iteration(void) {
210       if(localBarrier == 1 && iterations > WARM_ITER) {
211         _TRACE_BG_TLINE_END(&ackLogs[readyToSend]);
212         readyToSend++;
213       }
214
215       if(readyToSend == 5) {
216         if(thisIndex.x == 0 && thisIndex.y == 0  && (globalBarrier == 0 && iterations > WARM_ITER)) {
217           CkPrintf("Start of iteration %d\n", iterations);
218           BgPrintf("BgPrint> Start of iteration at %f\n");
219         }
220
221         if(localBarrier == 1 && iterations > WARM_ITER) {
222           void *curLog = NULL;
223           _TRACE_BG_END_EXECUTE(1);
224           _TRACE_BG_BEGIN_EXECUTE_NOMSG("start next iteration", &curLog);
225           for(int i=0; i<5; i++)
226             _TRACE_BG_ADD_BACKWARD_DEP(ackLogs[i]);
227         }
228         if(localBarrier == 1 && iterations >= WARM_ITER)  readyToSend = 0;
229         iterations++;
230
231         // Copy left column and right column into temporary arrays
232         double *left_edge = new double[blockDimX];
233         double *right_edge = new double[blockDimX];
234
235         for(int i=0; i<blockDimX; i++){
236             left_edge[i] = temperature[i+1][1];
237             right_edge[i] = temperature[i+1][blockDimY];
238         }
239
240         int x = thisIndex.x;
241         int y = thisIndex.y;
242
243         // Send my left edge
244         thisProxy(x, wrap_y(y-1)).receiveGhosts(RIGHT, blockDimX, left_edge);
245         // Send my right edge
246         thisProxy(x, wrap_y(y+1)).receiveGhosts(LEFT, blockDimX, right_edge);
247         // Send my top edge
248         thisProxy(wrap_x(x-1), y).receiveGhosts(BOTTOM, blockDimY, &temperature[1][1]);
249         // Send my bottom edge
250         thisProxy(wrap_x(x+1), y).receiveGhosts(TOP, blockDimY, &temperature[blockDimX][1]);
251
252         delete [] right_edge;
253         delete [] left_edge;
254
255       }
256     }
257
258     void receiveGhosts(int dir, int size, double gh[]) {
259       int i, j;
260       _TRACE_BG_TLINE_END(&sendLogs[dir-1]);
261
262       switch(dir) {
263         case LEFT:
264           arrived_left++;
265           for(i=0; i<size; i++)
266             temperature[i+1][0] = gh[i];
267           break;
268         case RIGHT:
269           arrived_right++;
270           for(i=0; i<size; i++)
271             temperature[i+1][blockDimY+1] = gh[i];
272           break;
273         case TOP:
274           arrived_top++;
275           for(j=0; j<size; j++)
276             temperature[0][j+1] = gh[j];
277           break;
278         case BOTTOM:
279           arrived_bottom++;
280           for(j=0; j<size; j++)
281             temperature[blockDimX+1][j+1] = gh[j];
282           break;
283         default:
284           CkAbort("ERROR\n");
285       }
286       check_and_compute();
287     }
288
289     void check_and_compute() {
290       double error = 0.0, max_error = 0.0;
291       void *curLog = NULL;
292
293       if (arrived_left >=1 && arrived_right >=1 && arrived_top >=1 && arrived_bottom >= 1) {
294         arrived_left--;
295         arrived_right--;
296         arrived_top--;
297         arrived_bottom--;
298
299         _TRACE_BG_END_EXECUTE(1);
300         _TRACE_BG_BEGIN_EXECUTE_NOMSG("start computation", &curLog);
301         for(int i=0; i<4; i++)
302           _TRACE_BG_ADD_BACKWARD_DEP(sendLogs[i]);
303
304         if(localBarrier == 1 && iterations > WARM_ITER) {
305           int x = thisIndex.x;
306           int y = thisIndex.y;
307           thisProxy(x, wrap_y(y-1)).begin_iteration();
308           thisProxy(x, wrap_y(y+1)).begin_iteration();
309           thisProxy(wrap_x(x-1), y).begin_iteration();
310           thisProxy(wrap_x(x+1), y).begin_iteration();
311         }
312
313         compute_kernel();       
314
315         for(int i=1; i<blockDimX+1; i++) {
316           for(int j=1; j<blockDimY+1; j++) {
317             error = fabs(new_temperature[i][j] - temperature[i][j]);
318             if(error > max_error) {
319               max_error = error;
320             }
321           }
322         }
323
324         double **tmp;
325         tmp = temperature;
326         temperature = new_temperature;
327         new_temperature = tmp;
328
329         constrainBC();
330
331         if(globalBarrier == 1 || (globalBarrier==0 && (iterations <= WARM_ITER || iterations >= MAX_ITER))) {
332           contribute(sizeof(double), &max_error, CkReduction::max_double,
333               CkCallback(CkIndex_Main::report(NULL), mainProxy));
334         } else {
335           begin_iteration();
336         }
337       }
338     }
339
340     // Check to see if we have received all neighbor values yet
341     // If all neighbor values have been received, we update our values and proceed
342     void compute_kernel()
343     {
344       for(int i=1; i<blockDimX+1; i++) {
345         for(int j=1; j<blockDimY+1; j++) {
346           // update my value based on the surrounding values
347           new_temperature[i][j] = (temperature[i-1][j]+temperature[i+1][j]+temperature[i][j-1]+temperature[i][j+1]+temperature[i][j]) * 0.2;
348         }
349       }
350     }
351
352     // Enforce some boundary conditions
353     void constrainBC()
354     {
355      if(thisIndex.y == 0 && thisIndex.x < num_chare_x/2) {
356         for(int i=1; i<=blockDimX; i++)
357           temperature[i][1] = 1.0;
358       }
359
360       if(thisIndex.x == num_chare_x-1 && thisIndex.y >= num_chare_y/2) {
361         for(int j=1; j<=blockDimY; j++)
362           temperature[blockDimX][j] = 0.0;
363       }
364     }
365
366 };
367
368 /** \class JacobiMap
369  *
370  */
371
372 class JacobiMap : public CkArrayMap {
373   public:
374     int X, Y;
375     int **mapping;
376
377     JacobiMap(int x, int y) {
378       X = x; Y = y;
379       int i, j;
380       mapping = new int*[x];
381       for (i=0; i<x; i++)
382         mapping[i] = new int[y];
383
384       TopoManager tmgr;
385       // we are assuming that the no. of chares in each dimension is a 
386       // multiple of the torus dimension
387       int dimNX = tmgr.getDimNX();
388       int dimNY = tmgr.getDimNY();
389       int dimNZ = tmgr.getDimNZ();
390       int dimNT = tmgr.getDimNT();
391
392       int numCharesPerPeX = x / dimNX;
393       int numCharesPerPeY = numCharesPerPeX / dimNY;
394       int numCharesPerPeZ = y / dimNZ;
395       int numCharesPerPeT = numCharesPerPeZ / dimNT;
396
397       if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d %d\n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ, numCharesPerPeT);
398
399       for(int i=0; i<dimNX; i++)
400         for(int j=0; j<dimNY; j++)
401           for(int k=0; k<dimNZ; k++)
402             for(int l=0; l<dimNT; l++)
403               for(int ci = i*numCharesPerPeX+j*numCharesPerPeY; ci < i*numCharesPerPeX+(j+1)*numCharesPerPeY; ci++)
404                 for(int cj = k*numCharesPerPeZ+l*numCharesPerPeT; cj < k*numCharesPerPeZ+(l+1)*numCharesPerPeT; cj++)
405                   mapping[ci][cj] = tmgr.coordinatesToRank(i, j, k, l);
406
407       /*if(CkMyPe() == 0) {
408         for(int ci=0; ci<x; ci++) {
409           for(int cj=0; cj<y; cj++) {
410             CkPrintf("%d ", mapping[ci][cj]);
411           }
412           CkPrintf("\n");
413         }
414       }*/
415
416     }
417
418     ~JacobiMap() {
419       for (int i=0; i<X; i++)
420         delete [] mapping[i];
421       delete [] mapping;
422     }
423
424     int procNum(int, const CkArrayIndex &idx) {
425       int *index = (int *)idx.data();
426       return mapping[index[0]][index[1]];
427     }
428 };
429
430 #include "jacobi2d.def.h"