examples/charm++: Always use CkWallTimer, not CmiWallTimer (tested)
[charm.git] / examples / charm++ / jacobi2d-sdag / jacobi2d.C
1 /** \file jacobi2d.C
2  *  Author: Eric Bohm and Abhinav S Bhatele
3  *
4  *  This is jacobi3d-sdag cut down to 2d and fixed to be a correct
5  *  implementation of the finite difference method by Eric Bohm.
6  *
7  *  Date Created: Dec 7th, * 2010
8  *
9  */
10
11 #include "jacobi2d.decl.h"
12
13 // See README for documentation
14
15 /*readonly*/ CProxy_Main mainProxy;
16 /*readonly*/ int arrayDimX;
17 /*readonly*/ int arrayDimY;
18 /*readonly*/ int blockDimX;
19 /*readonly*/ int blockDimY;
20
21 // specify the number of worker chares in each dimension
22 /*readonly*/ int num_chare_x;
23 /*readonly*/ int num_chare_y;
24
25 /*readonly*/ int maxiterations;
26
27 static unsigned long next = 1;
28
29
30
31 #define index(a, b)     ( (b)*(blockDimX+2) + (a) )
32
33 #define MAX_ITER                100
34 #define WARM_ITER               2
35 #define LEFT                    1
36 #define RIGHT                   2
37 #define TOP                     3
38 #define BOTTOM                  4
39 #define DIVIDEBY5               0.2
40 const double THRESHHOLD =  0.001;
41
42 double startTime;
43 double endTime;
44
45 /** \class Main
46  *
47  */
48 class Main : public CBase_Main {
49 public:
50   CProxy_Jacobi array;
51   int iterations;
52   double maxdifference;
53   Main_SDAG_CODE;
54   Main(CkArgMsg* m) {
55     __sdag_init();
56     if ( (m->argc < 3) || (m->argc > 6)) {
57       CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
58       CkPrintf("OR %s [array_size] [block_size] maxiterations\n", m->argv[0]);
59       CkPrintf("OR %s [array_size_X] [array_size_Y] [block_size_X] [block_size_Y] \n", m->argv[0]);
60       CkPrintf("OR %s [array_size_X] [array_size_Y] [block_size_X] [block_size_Y] maxiterations\n", m->argv[0]);
61       CkAbort("Abort");
62     }
63
64     // set iteration counter to zero
65     iterations = 0;
66     // store the main proxy
67     mainProxy = thisProxy;
68         
69     if(m->argc <= 4) {
70       arrayDimX = arrayDimY = atoi(m->argv[1]);
71       blockDimX = blockDimY = atoi(m->argv[2]);
72     }
73     else if (m->argc >= 5) {
74       arrayDimX = atoi(m->argv[1]);
75       arrayDimY = atoi(m->argv[2]);
76       blockDimX = atoi(m->argv[3]); 
77       blockDimY = atoi(m->argv[4]); 
78     }
79     maxiterations=MAX_ITER;
80     if(m->argc==4)
81       maxiterations=atoi(m->argv[3]); 
82     if(m->argc==6)
83       maxiterations=atoi(m->argv[5]); 
84       
85     if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
86       CkAbort("array_size_X % block_size_X != 0!");
87     if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
88       CkAbort("array_size_Y % block_size_Y != 0!");
89
90     num_chare_x = arrayDimX / blockDimX;
91     num_chare_y = arrayDimY / blockDimY;
92
93     // print info
94     CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
95     CkPrintf("Running Jacobi on %d processors with (%d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y);
96     CkPrintf("Array Dimensions: %d %d\n", arrayDimX, arrayDimY);
97     CkPrintf("Block Dimensions: %d %d\n", blockDimX, blockDimY);
98     CkPrintf("max iterations %d\n", maxiterations);
99     CkPrintf("Threshhold %.10g\n", THRESHHOLD);
100       
101     // NOTE: boundary conditions must be set based on values
102       
103     // make proxy and populate array in one call
104     array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y);
105
106     // initiate computation
107     thisProxy.run();
108   }
109
110   void done(bool success) {
111       if(success)
112         CkPrintf("Difference %.10g Satisfied Threshhold %.10g in %d Iterations\n", maxdifference,THRESHHOLD,iterations);
113       else
114         CkPrintf("Completed %d Iterations , Difference %lf fails threshhold\n", iterations,maxdifference);
115       endTime = CkWallTimer();
116       CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(maxiterations-1-WARM_ITER));
117       CkExit();
118   }
119 };
120
121 /** \class Jacobi
122  *
123  */
124
125 class Jacobi: public CBase_Jacobi {
126   Jacobi_SDAG_CODE
127
128   public:
129   double *temperature;
130   double *new_temperature;
131   int imsg;
132   int iterations;
133   int numExpected;
134   int istart,ifinish,jstart,jfinish;
135   double maxdifference;
136   bool leftBound, rightBound, topBound, bottomBound;
137   // Constructor, initialize values
138   Jacobi() {
139     __sdag_init();
140     usesAtSync=CmiTrue;
141
142     int i, j;
143     // allocate a two dimensional array
144     temperature = new double[(blockDimX+2) * (blockDimY+2)];
145     new_temperature = new double[(blockDimX+2) * (blockDimY+2)];
146
147     for(i=0; i<blockDimX+2; ++i) {
148       for(j=0; j<blockDimY+2; ++j) {
149         temperature[index(i, j)] = 0.;
150       } 
151     }
152     imsg=0;
153     iterations = 0;
154     numExpected=0;
155     maxdifference=0.;
156     // determine border conditions
157     leftBound=rightBound=topBound=bottomBound=false;
158     istart=jstart=1;
159     ifinish=blockDimX+1;
160     jfinish=blockDimY+1;
161
162     if(thisIndex.x==0)
163       {
164         leftBound=true;
165         istart++;
166       }
167     else
168       numExpected++;
169
170     if(thisIndex.x==num_chare_x-1)
171       {
172         rightBound=true;
173         ifinish--;
174       }
175     else
176       numExpected++;
177
178     if(thisIndex.y==0)
179       {
180         topBound=true;
181         jstart++;
182       }
183     else
184       numExpected++;
185
186     if(thisIndex.y==num_chare_y-1)
187       {
188         bottomBound=true;
189         jfinish--;
190       }
191     else
192       numExpected++;
193     constrainBC();
194   }
195
196   void pup(PUP::er &p)
197   {
198     CBase_Jacobi::pup(p);
199     __sdag_pup(p);
200     p|imsg;
201     p|iterations;
202     p|numExpected;
203     p|maxdifference;
204     p|istart; p|ifinish; p|jstart; p|jfinish;
205     p|leftBound; p|rightBound; p|topBound; p|bottomBound;
206     
207     size_t size = (blockDimX+2) * (blockDimY+2);
208     if (p.isUnpacking()) {
209       temperature = new double[size];
210       new_temperature = new double[size];
211     }
212     p(temperature, size);
213     p(new_temperature, size);
214   }
215
216   Jacobi(CkMigrateMessage* m) {__sdag_init();}
217
218   ~Jacobi() { 
219     delete [] temperature; 
220     delete [] new_temperature; 
221   }
222
223   // Send ghost faces to the six neighbors
224   void begin_iteration(void) {
225     AtSync();
226     if (thisIndex.x == 0 && thisIndex.y == 0) {
227       CkPrintf("Start of iteration %d\n", iterations);
228     }
229     iterations++;
230
231     if(!leftBound)
232       {
233         double *leftGhost =  new double[blockDimY];
234         for(int j=0; j<blockDimY; ++j) 
235           leftGhost[j] = temperature[index(1, j+1)];
236         thisProxy(thisIndex.x-1, thisIndex.y)
237           .receiveGhosts(iterations, RIGHT, blockDimY, leftGhost);
238         delete [] leftGhost;
239       }
240     if(!rightBound)
241       {
242         double *rightGhost =  new double[blockDimY];
243         for(int j=0; j<blockDimY; ++j) 
244           rightGhost[j] = temperature[index(blockDimX, j+1)];
245         thisProxy(thisIndex.x+1, thisIndex.y)
246           .receiveGhosts(iterations, LEFT, blockDimY, rightGhost);
247         delete [] rightGhost;
248       }
249     if(!topBound)
250       {
251         double *topGhost =  new double[blockDimX];
252         for(int i=0; i<blockDimX; ++i) 
253           topGhost[i] = temperature[index(i+1, 1)];
254         thisProxy(thisIndex.x, thisIndex.y-1)
255           .receiveGhosts(iterations, BOTTOM, blockDimX, topGhost);
256         delete [] topGhost;
257       }
258     if(!bottomBound)
259       {
260         double *bottomGhost =  new double[blockDimX];
261         for(int i=0; i<blockDimX; ++i) 
262           bottomGhost[i] = temperature[index(i+1, blockDimY)];
263         thisProxy(thisIndex.x, thisIndex.y+1)
264           .receiveGhosts(iterations, TOP, blockDimX, bottomGhost);
265         delete [] bottomGhost;
266       }
267   }
268
269   void processGhosts(int dir, int size, double gh[]) {
270     switch(dir) {
271     case LEFT:
272       for(int j=0; j<size; ++j) {
273         temperature[index(0, j+1)] = gh[j];
274       }
275       break;
276     case RIGHT:
277       for(int j=0; j<size; ++j) {
278         temperature[index(blockDimX+1, j+1)] = gh[j];
279       }
280       break;
281     case TOP:
282       for(int i=0; i<size; ++i) {
283         temperature[index(i+1, 0)] = gh[i];
284       }
285       break;
286     case BOTTOM:
287       for(int i=0; i<size; ++i) {
288         temperature[index(i+1, blockDimY+1)] = gh[i];
289       }
290       break;
291     default:
292       CkAbort("ERROR\n");
293     }
294   }
295
296
297   void check_and_compute(CkCallback cb, int numSteps) {
298     compute_kernel();
299
300     double *tmp;
301     tmp = temperature;
302     temperature = new_temperature;
303     new_temperature = tmp;
304     constrainBC();
305
306     contribute(sizeof(double), &maxdifference, CkReduction::max_double, cb);
307     if (numSteps > 1)
308       doStep(cb, numSteps-1);
309   }
310
311
312   // When all neighbor values have been received, 
313   // we update our values and proceed
314   void compute_kernel() {
315     double temperatureIth=0.;
316     double difference=0.;
317     maxdifference=0.;
318 #pragma unroll    
319     for(int i=istart; i<ifinish; ++i) {
320       for(int j=jstart; j<jfinish; ++j) {
321         // calculate discrete mean value property 5 pt stencil
322         temperatureIth=(temperature[index(i, j)] 
323                         + temperature[index(i-1, j)] 
324                         +  temperature[index(i+1, j)]
325                         +  temperature[index(i, j-1)]
326                         +  temperature[index(i, j+1)]) * 0.2;
327
328         // update relative error
329         difference=temperatureIth - temperature[index(i, j)];
330         // fix sign without fabs overhead
331         if(difference<0) difference*=-1.0; 
332         maxdifference=(maxdifference>difference) ? maxdifference : difference;
333         new_temperature[index(i, j)] = temperatureIth;
334       }
335     }
336   }
337
338   // Enforce some boundary conditions
339   void constrainBC() {
340     if(topBound)
341       for(int i=0; i<blockDimX+2; ++i)
342         temperature[index(i, 1)] = 1.0;
343     if(leftBound)
344       for(int j=0; j<blockDimY+2; ++j)
345         temperature[index(1, j)] = 1.0;
346
347     if(bottomBound)
348       for(int i=0; i<blockDimX+2; ++i)
349         temperature[index(i, blockDimY)] = 0.;
350     if(rightBound)
351       for(int j=0; j<blockDimY+2; ++j)
352         temperature[index(blockDimX, j)] = 0.;
353
354   }
355   // for debugging
356  void dumpMatrix(double *matrix)
357   {
358     CkPrintf("[%d,%d]\n",thisIndex.x, thisIndex.y);
359     for(int i=0; i<blockDimX+2;++i)
360       {
361         for(int j=0; j<blockDimY+2;++j)
362           {
363             CkPrintf("%0.3lf ",matrix[index(i,j)]);
364           }
365         CkPrintf("\n");
366       }
367   }
368 };
369
370
371 #include "jacobi2d.def.h"