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