add Jacobi3d example without ft and with correct boundary constrain so that parallel...
[charm.git] / examples / charm++ / jacobi3d-sdag-constrain / jacobi3d.C
1 #define MAX_ITER        200     
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 #include "jacobi3d.decl.h"
12
13 /*readonly*/ CProxy_Main mainProxy;
14 /*readonly*/ int arrayDimX;
15 /*readonly*/ int arrayDimY;
16 /*readonly*/ int arrayDimZ;
17 /*readonly*/ int blockDimX;
18 /*readonly*/ int blockDimY;
19 /*readonly*/ int blockDimZ;
20
21 // specify the number of worker chares in each dimension
22 /*readonly*/ int num_chare_x;
23 /*readonly*/ int num_chare_y;
24 /*readonly*/ int num_chare_z;
25 /*readonly*/ int max_iter;
26
27 #define wrapX(a)        (((a)+num_chare_x)%num_chare_x)
28 #define wrapY(a)        (((a)+num_chare_y)%num_chare_y)
29 #define wrapZ(a)        (((a)+num_chare_z)%num_chare_z)
30
31 #define index(a,b,c)    ((a)+(b)*(blockDimX+2)+(c)*(blockDimX+2)*(blockDimY+2))
32
33 double startTime;
34 double endTime;
35
36 /** \class Main
37  *
38  */
39 class Main : public CBase_Main {
40 public:
41   CProxy_Jacobi array;
42   int iterations;
43
44   Main(CkArgMsg* m) {
45     if ( (m->argc != 3) && (m->argc != 7) && (m->argc != 4) && (m->argc != 8)) {
46       CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
47       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]);
48       CkAbort("Abort");
49     }
50
51     // set iteration counter to zero
52     iterations = 0;
53     max_iter = MAX_ITER;
54     // store the main proxy
55     mainProxy = thisProxy;
56         
57     if(m->argc <5 ) {
58       arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
59       blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]); 
60       if(m->argc == 4)
61           max_iter =  atoi(m->argv[3]);
62     }
63     else if (m->argc <9) {
64       arrayDimX = atoi(m->argv[1]);
65       arrayDimY = atoi(m->argv[2]);
66       arrayDimZ = atoi(m->argv[3]);
67       blockDimX = atoi(m->argv[4]); 
68       blockDimY = atoi(m->argv[5]); 
69       blockDimZ = atoi(m->argv[6]);
70       if(m->argc == 8)
71           max_iter =  atoi(m->argv[7]);
72     }
73
74     if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
75       CkAbort("array_size_X % block_size_X != 0!");
76     if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
77       CkAbort("array_size_Y % block_size_Y != 0!");
78     if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
79       CkAbort("array_size_Z % block_size_Z != 0!");
80
81     num_chare_x = arrayDimX / blockDimX;
82     num_chare_y = arrayDimY / blockDimY;
83     num_chare_z = arrayDimZ / blockDimZ;
84
85     // print info
86     CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
87     CkPrintf("Running Jacobi on %d processors with (%d, %d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y, num_chare_z);
88     CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
89     CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
90
91     // Create new array of worker chares
92     array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y, num_chare_z);
93
94     //Start the computation
95     array.run();
96     startTime = CkWallTimer();
97   }
98
99   void commonExit(int iter)
100   {
101     endTime = CkWallTimer();
102     CkPrintf("Time elapsed per iteration: %f total time %f \n", (endTime - startTime) / iter, (endTime-startTime));
103     CkExit();
104
105   }
106   void doneConverge(int done_iters) {
107       CkPrintf(" finished due to convergence %d \n", done_iters); 
108       commonExit(done_iters);
109   }
110   void doneIter(double error)
111   {
112       CkPrintf(" finished due to maximum iterations %d with error  %f \n", max_iter, error); 
113       commonExit(max_iter);
114   }
115 };
116
117 /** \class Jacobi
118  *
119  */
120
121 class Jacobi: public CBase_Jacobi {
122   Jacobi_SDAG_CODE
123
124 public:
125   int iterations;
126   int neighbors;
127   int remoteCount;
128   double error;
129   double *temperature;
130   double *new_temperature;
131   bool converged;
132
133   // Constructor, initialize values
134   Jacobi() {
135     converged = false;
136     neighbors = 6;
137     if(thisIndex.x == 0) 
138         neighbors--;
139     if( thisIndex.x== num_chare_x-1)
140         neighbors--;
141     if(thisIndex.y == 0) 
142         neighbors--;
143     if( thisIndex.y== num_chare_y-1)
144         neighbors--;
145     if(thisIndex.z == 0) 
146         neighbors--;
147     if( thisIndex.z== num_chare_z-1)
148           neighbors--;
149
150     // allocate a three dimensional array
151     temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
152     new_temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
153
154     for(int k=0; k<blockDimZ+2; ++k)
155       for(int j=0; j<blockDimY+2; ++j)
156         for(int i=0; i<blockDimX+2; ++i)
157           new_temperature[index(i, j, k)] = temperature[index(i, j, k)] = 0.0;
158     //print();
159     iterations = 0;
160     constrainBC();
161     //print();
162   }
163
164   void pup(PUP::er &p)
165   {
166     CBase_Jacobi::pup(p);
167     __sdag_pup(p);
168     p|iterations;
169     p|neighbors;
170
171     size_t size = (blockDimX+2) * (blockDimY+2) * (blockDimZ+2);
172     if (p.isUnpacking()) {
173       temperature = new double[size];
174       new_temperature = new double[size];
175     }
176     p(temperature, size);
177     p(new_temperature, size);
178   }
179
180   Jacobi(CkMigrateMessage* m) { }
181
182   ~Jacobi() { 
183     delete [] temperature; 
184     delete [] new_temperature; 
185   }
186
187   // Send ghost faces to the six neighbors
188   void begin_iteration(void) {
189     // Copy different faces into messages
190     double *leftGhost =  new double[blockDimY*blockDimZ];
191     double *rightGhost =  new double[blockDimY*blockDimZ];
192     double *topGhost =  new double[blockDimX*blockDimZ];
193     double *bottomGhost =  new double[blockDimX*blockDimZ];
194     double *frontGhost =  new double[blockDimX*blockDimY];
195     double *backGhost =  new double[blockDimX*blockDimY];
196     for(int k=0; k<blockDimZ; ++k)
197       for(int j=0; j<blockDimY; ++j) {
198         leftGhost[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
199         rightGhost[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
200       }
201
202     for(int k=0; k<blockDimZ; ++k)
203       for(int i=0; i<blockDimX; ++i) {
204         topGhost[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
205         bottomGhost[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
206       }
207
208     for(int j=0; j<blockDimY; ++j)
209       for(int i=0; i<blockDimX; ++i) {
210         frontGhost[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
211         backGhost[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
212       }
213
214     int x = thisIndex.x, y = thisIndex.y, z = thisIndex.z;
215     if(thisIndex.x>0)
216         thisProxy(wrapX(x-1),y,z).updateGhosts(iterations, RIGHT,  blockDimY, blockDimZ, rightGhost);
217     if(thisIndex.x<num_chare_x-1)
218         thisProxy(wrapX(x+1),y,z).updateGhosts(iterations, LEFT,   blockDimY, blockDimZ, leftGhost);
219     if(thisIndex.y>0)
220         thisProxy(x,wrapY(y-1),z).updateGhosts(iterations, TOP,    blockDimX, blockDimZ, topGhost);
221     if(thisIndex.y<num_chare_y-1)
222         thisProxy(x,wrapY(y+1),z).updateGhosts(iterations, BOTTOM, blockDimX, blockDimZ, bottomGhost);
223     if(thisIndex.z>0)
224         thisProxy(x,y,wrapZ(z-1)).updateGhosts(iterations, BACK,   blockDimX, blockDimY, backGhost);
225     if(thisIndex.z<num_chare_z-1)
226         thisProxy(x,y,wrapZ(z+1)).updateGhosts(iterations, FRONT,  blockDimX, blockDimY, frontGhost);
227
228     delete [] leftGhost;
229     delete [] rightGhost;
230     delete [] bottomGhost;
231     delete [] topGhost;
232     delete [] frontGhost;
233     delete [] backGhost;
234   }
235
236   void updateBoundary(int dir, int height, int width, double* gh) {
237     switch(dir) {
238     case LEFT:
239       for(int k=0; k<width; ++k)
240         for(int j=0; j<height; ++j) {
241           temperature[index(0, j+1, k+1)] = gh[k*height+j];
242         }
243       break;
244     case RIGHT:
245       for(int k=0; k<width; ++k)
246         for(int j=0; j<height; ++j) {
247           temperature[index(blockDimX+1, j+1, k+1)] = gh[k*height+j];
248         }
249       break;
250     case BOTTOM:
251       for(int k=0; k<width; ++k)
252         for(int i=0; i<height; ++i) {
253           temperature[index(i+1, 0, k+1)] = gh[k*height+i];
254         }
255       break;
256     case TOP:
257       for(int k=0; k<width; ++k)
258         for(int i=0; i<height; ++i) {
259           temperature[index(i+1, blockDimY+1, k+1)] = gh[k*height+i];
260         }
261       break;
262     case FRONT:
263       for(int j=0; j<width; ++j)
264         for(int i=0; i<height; ++i) {
265           temperature[index(i+1, j+1, 0)] = gh[j*height+i];
266         }
267       break;
268     case BACK:
269       for(int j=0; j<width; ++j)
270         for(int i=0; i<height; ++i) {
271           temperature[index(i+1, j+1, blockDimZ+1)] = gh[j*height+i];
272         }
273       break;
274     default:
275       CkAbort("ERROR\n");
276     }
277   }
278
279   // Check to see if we have received all neighbor values yet
280   // If all neighbor values have been received, we update our values and proceed
281   double computeKernel() {
282     double error = 0.0, max_error = 0.0;
283     for(int k=1; k<blockDimZ+1; ++k)
284       for(int j=1; j<blockDimY+1; ++j)
285         for(int i=1; i<blockDimX+1; ++i) {
286           // update my value based on the surrounding values
287           new_temperature[index(i, j, k)] = (temperature[index(i-1, j, k)] 
288                                              +  temperature[index(i+1, j, k)]
289                                              +  temperature[index(i, j-1, k)]
290                                              +  temperature[index(i, j+1, k)]
291                                              +  temperature[index(i, j, k-1)]
292                                              +  temperature[index(i, j, k+1)]
293                                              +  temperature[index(i, j, k)] ) * DIVIDEBY7;
294           error = fabs(new_temperature[index(i,j,k)] - temperature[index(i,j,k)]);
295           if (error > max_error) {
296             max_error = error;
297           }
298         } // end for
299     
300     double *tmp;
301     tmp = temperature;
302     temperature = new_temperature;
303     new_temperature = tmp;
304
305     //constrainBC();
306
307     return max_error;
308   }
309
310   void print()
311   {
312
313     for(int k=1; k<blockDimZ+2; ++k)
314       for(int j=1; j<blockDimY+2; ++j)
315         for(int i=1; i<blockDimX+2; ++i)
316           CkPrintf(" -%d:%d:%d %f ", k,j,i, temperature[index(k, j, i)]);
317     CkPrintf("--------------------------------\n");
318   }
319   // Enforce some boundary conditions
320   void constrainBC() {
321     // // Heat right, left
322     if(thisIndex.x == 0 )
323         for(int j=0; j<blockDimY+2; ++j)
324             for(int k=0; k<blockDimZ+2; ++k)
325             {   
326                 new_temperature[index(0, j, k)] = temperature[index(0, j, k)] = 255.0;
327             }
328     if(thisIndex.y == 0 )
329         for(int j=0; j<blockDimX+2; ++j)
330             for(int k=0; k<blockDimZ+2; ++k)
331             {   
332                 new_temperature[index(j,0, k)]  = temperature[index(j,0, k)] = 255.0;
333             }
334     if(thisIndex.z == 0 )
335         for(int j=0; j<blockDimX+2; ++j)
336             for(int k=0; k<blockDimY+2; ++k)
337             {   
338                 new_temperature[index(j, k, 0)] = temperature[index(j, k, 0)] = 255.0;
339             }
340
341   }
342 };
343
344 #include "jacobi3d.def.h"