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