fix bug and clean up jacobi examples
[charm.git] / examples / charm++ / jacobi3d-2d-decomposition / jacobi3d.C
1 #define MAX_ITER        100     
2 #define LEFT                    1
3 #define RIGHT                   2
4 #define TOP                     3
5 #define BOTTOM                  4
6 #define FRONT                   5
7 #define BACK                    6
8 #define DIVIDEBY7               0.14285714285714285714
9 #define DELTA                   0.01
10
11 const double THRESHOLD   =  0.004;
12 #include "jacobi3d.decl.h"
13
14 /*readonly*/ CProxy_Main mainProxy;
15 /*readonly*/ int arrayDimX;
16 /*readonly*/ int arrayDimY;
17 /*readonly*/ int arrayDimZ;
18 /*readonly*/ int blockDimX;
19 /*readonly*/ int blockDimY;
20 /*readonly*/ int blockDimZ;
21
22 // specify the number of worker chares in each dimension
23 /*readonly*/ int numChareX;
24 /*readonly*/ int numChareY;
25 /*readonly*/ int numChareZ;
26 /*readonly*/ int maxiterations;
27
28 class Main : public CBase_Main {
29   double startTime;
30   double endTime;
31 public:
32   CProxy_Jacobi array;
33   int iterations;
34
35   Main(CkArgMsg* m) {
36     if ( (m->argc<3) || (m->argc>8) ) {
37       CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
38       CkPrintf("OR %s [array_size] [block_size] maxiterations\n", m->argv[0]);
39       CkPrintf("OR %s [array_size_X] [array_size_Y] [array_size_Z] [block_size_X] [block_size_Y] [block_size_Z]\n", m->argv[0]);
40       CkPrintf("OR %s [array_size_X] [array_size_Y] [array_size_Z] [block_size_X] [block_size_Y] [block_size_Z] maxiterations\n", m->argv[0]);
41       CkAbort("Abort");
42     }
43
44     iterations = 0;
45     // store the main proxy
46     mainProxy = thisProxy;
47
48     if(m->argc <=4 ) {
49       arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
50       blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]); 
51     }
52     else if (m->argc <=8) {
53       arrayDimX = atoi(m->argv[1]);
54       arrayDimY = atoi(m->argv[2]);
55       arrayDimZ = atoi(m->argv[3]);
56       blockDimX = atoi(m->argv[4]); 
57       blockDimY = atoi(m->argv[5]); 
58       blockDimZ = atoi(m->argv[6]);
59     }
60
61     maxiterations = MAX_ITER;
62     if(m->argc==4)
63       maxiterations = atoi(m->argv[3]); 
64     if(m->argc==8)
65       maxiterations = atoi(m->argv[7]); 
66
67     if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
68       CkAbort("array_size_X % block_size_X != 0!");
69     if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
70       CkAbort("array_size_Y % block_size_Y != 0!");
71     if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
72       CkAbort("array_size_Z % block_size_Z != 0!");
73
74     numChareX = arrayDimX / blockDimX;
75     numChareY = arrayDimY / blockDimY;
76     numChareZ = arrayDimZ / blockDimZ;
77
78     // print info
79     CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
80     CkPrintf("Running Jacobi on %d processors with (%d, %d, %d) chares\n", CkNumPes(), numChareX, numChareY, numChareZ);
81     CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
82     CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
83
84     // Create new array of worker chares
85     array = CProxy_Jacobi::ckNew(numChareX, numChareY, numChareZ);
86
87     //Start the computation
88     array.run();
89     startTime = CkWallTimer();
90   }
91
92   void done(int totalIter)
93   {
94     if(totalIter >= maxiterations)
95       CkPrintf("Finish due to max iterations %d, total time %.3f seconds. \n", totalIter, CkWallTimer()-startTime); 
96     else
97       CkPrintf("Finish due to convergence, iterations %d, total time %.3f seconds. \n", totalIter, CkWallTimer()-startTime); 
98     CkExit();
99   }
100 };
101
102
103 class Jacobi: public CBase_Jacobi {
104   Jacobi_SDAG_CODE
105
106 public:
107     int iterations;
108     int neighbors;
109     int remoteCount;
110     double error;
111     double ***temperature;
112     double ***new_temperature;
113     int converged;
114     bool leftBound, rightBound, topBound, bottomBound, frontBound, backBound;
115     int istart, jstart, kstart, ifinish, jfinish, kfinish;
116     double max_error;
117
118     // Constructor, initialize values
119     Jacobi() {
120       converged = 0;
121       neighbors = 0;
122       istart=jstart=kstart=1;
123       ifinish=blockDimX+1;
124       jfinish=blockDimY+1;
125       kfinish=blockDimZ+1;
126
127       leftBound = rightBound = topBound = bottomBound = frontBound = backBound = false;
128       if(thisIndex.x == 0) {
129         leftBound = true;
130         istart++;
131       }
132       else
133         neighbors++;
134       if(thisIndex.x == numChareX-1) {
135         rightBound = true;
136         ifinish--;
137       }
138       else
139         neighbors++;
140
141       if(thisIndex.y == 0) {
142         bottomBound = true;
143         jstart++;
144       }
145       else
146         neighbors++;
147
148       if(thisIndex.y == numChareY-1) {
149         topBound = true;
150         jfinish--;
151       }
152       else
153         neighbors++;
154
155       if(thisIndex.z == 0) {
156         backBound = true;
157         kstart++;
158       }
159       else
160         neighbors++;
161
162       if(thisIndex.z == numChareZ-1) {
163         frontBound = true;
164         kfinish--;
165       }
166       else
167         neighbors++;
168
169       // allocate a three dimensional array
170       temperature = new double**[blockDimX+2];
171       new_temperature = new double**[blockDimX+2];
172       for(int i=0; i<blockDimX+2; i++)
173       {
174         temperature[i] = new double*[blockDimY+2];
175         new_temperature[i] = new double*[blockDimY+2];
176         for(int j=0; j<blockDimY+2; j++)
177         {
178           temperature[i][j] = new double[blockDimZ+2];
179           new_temperature[i][j] = new double[blockDimZ+2];
180         }
181       }
182
183       for(int k=0; k<blockDimZ+2; ++k)
184         for(int j=0; j<blockDimY+2; ++j)
185           for(int i=0; i<blockDimX+2; ++i)
186             new_temperature[i][j][k] = temperature[i][j][k] = 0.0;
187       iterations = 0;
188       constrainBC();
189     }
190
191     void pup(PUP::er &p)
192     {
193       CBase_Jacobi::pup(p);
194       __sdag_pup(p);
195       p|iterations;
196       p|neighbors;
197
198       if (p.isUnpacking()) {
199         // allocate a three dimensional array
200         temperature = new double**[blockDimX+2];
201         new_temperature = new double**[blockDimX+2];
202         for(int i=0; i<blockDimX+2; i++)
203         {
204           temperature[i] = new double*[blockDimY+2];
205           new_temperature[i] = new double*[blockDimY+2];
206           for(int j=0; j<blockDimY+2; j++)
207           {
208             temperature[i][j] = new double[blockDimZ+2];
209             new_temperature[i][j] = new double[blockDimZ+2];
210           }
211         }
212       }
213
214       for(int k=0; k<blockDimZ+2; ++k)
215         for(int j=0; j<blockDimY+2; ++j)
216           for(int i=0; i<blockDimX+2; ++i)
217           {
218             p|new_temperature[i][j][k];
219             p|temperature[i][j][k];
220           }
221       iterations = 0;
222     }
223
224     Jacobi(CkMigrateMessage* m) { }
225
226     ~Jacobi() { 
227       delete [] temperature; 
228       delete [] new_temperature; 
229     }
230
231     // Send ghost faces to the six neighbors
232     void begin_iteration(void) {
233       iterations++;
234       // Copy different faces into messages
235       double *leftGhost =  new double[blockDimY*blockDimZ];
236       double *rightGhost =  new double[blockDimY*blockDimZ];
237       double *topGhost =  new double[blockDimX*blockDimZ];
238       double *bottomGhost =  new double[blockDimX*blockDimZ];
239       double *frontGhost =  new double[blockDimX*blockDimY];
240       double *backGhost =  new double[blockDimX*blockDimY];
241       for(int k=0; k<blockDimZ; ++k)
242         for(int j=0; j<blockDimY; ++j) {
243           leftGhost[k*blockDimY+j] = temperature[1][j+1][k+1];
244           rightGhost[k*blockDimY+j] = temperature[blockDimX][j+1][k+1];
245         }
246
247       for(int k=0; k<blockDimZ; ++k)
248         for(int i=0; i<blockDimX; ++i) {
249           bottomGhost[k*blockDimX+i] = temperature[i+1][1][k+1];
250           topGhost[k*blockDimX+i] = temperature[i+1][blockDimY][k+1];
251         }
252
253       for(int j=0; j<blockDimY; ++j)
254         for(int i=0; i<blockDimX; ++i) {
255           backGhost[j*blockDimX+i] = temperature[i+1][j+1][1];
256           frontGhost[j*blockDimX+i] = temperature[i+1][j+1][blockDimZ];
257         }
258
259       int x = thisIndex.x, y = thisIndex.y, z = thisIndex.z;
260       if(!leftBound)
261         thisProxy(x-1,y,z).receiveGhosts(iterations, RIGHT,  blockDimY, blockDimZ, leftGhost);
262       if(!rightBound)
263         thisProxy(x+1,y,z).receiveGhosts(iterations, LEFT,   blockDimY, blockDimZ, rightGhost);
264       if(!topBound)
265         thisProxy(x,y+1,z).receiveGhosts(iterations, BOTTOM,    blockDimX, blockDimZ, topGhost);
266       if(!bottomBound)
267         thisProxy(x,y-1,z).receiveGhosts(iterations, TOP, blockDimX, blockDimZ, bottomGhost);
268       if(!frontBound)
269         thisProxy(x,y,z+1).receiveGhosts(iterations, BACK,   blockDimX, blockDimY, frontGhost);
270       if(!backBound)
271         thisProxy(x,y,z-1).receiveGhosts(iterations, FRONT,  blockDimX, blockDimY, backGhost);
272
273       delete [] leftGhost;
274       delete [] rightGhost;
275       delete [] bottomGhost;
276       delete [] topGhost;
277       delete [] frontGhost;
278       delete [] backGhost;
279     }
280
281     void processGhosts(int dir, int height, int width, double* gh) {
282       switch(dir) {
283       case LEFT:
284         for(int k=0; k<width; ++k)
285           for(int j=0; j<height; ++j) {
286             temperature[0][j+1][k+1] = gh[k*height+j];
287           }
288         break;
289       case RIGHT:
290         for(int k=0; k<width; ++k)
291           for(int j=0; j<height; ++j) {
292             temperature[blockDimX+1][j+1][k+1] = gh[k*height+j];
293           }
294         break;
295       case BOTTOM:
296         for(int k=0; k<width; ++k)
297           for(int i=0; i<height; ++i) {
298             temperature[i+1][0][k+1] = gh[k*height+i];
299           }
300         break;
301       case TOP:
302         for(int k=0; k<width; ++k)
303           for(int i=0; i<height; ++i) {
304             temperature[i+1][blockDimY+1][k+1] = gh[k*height+i];
305           }
306         break;
307       case FRONT:
308         for(int j=0; j<width; ++j)
309           for(int i=0; i<height; ++i) {
310             temperature[i+1][j+1][blockDimZ+1] = gh[j*height+i];
311           }
312         break;
313       case BACK:
314         for(int j=0; j<width; ++j)
315           for(int i=0; i<height; ++i) {
316             temperature[i+1][j+1][0] = gh[j*height+i];
317           }
318         break;
319       default:
320         CkAbort("ERROR\n");
321       }
322     }
323
324     // Check to see if we have received all neighbor values yet
325     void check_and_compute() {
326       double error = 0.0;
327       max_error = 0.0;
328       for(int i=istart; i<ifinish; ++i) 
329         for(int j=jstart; j<jfinish; ++j)
330           for(int k=kstart; k<kfinish; ++k){
331             // update my value based on the surrounding values
332             new_temperature[i][j][k] = (temperature[i-1][j][k] 
333               +  temperature[i+1][j][k]
334               +  temperature[i][j-1][k]
335               +  temperature[i][j+1][k]
336               +  temperature[i][j][k-1]
337               +  temperature[i][j][k+1]
338               +  temperature[i][j][k] ) * DIVIDEBY7;
339             error = fabs(new_temperature[i][j][k] - temperature[i][j][k]);
340             if (error > max_error) {
341               max_error = error;
342             }
343           } // end for
344
345       double ***tmp;
346       tmp = temperature;
347       temperature = new_temperature;
348       new_temperature = tmp;
349       //dumpMatrix();
350     }
351
352     void dumpMatrix()
353     {
354
355       if(thisIndex.x + thisIndex.y + thisIndex.z == 0)
356       {
357         CkPrintf("\n\n[%d:%d:%d]\n", thisIndex.x, thisIndex.y, thisIndex.z);
358         for(int i=0; i<blockDimX+2; ++i){
359           for(int j=0; j<blockDimY+2; ++j){
360             for(int k=0; k<blockDimZ+2; ++k){
361               CkPrintf(" [%d:%d:%d %.3f] ", i,j,k, temperature[i][j][k]);
362             }
363             CkPrintf("\n");
364           }
365           CkPrintf("\n");
366         }
367         CkPrintf("\n\n");
368       }
369     }
370     // Enforce some boundary conditions
371     void constrainBC() {
372       if(leftBound)
373         for(int j=0; j<blockDimY+2; ++j)
374           for(int k=0; k<blockDimZ+2; ++k)
375           {   
376             new_temperature[1][j][k] = temperature[1][j][k] = 1.0;
377           }
378       if(rightBound)
379         for(int j=0; j<blockDimY+2; ++j)
380           for(int k=0; k<blockDimZ+2; ++k)
381           {   
382             new_temperature[blockDimX][j][k] = temperature[blockDimX][j][k] = 1.0;
383           }
384
385       if(topBound)
386         for(int i=0; i<blockDimX+2; ++i)
387           for(int k=0; k<blockDimZ+2; ++k)
388           {   
389             new_temperature[i][blockDimY][k]  = temperature[i][blockDimY][k] = 1.0;
390           }
391       if(bottomBound)
392         for(int i=0; i<blockDimX+2; ++i)
393           for(int k=0; k<blockDimZ+2; ++k)
394           {   
395             new_temperature[i][1][k]  = temperature[i][1][k] = 1.0;
396           }
397
398       if(frontBound)
399         for(int i=0; i<blockDimX+2; ++i)
400           for(int j=0; j<blockDimY+2; ++j)
401           {   
402             new_temperature[i][j][blockDimZ] = temperature[i][j][blockDimZ] = 1.0;
403           }
404
405       if(backBound)
406         for(int i=0; i<blockDimX+2; ++i)
407           for(int j=0; j<blockDimY+2; ++j)
408           {   
409             new_temperature[i][j][1] = temperature[i][j][1] = 1.0;
410           }
411     }
412 };
413
414 #include "jacobi3d.def.h"