Example Jacobi2D SDAG: clean up command-line argument handling
[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(CkArgMsg* m) {
54     if ( (m->argc < 3) || (m->argc > 6)) {
55       CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
56       CkPrintf("OR %s [array_size] [block_size] maxiterations\n", m->argv[0]);
57       CkPrintf("OR %s [array_size_X] [array_size_Y] [block_size_X] [block_size_Y] \n", m->argv[0]);
58       CkPrintf("OR %s [array_size_X] [array_size_Y] [block_size_X] [block_size_Y] maxiterations\n", m->argv[0]);
59       CkAbort("Abort");
60     }
61
62     // set iteration counter to zero
63     iterations = 0;
64     // store the main proxy
65     mainProxy = thisProxy;
66         
67     if(m->argc <= 4) {
68       arrayDimX = arrayDimY = atoi(m->argv[1]);
69       blockDimX = blockDimY = atoi(m->argv[2]);
70     }
71     else if (m->argc >= 5) {
72       arrayDimX = atoi(m->argv[1]);
73       arrayDimY = atoi(m->argv[2]);
74       blockDimX = atoi(m->argv[3]); 
75       blockDimY = atoi(m->argv[4]); 
76     }
77     maxiterations=MAX_ITER;
78     if(m->argc==4)
79       maxiterations=atoi(m->argv[3]); 
80     if(m->argc==6)
81       maxiterations=atoi(m->argv[5]); 
82       
83     if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
84       CkAbort("array_size_X % block_size_X != 0!");
85     if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
86       CkAbort("array_size_Y % block_size_Y != 0!");
87
88     num_chare_x = arrayDimX / blockDimX;
89     num_chare_y = arrayDimY / blockDimY;
90
91     // print info
92     CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
93     CkPrintf("Running Jacobi on %d processors with (%d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y);
94     CkPrintf("Array Dimensions: %d %d\n", arrayDimX, arrayDimY);
95     CkPrintf("Block Dimensions: %d %d\n", blockDimX, blockDimY);
96     CkPrintf("max iterations %d\n", maxiterations);
97     CkPrintf("Threshhold %.10g\n", THRESHHOLD);
98       
99     // NOTE: boundary conditions must be set based on values
100       
101     // make proxy and populate array in one call
102     array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y);
103       
104     // initiate computation
105     array.doStep();
106   }
107
108
109
110   // Start timer after a few iterations for performance stability
111   void warmupIter(CkReductionMsg *msg) {
112     iterations++;
113     if (iterations == WARM_ITER)
114       startTime = CmiWallTimer();
115     array.doStep();
116     delete msg;
117   }
118
119   // Each worker reports back to here when it completes an iteration
120   // Asynchronous check on threshhold satisfaction
121   void report(CkReductionMsg *msg) {
122
123     iterations++;
124     maxdifference=((double *) msg->getData())[0];
125     delete msg;
126     if ( maxdifference - THRESHHOLD<0) 
127       {
128         CkPrintf("boo Difference %.10g Satisfied Threshhold %.10g in %d Iterations\n", maxdifference,THRESHHOLD,iterations);
129         done(true);
130       }
131   }
132
133   // when iteration count met
134   void iterationsDone(CkReductionMsg *msg) {
135     iterations++;
136     delete msg;
137     done(false);
138   }
139
140   void done(bool success)
141     {
142       if(success)
143         CkPrintf("Difference %.10g Satisfied Threshhold %.10g in %d Iterations\n", maxdifference,THRESHHOLD,iterations);
144       else
145         CkPrintf("Completed %d Iterations , Difference %lf fails threshhold\n", iterations,maxdifference);
146       endTime = CmiWallTimer();
147       CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(maxiterations-1-WARM_ITER));
148       CkExit();
149     }
150
151 };
152
153 /** \class Jacobi
154  *
155  */
156
157 class Jacobi: public CBase_Jacobi {
158   Jacobi_SDAG_CODE
159
160   public:
161   double *temperature;
162   double *new_temperature;
163   int imsg;
164   int iterations;
165   int numExpected;
166   int istart,ifinish,jstart,jfinish;
167   double maxdifference;
168   bool leftBound, rightBound, topBound, bottomBound;
169   // Constructor, initialize values
170   Jacobi() {
171     __sdag_init();
172     usesAtSync=CmiTrue;
173
174     int i, j;
175     // allocate a two dimensional array
176     temperature = new double[(blockDimX+2) * (blockDimY+2)];
177     new_temperature = new double[(blockDimX+2) * (blockDimY+2)];
178
179     for(i=0; i<blockDimX+2; ++i) {
180       for(j=0; j<blockDimY+2; ++j) {
181         temperature[index(i, j)] = 0.;
182       } 
183     }
184     imsg=0;
185     iterations = 0;
186     numExpected=0;
187     maxdifference=0.;
188     // determine border conditions
189     leftBound=rightBound=topBound=bottomBound=false;
190     istart=jstart=1;
191     ifinish=blockDimX+1;
192     jfinish=blockDimY+1;
193
194     if(thisIndex.x==0)
195       {
196         leftBound=true;
197         istart++;
198       }
199     else
200       numExpected++;
201
202     if(thisIndex.x==num_chare_x-1)
203       {
204         rightBound=true;
205         ifinish--;
206       }
207     else
208       numExpected++;
209
210     if(thisIndex.y==0)
211       {
212         topBound=true;
213         jstart++;
214       }
215     else
216       numExpected++;
217
218     if(thisIndex.y==num_chare_y-1)
219       {
220         bottomBound=true;
221         jfinish--;
222       }
223     else
224       numExpected++;
225     constrainBC();
226   }
227
228   void pup(PUP::er &p)
229   {
230     CBase_Jacobi::pup(p);
231     __sdag_pup(p);
232     p|imsg;
233     p|iterations;
234     p|numExpected;
235     p|maxdifference;
236     p|istart; p|ifinish; p|jstart; p|jfinish;
237     p|leftBound; p|rightBound; p|topBound; p|bottomBound;
238     
239     size_t size = (blockDimX+2) * (blockDimY+2);
240     if (p.isUnpacking()) {
241       temperature = new double[size];
242       new_temperature = new double[size];
243     }
244     p(temperature, size);
245     p(new_temperature, size);
246   }
247
248   Jacobi(CkMigrateMessage* m) {__sdag_init();}
249
250   ~Jacobi() { 
251     delete [] temperature; 
252     delete [] new_temperature; 
253   }
254
255   // Send ghost faces to the six neighbors
256   void begin_iteration(void) {
257     AtSync();
258     if (thisIndex.x == 0 && thisIndex.y == 0) {
259       CkPrintf("Start of iteration %d\n", iterations);
260     }
261     iterations++;
262
263     if(!leftBound)
264       {
265         double *leftGhost =  new double[blockDimY];
266         for(int j=0; j<blockDimY; ++j) 
267           leftGhost[j] = temperature[index(1, j+1)];
268         thisProxy(thisIndex.x-1, thisIndex.y)
269           .receiveGhosts(iterations, RIGHT, blockDimY, leftGhost);
270         delete [] leftGhost;
271       }
272     if(!rightBound)
273       {
274         double *rightGhost =  new double[blockDimY];
275         for(int j=0; j<blockDimY; ++j) 
276           rightGhost[j] = temperature[index(blockDimX, j+1)];
277         thisProxy(thisIndex.x+1, thisIndex.y)
278           .receiveGhosts(iterations, LEFT, blockDimY, rightGhost);
279         delete [] rightGhost;
280       }
281     if(!topBound)
282       {
283         double *topGhost =  new double[blockDimX];
284         for(int i=0; i<blockDimX; ++i) 
285           topGhost[i] = temperature[index(i+1, 1)];
286         thisProxy(thisIndex.x, thisIndex.y-1)
287           .receiveGhosts(iterations, BOTTOM, blockDimX, topGhost);
288         delete [] topGhost;
289       }
290     if(!bottomBound)
291       {
292         double *bottomGhost =  new double[blockDimX];
293         for(int i=0; i<blockDimX; ++i) 
294           bottomGhost[i] = temperature[index(i+1, blockDimY)];
295         thisProxy(thisIndex.x, thisIndex.y+1)
296           .receiveGhosts(iterations, TOP, blockDimX, bottomGhost);
297         delete [] bottomGhost;
298       }
299   }
300
301   void processGhosts(int dir, int size, double gh[]) {
302     switch(dir) {
303     case LEFT:
304       for(int j=0; j<size; ++j) {
305         temperature[index(0, j+1)] = gh[j];
306       }
307       break;
308     case RIGHT:
309       for(int j=0; j<size; ++j) {
310         temperature[index(blockDimX+1, j+1)] = gh[j];
311       }
312       break;
313     case TOP:
314       for(int i=0; i<size; ++i) {
315         temperature[index(i+1, 0)] = gh[i];
316       }
317       break;
318     case BOTTOM:
319       for(int i=0; i<size; ++i) {
320         temperature[index(i+1, blockDimY+1)] = gh[i];
321       }
322       break;
323     default:
324       CkAbort("ERROR\n");
325     }
326   }
327
328
329   void check_and_compute() {
330     compute_kernel();
331
332     double *tmp;
333     tmp = temperature;
334     temperature = new_temperature;
335     new_temperature = tmp;
336     constrainBC();
337     if (iterations <= WARM_ITER)
338       contribute(sizeof(double), &maxdifference, CkReduction::max_double, CkCallback(CkIndex_Main::warmupIter(NULL), mainProxy));
339     else if(iterations >= maxiterations)
340       contribute(sizeof(double), &maxdifference, CkReduction::max_double, CkCallback(CkIndex_Main::iterationsDone(NULL), mainProxy));
341     else
342       {  // nonblocking report 
343         contribute(sizeof(double), &maxdifference, CkReduction::max_double, CkCallback(CkIndex_Main::report(NULL), mainProxy));
344         doStep();
345       }
346   }
347
348
349   // When all neighbor values have been received, 
350   // we update our values and proceed
351   void compute_kernel() {
352     double temperatureIth=0.;
353     double difference=0.;
354     maxdifference=0.;
355 #pragma unroll    
356     for(int i=istart; i<ifinish; ++i) {
357       for(int j=jstart; j<jfinish; ++j) {
358         // calculate discrete mean value property 5 pt stencil
359         temperatureIth=(temperature[index(i, j)] 
360                         + temperature[index(i-1, j)] 
361                         +  temperature[index(i+1, j)]
362                         +  temperature[index(i, j-1)]
363                         +  temperature[index(i, j+1)]) * 0.2;
364
365         // update relative error
366         difference=temperatureIth - temperature[index(i, j)];
367         // fix sign without fabs overhead
368         if(difference<0) difference*=-1.0; 
369         maxdifference=(maxdifference>difference) ? maxdifference : difference;
370         new_temperature[index(i, j)] = temperatureIth;
371       }
372     }
373   }
374
375   // Enforce some boundary conditions
376   void constrainBC() {
377     if(topBound)
378       for(int i=0; i<blockDimX+2; ++i)
379         temperature[index(i, 1)] = 1.0;
380     if(leftBound)
381       for(int j=0; j<blockDimY+2; ++j)
382         temperature[index(1, j)] = 1.0;
383
384     if(bottomBound)
385       for(int i=0; i<blockDimX+2; ++i)
386         temperature[index(i, blockDimY)] = 0.;
387     if(rightBound)
388       for(int j=0; j<blockDimY+2; ++j)
389         temperature[index(blockDimX, j)] = 0.;
390
391   }
392   // for debugging
393  void dumpMatrix(double *matrix)
394   {
395     CkPrintf("[%d,%d]\n",thisIndex.x, thisIndex.y);
396     for(int i=0; i<blockDimX+2;++i)
397       {
398         for(int j=0; j<blockDimY+2;++j)
399           {
400             CkPrintf("%0.3lf ",matrix[index(i,j)]);
401           }
402         CkPrintf("\n");
403       }
404   }
405 };
406
407
408 #include "jacobi2d.def.h"