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