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