stencil3d ldb example: moved barrier from main to the array
[charm.git] / examples / charm++ / load_balancing / stencil3d / stencil3d.C
1 /** \file stencil3d.C
2  *  Author: Abhinav S Bhatele
3  *  Date Created: December 28th, 2010
4  *
5  *  This example is written to be used with periodic measurement-based load
6  *  balancers at sync.
7  *
8  *
9  *
10  *            *****************
11  *         *               *  *
12  *   ^  *****************     *
13  *   |  *               *     *
14  *   |  *               *     *
15  *   |  *               *     *
16  *   Y  *               *     *
17  *   |  *               *     *
18  *   |  *               *     *
19  *   |  *               *  * 
20  *   ~  *****************    Z
21  *      <------ X ------> 
22  *
23  *   X: left, right --> wrap_x
24  *   Y: top, bottom --> wrap_y
25  *   Z: front, back --> wrap_z
26  */
27
28 #include "stencil3d.decl.h"
29 #include "TopoManager.h"
30
31 /*readonly*/ CProxy_Main mainProxy;
32 /*readonly*/ int arrayDimX;
33 /*readonly*/ int arrayDimY;
34 /*readonly*/ int arrayDimZ;
35 /*readonly*/ int blockDimX;
36 /*readonly*/ int blockDimY;
37 /*readonly*/ int blockDimZ;
38
39 // specify the number of worker chares in each dimension
40 /*readonly*/ int num_chare_x;
41 /*readonly*/ int num_chare_y;
42 /*readonly*/ int num_chare_z;
43
44 static unsigned long next = 1;
45
46 int myrand(int numpes) {
47   next = next * 1103515245 + 12345;
48   return((unsigned)(next/65536) % numpes);
49 }
50
51 // We want to wrap entries around, and because mod operator % 
52 // sometimes misbehaves on negative values. -1 maps to the highest value.
53 #define wrap_x(a)       (((a)+num_chare_x)%num_chare_x)
54 #define wrap_y(a)       (((a)+num_chare_y)%num_chare_y)
55 #define wrap_z(a)       (((a)+num_chare_z)%num_chare_z)
56
57 #define index(a,b,c)    ((a)+(b)*(blockDimX+2)+(c)*(blockDimX+2)*(blockDimY+2))
58
59 #define MAX_ITER        26
60 #define LBPERIOD        5
61 #define LEFT            1
62 #define RIGHT           2
63 #define TOP             3
64 #define BOTTOM          4
65 #define FRONT           5
66 #define BACK            6
67 #define DIVIDEBY7       0.14285714285714285714
68
69 double startTime;
70 double endTime;
71
72 /** \class Main
73  *
74  */
75 class Main : public CBase_Main {
76   public:
77     CProxy_Stencil array;
78
79     Main(CkArgMsg* m) {
80       if ( (m->argc != 3) && (m->argc != 7) ) {
81         CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
82         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]);
83         CkAbort("Abort");
84       }
85
86       // store the main proxy
87       mainProxy = thisProxy;
88         
89       if(m->argc == 3) {
90         arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
91         blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]); 
92       }
93       else if (m->argc == 7) {
94         arrayDimX = atoi(m->argv[1]);
95         arrayDimY = atoi(m->argv[2]);
96         arrayDimZ = atoi(m->argv[3]);
97         blockDimX = atoi(m->argv[4]); 
98         blockDimY = atoi(m->argv[5]); 
99         blockDimZ = atoi(m->argv[6]);
100       }
101
102       if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
103         CkAbort("array_size_X % block_size_X != 0!");
104       if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
105         CkAbort("array_size_Y % block_size_Y != 0!");
106       if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
107         CkAbort("array_size_Z % block_size_Z != 0!");
108
109       num_chare_x = arrayDimX / blockDimX;
110       num_chare_y = arrayDimY / blockDimY;
111       num_chare_z = arrayDimZ / blockDimZ;
112
113       // print info
114       CkPrintf("\nSTENCIL COMPUTATION WITH BARRIERS\n");
115       CkPrintf("Running Stencil on %d processors with (%d, %d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y, num_chare_z);
116       CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
117       CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
118
119       // Create new array of worker chares
120       array = CProxy_Stencil::ckNew(num_chare_x, num_chare_y, num_chare_z);
121
122       //Start the computation
123       startTime = CmiWallTimer();
124       array.doStep();
125     }
126
127     // Each worker reports back to here when it completes an iteration
128     void report() {
129       CkExit();
130     }
131 };
132
133 /** \class Stencil
134  *
135  */
136
137 class Stencil: public CBase_Stencil {
138   Stencil_SDAG_CODE
139
140   public:
141     int iterations;
142     int imsg;
143
144     double *temperature;
145     double *new_temperature;
146
147     // Constructor, initialize values
148     Stencil() {
149       __sdag_init();
150       usesAtSync=CmiTrue;
151
152       int i, j, k;
153       // allocate a three dimensional array
154       temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
155       new_temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
156
157       for(k=0; k<blockDimZ+2; ++k)
158         for(j=0; j<blockDimY+2; ++j)
159           for(i=0; i<blockDimX+2; ++i)
160             temperature[index(i, j, k)] = 0.0;
161
162       iterations = 0;
163       imsg = 0;
164       constrainBC();
165     }
166
167     void pup(PUP::er &p)
168     {
169       CBase_Stencil::pup(p);
170       __sdag_pup(p);
171       p|iterations;
172       p|imsg;
173
174       size_t size = (blockDimX+2) * (blockDimY+2) * (blockDimZ+2);
175       if (p.isUnpacking()) {
176         temperature = new double[size];
177         new_temperature = new double[size];
178       }
179       p(temperature, size);
180       p(new_temperature, size);
181     }
182
183     Stencil(CkMigrateMessage* m) { __sdag_init(); }
184
185     ~Stencil() { 
186       delete [] temperature; 
187       delete [] new_temperature; 
188     }
189
190     // Send ghost faces to the six neighbors
191     void begin_iteration(void) {
192       iterations++;
193
194       // Copy different faces into messages
195       double *leftGhost =  new double[blockDimY*blockDimZ];
196       double *rightGhost =  new double[blockDimY*blockDimZ];
197       double *topGhost =  new double[blockDimX*blockDimZ];
198       double *bottomGhost =  new double[blockDimX*blockDimZ];
199       double *frontGhost =  new double[blockDimX*blockDimY];
200       double *backGhost =  new double[blockDimX*blockDimY];
201
202       for(int k=0; k<blockDimZ; ++k)
203         for(int j=0; j<blockDimY; ++j) {
204           leftGhost[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
205           rightGhost[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
206         }
207
208       for(int k=0; k<blockDimZ; ++k)
209         for(int i=0; i<blockDimX; ++i) {
210           topGhost[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
211           bottomGhost[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
212         }
213
214       for(int j=0; j<blockDimY; ++j)
215         for(int i=0; i<blockDimX; ++i) {
216           frontGhost[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
217           backGhost[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
218         }
219
220       // Send my left face
221       thisProxy(wrap_x(thisIndex.x-1), thisIndex.y, thisIndex.z)
222           .receiveGhosts(iterations, RIGHT, blockDimY, blockDimZ, leftGhost);
223       // Send my right face
224       thisProxy(wrap_x(thisIndex.x+1), thisIndex.y, thisIndex.z)
225           .receiveGhosts(iterations, LEFT, blockDimY, blockDimZ, rightGhost);
226       // Send my bottom face
227       thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z)
228           .receiveGhosts(iterations, TOP, blockDimX, blockDimZ, bottomGhost);
229       // Send my top face
230       thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z)
231           .receiveGhosts(iterations, BOTTOM, blockDimX, blockDimZ, topGhost);
232       // Send my front face
233       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z-1))
234           .receiveGhosts(iterations, BACK, blockDimX, blockDimY, frontGhost);
235       // Send my back face
236       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z+1))
237           .receiveGhosts(iterations, FRONT, blockDimX, blockDimY, backGhost);
238
239       delete[] leftGhost;
240       delete[] rightGhost;
241       delete[] bottomGhost;
242       delete[] topGhost;
243       delete[] frontGhost;
244       delete[] backGhost;
245     }
246
247     void processGhosts(int dir, int height, int width, double gh[]) {
248       switch(dir) {
249         case LEFT:
250           for(int k=0; k<width; ++k)
251             for(int j=0; j<height; ++j) {
252               temperature[index(0, j+1, k+1)] = gh[k*height+j];
253             }
254           break;
255         case RIGHT:
256           for(int k=0; k<width; ++k)
257             for(int j=0; j<height; ++j) {
258               temperature[index(blockDimX+1, j+1, k+1)] = gh[k*height+j];
259             }
260           break;
261         case BOTTOM:
262           for(int k=0; k<width; ++k)
263             for(int i=0; i<height; ++i) {
264               temperature[index(i+1, 0, k+1)] = gh[k*height+i];
265             }
266           break;
267         case TOP:
268           for(int k=0; k<width; ++k)
269             for(int i=0; i<height; ++i) {
270               temperature[index(i+1, blockDimY+1, k+1)] = gh[k*height+i];
271             }
272           break;
273         case FRONT:
274           for(int j=0; j<width; ++j)
275             for(int i=0; i<height; ++i) {
276               temperature[index(i+1, j+1, 0)] = gh[j*height+i];
277             }
278           break;
279         case BACK:
280           for(int j=0; j<width; ++j)
281             for(int i=0; i<height; ++i) {
282               temperature[index(i+1, j+1, blockDimZ+1)] = gh[j*height+i];
283             }
284           break;
285         default:
286           CkAbort("ERROR\n");
287       }
288     }
289
290
291     void check_and_compute() {
292       compute_kernel();
293
294       // calculate error
295       // not being done right now since we are doing a fixed no. of iterations
296
297       double *tmp;
298       tmp = temperature;
299       temperature = new_temperature;
300       new_temperature = tmp;
301
302       constrainBC();
303
304       if(thisIndex.x == 0 && thisIndex.y == 0 && thisIndex.z == 0) {
305         endTime = CmiWallTimer();
306         CkPrintf("[%d] Time per iteration: %f %f\n", iterations, (endTime - startTime), endTime);
307       }
308
309       if(iterations == MAX_ITER)
310         contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
311       else {
312         startTime = CmiWallTimer();
313         if(iterations % LBPERIOD == 0)
314           AtSync();
315         else
316           contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Stencil::doStep(), thisProxy));
317       }
318     }
319
320     // Check to see if we have received all neighbor values yet
321     // If all neighbor values have been received, we update our values and proceed
322     void compute_kernel() {
323 #pragma unroll    
324       for(int k=1; k<blockDimZ+1; ++k)
325         for(int j=1; j<blockDimY+1; ++j)
326           for(int i=1; i<blockDimX+1; ++i) {
327             // update my value based on the surrounding values
328             new_temperature[index(i, j, k)] = (temperature[index(i-1, j, k)] 
329                                             +  temperature[index(i+1, j, k)]
330                                             +  temperature[index(i, j-1, k)]
331                                             +  temperature[index(i, j+1, k)]
332                                             +  temperature[index(i, j, k-1)]
333                                             +  temperature[index(i, j, k+1)]
334                                             +  temperature[index(i, j, k)] ) * DIVIDEBY7;
335           } // end for
336     }
337
338     // Enforce some boundary conditions
339     void constrainBC() {
340       // Heat left, top and front faces of each chare's block
341       for(int k=1; k<blockDimZ+1; ++k)
342         for(int i=1; i<blockDimX+1; ++i)
343           temperature[index(i, 1, k)] = 255.0;
344       for(int k=1; k<blockDimZ+1; ++k)
345         for(int j=1; j<blockDimY+1; ++j)
346           temperature[index(1, j, k)] = 255.0;
347       for(int j=1; j<blockDimY+1; ++j)
348         for(int i=1; i<blockDimX+1; ++i)
349           temperature[index(i, j, 1)] = 255.0;
350     }
351
352     void ResumeFromSync() {
353       doStep();
354     }
355 };
356
357 #include "stencil3d.def.h"