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