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