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