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