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