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