4b0b88a1de4f84421b1837471965533e0bb96c77
[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(CkCallback(CkIndex_Main::warmupIter(NULL), thisProxy), 1);
106   }
107
108   // Start timer after a few iterations for performance stability
109   void warmupIter(CkReductionMsg *msg) {
110     iterations++;
111     if (iterations == WARM_ITER) {
112       startTime = CmiWallTimer();
113       array.doStep(CkCallback(CkIndex_Main::report(NULL), thisProxy),
114                    maxiterations);
115     } else {
116       array.doStep(CkCallback(CkIndex_Main::warmupIter(NULL), thisProxy),
117                    1);
118     }
119     delete msg;
120   }
121
122   // Each worker reports back to here when it completes an iteration
123   // Asynchronous check on threshhold satisfaction
124   void report(CkReductionMsg *msg) {
125     iterations++;
126     maxdifference=((double *) msg->getData())[0];
127     delete msg;
128     if ( maxdifference - THRESHHOLD<0) 
129       done(true);
130     if (iterations > maxiterations)
131       done(false);
132   }
133
134   void done(bool success) {
135       if(success)
136         CkPrintf("Difference %.10g Satisfied Threshhold %.10g in %d Iterations\n", maxdifference,THRESHHOLD,iterations);
137       else
138         CkPrintf("Completed %d Iterations , Difference %lf fails threshhold\n", iterations,maxdifference);
139       endTime = CmiWallTimer();
140       CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(maxiterations-1-WARM_ITER));
141       CkExit();
142   }
143 };
144
145 /** \class Jacobi
146  *
147  */
148
149 class Jacobi: public CBase_Jacobi {
150   Jacobi_SDAG_CODE
151
152   public:
153   double *temperature;
154   double *new_temperature;
155   int imsg;
156   int iterations;
157   int numExpected;
158   int istart,ifinish,jstart,jfinish;
159   double maxdifference;
160   bool leftBound, rightBound, topBound, bottomBound;
161   // Constructor, initialize values
162   Jacobi() {
163     __sdag_init();
164     usesAtSync=CmiTrue;
165
166     int i, j;
167     // allocate a two dimensional array
168     temperature = new double[(blockDimX+2) * (blockDimY+2)];
169     new_temperature = new double[(blockDimX+2) * (blockDimY+2)];
170
171     for(i=0; i<blockDimX+2; ++i) {
172       for(j=0; j<blockDimY+2; ++j) {
173         temperature[index(i, j)] = 0.;
174       } 
175     }
176     imsg=0;
177     iterations = 0;
178     numExpected=0;
179     maxdifference=0.;
180     // determine border conditions
181     leftBound=rightBound=topBound=bottomBound=false;
182     istart=jstart=1;
183     ifinish=blockDimX+1;
184     jfinish=blockDimY+1;
185
186     if(thisIndex.x==0)
187       {
188         leftBound=true;
189         istart++;
190       }
191     else
192       numExpected++;
193
194     if(thisIndex.x==num_chare_x-1)
195       {
196         rightBound=true;
197         ifinish--;
198       }
199     else
200       numExpected++;
201
202     if(thisIndex.y==0)
203       {
204         topBound=true;
205         jstart++;
206       }
207     else
208       numExpected++;
209
210     if(thisIndex.y==num_chare_y-1)
211       {
212         bottomBound=true;
213         jfinish--;
214       }
215     else
216       numExpected++;
217     constrainBC();
218   }
219
220   void pup(PUP::er &p)
221   {
222     CBase_Jacobi::pup(p);
223     __sdag_pup(p);
224     p|imsg;
225     p|iterations;
226     p|numExpected;
227     p|maxdifference;
228     p|istart; p|ifinish; p|jstart; p|jfinish;
229     p|leftBound; p|rightBound; p|topBound; p|bottomBound;
230     
231     size_t size = (blockDimX+2) * (blockDimY+2);
232     if (p.isUnpacking()) {
233       temperature = new double[size];
234       new_temperature = new double[size];
235     }
236     p(temperature, size);
237     p(new_temperature, size);
238   }
239
240   Jacobi(CkMigrateMessage* m) {__sdag_init();}
241
242   ~Jacobi() { 
243     delete [] temperature; 
244     delete [] new_temperature; 
245   }
246
247   // Send ghost faces to the six neighbors
248   void begin_iteration(void) {
249     AtSync();
250     if (thisIndex.x == 0 && thisIndex.y == 0) {
251       CkPrintf("Start of iteration %d\n", iterations);
252     }
253     iterations++;
254
255     if(!leftBound)
256       {
257         double *leftGhost =  new double[blockDimY];
258         for(int j=0; j<blockDimY; ++j) 
259           leftGhost[j] = temperature[index(1, j+1)];
260         thisProxy(thisIndex.x-1, thisIndex.y)
261           .receiveGhosts(iterations, RIGHT, blockDimY, leftGhost);
262         delete [] leftGhost;
263       }
264     if(!rightBound)
265       {
266         double *rightGhost =  new double[blockDimY];
267         for(int j=0; j<blockDimY; ++j) 
268           rightGhost[j] = temperature[index(blockDimX, j+1)];
269         thisProxy(thisIndex.x+1, thisIndex.y)
270           .receiveGhosts(iterations, LEFT, blockDimY, rightGhost);
271         delete [] rightGhost;
272       }
273     if(!topBound)
274       {
275         double *topGhost =  new double[blockDimX];
276         for(int i=0; i<blockDimX; ++i) 
277           topGhost[i] = temperature[index(i+1, 1)];
278         thisProxy(thisIndex.x, thisIndex.y-1)
279           .receiveGhosts(iterations, BOTTOM, blockDimX, topGhost);
280         delete [] topGhost;
281       }
282     if(!bottomBound)
283       {
284         double *bottomGhost =  new double[blockDimX];
285         for(int i=0; i<blockDimX; ++i) 
286           bottomGhost[i] = temperature[index(i+1, blockDimY)];
287         thisProxy(thisIndex.x, thisIndex.y+1)
288           .receiveGhosts(iterations, TOP, blockDimX, bottomGhost);
289         delete [] bottomGhost;
290       }
291   }
292
293   void processGhosts(int dir, int size, double gh[]) {
294     switch(dir) {
295     case LEFT:
296       for(int j=0; j<size; ++j) {
297         temperature[index(0, j+1)] = gh[j];
298       }
299       break;
300     case RIGHT:
301       for(int j=0; j<size; ++j) {
302         temperature[index(blockDimX+1, j+1)] = gh[j];
303       }
304       break;
305     case TOP:
306       for(int i=0; i<size; ++i) {
307         temperature[index(i+1, 0)] = gh[i];
308       }
309       break;
310     case BOTTOM:
311       for(int i=0; i<size; ++i) {
312         temperature[index(i+1, blockDimY+1)] = gh[i];
313       }
314       break;
315     default:
316       CkAbort("ERROR\n");
317     }
318   }
319
320
321   void check_and_compute(CkCallback cb, int numSteps) {
322     compute_kernel();
323
324     double *tmp;
325     tmp = temperature;
326     temperature = new_temperature;
327     new_temperature = tmp;
328     constrainBC();
329
330     contribute(sizeof(double), &maxdifference, CkReduction::max_double, cb);
331     if (numSteps > 1)
332       doStep(cb, numSteps-1);
333   }
334
335
336   // When all neighbor values have been received, 
337   // we update our values and proceed
338   void compute_kernel() {
339     double temperatureIth=0.;
340     double difference=0.;
341     maxdifference=0.;
342 #pragma unroll    
343     for(int i=istart; i<ifinish; ++i) {
344       for(int j=jstart; j<jfinish; ++j) {
345         // calculate discrete mean value property 5 pt stencil
346         temperatureIth=(temperature[index(i, j)] 
347                         + temperature[index(i-1, j)] 
348                         +  temperature[index(i+1, j)]
349                         +  temperature[index(i, j-1)]
350                         +  temperature[index(i, j+1)]) * 0.2;
351
352         // update relative error
353         difference=temperatureIth - temperature[index(i, j)];
354         // fix sign without fabs overhead
355         if(difference<0) difference*=-1.0; 
356         maxdifference=(maxdifference>difference) ? maxdifference : difference;
357         new_temperature[index(i, j)] = temperatureIth;
358       }
359     }
360   }
361
362   // Enforce some boundary conditions
363   void constrainBC() {
364     if(topBound)
365       for(int i=0; i<blockDimX+2; ++i)
366         temperature[index(i, 1)] = 1.0;
367     if(leftBound)
368       for(int j=0; j<blockDimY+2; ++j)
369         temperature[index(1, j)] = 1.0;
370
371     if(bottomBound)
372       for(int i=0; i<blockDimX+2; ++i)
373         temperature[index(i, blockDimY)] = 0.;
374     if(rightBound)
375       for(int j=0; j<blockDimY+2; ++j)
376         temperature[index(blockDimX, j)] = 0.;
377
378   }
379   // for debugging
380  void dumpMatrix(double *matrix)
381   {
382     CkPrintf("[%d,%d]\n",thisIndex.x, thisIndex.y);
383     for(int i=0; i<blockDimX+2;++i)
384       {
385         for(int j=0; j<blockDimY+2;++j)
386           {
387             CkPrintf("%0.3lf ",matrix[index(i,j)]);
388           }
389         CkPrintf("\n");
390       }
391   }
392 };
393
394
395 #include "jacobi2d.def.h"