examples/charm++: Always use CkWallTimer, not CmiWallTimer (tested)
[charm.git] / examples / charm++ / gaussSeidel3D / gaussSeidel3d.C
1
2 #include "gaussSeidel3d.decl.h"
3 #include "TopoManager.h"
4
5 // See README for documentation
6
7 /*readonly*/ CProxy_Main mainProxy;
8 /*readonly*/ int arrayDimX;
9 /*readonly*/ int arrayDimY;
10 /*readonly*/ int arrayDimZ;
11 /*readonly*/ int blockDimX;
12 /*readonly*/ int blockDimY;
13 /*readonly*/ int blockDimZ;
14
15 // specify the number of worker chares in each dimension
16 /*readonly*/ int num_chare_x;
17 /*readonly*/ int num_chare_y;
18 /*readonly*/ int num_chare_z;
19
20 /*readonly*/ int globalBarrier;
21
22 static unsigned long next = 1;
23
24 int myrand(int numpes) {
25   next = next * 1103515245 + 12345;
26   return((unsigned)(next/65536) % numpes);
27 }
28
29
30 #define USE_3D_ARRAYS           0
31 #if USE_3D_ARRAYS
32 #define index(a, b, c)  a][b][c 
33 #else
34 #define index(a, b, c)  (a*(blockDimY+2)*(blockDimZ+2) + b*(blockDimZ+2) + c)
35 #endif
36
37 #define MAX_ITER                10
38 #define LEFT                    1
39 #define RIGHT                   2
40 #define TOP                     3
41 #define BOTTOM                  4
42 #define FRONT                   5
43 #define BACK                    6
44 #define DIVIDEBY7               0.14285714285714285714
45
46
47 char * dirstring(int dir){
48   switch(dir){
49   case LEFT:
50     return "LEFT";
51   case RIGHT:
52     return "RIGHT";
53   case TOP:
54     return "TOP";
55   case BOTTOM:
56     return "BOTTOM";
57   case FRONT:
58     return "FRONT";
59   case BACK:
60     return "BACK";
61   }
62 }
63
64
65 double startTime;
66 double endTime;
67
68 /** \class ghostMsg
69  *
70  */
71 class ghostMsg: public CMessage_ghostMsg {
72 public:
73   int dir;
74   int height;
75   int width;
76   double* gh;
77
78   ghostMsg(int _d, int _h, int _w) : dir(_d), height(_h), width(_w) {
79   }
80 };
81
82 /** \class Main
83  *
84  */
85 class Main : public CBase_Main {
86 public:
87   CProxy_GaussSeidel array;
88   int iterations;
89
90   Main(CkArgMsg* m) {
91     if ( (m->argc != 3) && (m->argc != 7) ) {
92       CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
93       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]);
94       CkAbort("Abort");
95     }
96
97     traceRegisterUserEvent("Begin Iteration ***", 1000);
98
99     // set iteration counter to zero
100     iterations = 0;
101       
102     // store the main proxy
103     mainProxy = thisProxy;
104       
105     if(m->argc == 3) {
106       arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
107       blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]); 
108     }
109     else if (m->argc == 7) {
110       arrayDimX = atoi(m->argv[1]);
111       arrayDimY = atoi(m->argv[2]);
112       arrayDimZ = atoi(m->argv[3]);
113       blockDimX = atoi(m->argv[4]); 
114       blockDimY = atoi(m->argv[5]); 
115       blockDimZ = atoi(m->argv[6]);
116     }
117     
118     if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
119       CkAbort("array_size_X % block_size_X != 0!");
120     if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
121       CkAbort("array_size_Y % block_size_Y != 0!");
122     if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
123       CkAbort("array_size_Z % block_size_Z != 0!");
124     
125     num_chare_x = arrayDimX / blockDimX;
126     num_chare_y = arrayDimY / blockDimY;
127     num_chare_z = arrayDimZ / blockDimZ;
128       
129     // print info
130     CkPrintf("\nSTENCIL Gauss Seidel COMPUTATION WITH NO BARRIERS\n");
131     CkPrintf("Running GaussSeidel on %d processors with (%d, %d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y, num_chare_z);
132     CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
133     CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
134
135     // Create new array of worker chares
136     array = CProxy_GaussSeidel::ckNew(num_chare_x, num_chare_y, num_chare_z); 
137
138     startTime = CkWallTimer();
139
140     //Start the computation
141     array.doStep();
142   }
143
144   // Each worker reports back to here when it completes an iteration
145   void report() {
146     iterations++;
147     endTime = CkWallTimer();
148     CkPrintf("Average time for each of the first %d iteration: %f\n", iterations ,(endTime - startTime)/(iterations));
149     if(iterations >= MAX_ITER){
150       CkExit();
151     } else {
152       array.doStep();     
153     }
154   }
155 };
156
157 /** \class GaussSeidel
158  *
159  */
160
161 class GaussSeidel: public CBase_GaussSeidel {
162   GaussSeidel_SDAG_CODE
163
164   public:
165   int iterations;
166   int imsg;
167
168 #if USE_3D_ARRAYS
169   double ***temperature;
170 #else
171   double *temperature;
172 #endif
173
174   // Constructor, initialize values
175   GaussSeidel() {
176     __sdag_init();
177
178     int i, j, k;
179     // allocate a three dimensional array
180 #if USE_3D_ARRAYS
181     temperature = new double**[blockDimX+2];
182     for (i=0; i<blockDimX+2; i++) {
183       temperature[i] = new double*[blockDimY+2];
184       for(j=0; j<blockDimY+2; j++) {
185         temperature[i][j] = new double[blockDimZ+2];
186       }
187     }
188 #else
189     temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
190 #endif
191
192     for(i=0; i<blockDimX+2; ++i) {
193       for(j=0; j<blockDimY+2; ++j) {
194         for(k=0; k<blockDimZ+2; ++k) {
195           temperature[index(i, j, k)] = 0.0;
196         }
197       } 
198     }
199
200     iterations = 0;
201     imsg = 0;
202     constrainBC();
203   }
204
205   GaussSeidel(CkMigrateMessage* m) {}
206
207   ~GaussSeidel() { 
208 #if USE_3D_ARRAYS
209     for (int i=0; i<blockDimX+2; i++) {
210       for(int j=0; j<blockDimY+2; j++) {
211         delete [] temperature[i][j];
212       }
213       delete [] temperature[i];
214     }
215     delete [] temperature; 
216 #else
217     delete [] temperature; 
218 #endif
219   }
220
221   // Send ghost faces to three of the neighbors containing the values computed during the previous step
222   void begin_iteration(void) {
223     //    CkPrintf("Elem %d,%d,%d Start of iteration %d\n", thisIndex.x, thisIndex.y, thisIndex.z, iterations);
224     if(thisIndex.x == 0 and thisIndex.y == 0 and thisIndex.z == 0)
225       traceUserEvent(1000);
226
227     // Copy different faces into messages
228       
229       
230     // Send my left face
231     if(thisIndex.x > 0){
232       ghostMsg *leftMsg = new (blockDimY*blockDimZ) ghostMsg(RIGHT, blockDimY, blockDimZ);
233       CkSetRefNum(leftMsg, iterations);
234       for(int j=0; j<blockDimY; ++j) 
235         for(int k=0; k<blockDimZ; ++k) 
236           leftMsg->gh[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
237       thisProxy(thisIndex.x-1, thisIndex.y, thisIndex.z).receiveGhostsPrevX(leftMsg);
238     }
239       
240     // Send my top face
241     if(thisIndex.y > 0){
242       ghostMsg *topMsg = new (blockDimX*blockDimZ) ghostMsg(BOTTOM, blockDimX, blockDimZ);
243       CkSetRefNum(topMsg, iterations);
244       for(int i=0; i<blockDimX; ++i) 
245         for(int k=0; k<blockDimZ; ++k)
246           topMsg->gh[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
247       thisProxy(thisIndex.x, thisIndex.y-1, thisIndex.z).receiveGhostsPrevY(topMsg);
248     }
249       
250     // Send my front face
251     if(thisIndex.z > 0){
252       ghostMsg *frontMsg = new (blockDimX*blockDimY) ghostMsg(BACK, blockDimX, blockDimY);
253       CkSetRefNum(frontMsg, iterations);
254       for(int i=0; i<blockDimX; ++i) 
255         for(int j=0; j<blockDimY; ++j) 
256           frontMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
257       thisProxy(thisIndex.x, thisIndex.y, thisIndex.z-1).receiveGhostsPrevZ(frontMsg);
258     }
259       
260       
261   }
262   
263   
264   
265   void processGhosts(ghostMsg *gmsg) {
266     //    CkPrintf("Elem %d,%d,%d processGhosts dir=%s\n", thisIndex.x, thisIndex.y, thisIndex.z, dirstring(gmsg->dir) );
267
268     int height = gmsg->height;
269     int width = gmsg->width;
270
271     switch(gmsg->dir) {
272     case LEFT:
273       for(int j=0; j<height; ++j) 
274         for(int k=0; k<width; ++k) {
275           temperature[index(0, j+1, k+1)] = gmsg->gh[k*height+j];
276         }
277       break;
278     case RIGHT:
279       for(int j=0; j<height; ++j) 
280         for(int k=0; k<width; ++k) {
281           temperature[index(blockDimX+1, j+1, k+1)] = gmsg->gh[k*height+j];
282         }
283       break;
284     case TOP:
285       for(int i=0; i<height; ++i) 
286         for(int k=0; k<width; ++k) {
287           temperature[index(i+1, 0, k+1)] = gmsg->gh[k*height+i];
288         }
289       break;
290     case BOTTOM:
291       for(int i=0; i<height; ++i) 
292         for(int k=0; k<width; ++k) {
293           temperature[index(i+1, blockDimY+1, k+1)] = gmsg->gh[k*height+i];
294         }
295       break;
296     case FRONT:
297       for(int i=0; i<height; ++i) 
298         for(int j=0; j<width; ++j) {
299           temperature[index(i+1, j+1, blockDimZ+1)] = gmsg->gh[j*height+i];
300         }
301       break;
302     case BACK:
303       for(int i=0; i<height; ++i) 
304         for(int j=0; j<width; ++j) {
305           temperature[index(i+1, j+1, 0)] = gmsg->gh[j*height+i];
306         }
307       break;
308     default:
309       CkAbort("ERROR\n");
310     }
311
312     delete gmsg;
313   }
314
315
316   void compute() {
317     compute_kernel();
318     constrainBC();
319
320
321     // Send my right face
322     if(thisIndex.x < num_chare_x-1){
323       ghostMsg *rightMsg = new (blockDimY*blockDimZ) ghostMsg(LEFT, blockDimY, blockDimZ);
324       CkSetRefNum(rightMsg, iterations);
325       for(int j=0; j<blockDimY; ++j) 
326         for(int k=0; k<blockDimZ; ++k) {
327           rightMsg->gh[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
328         }
329       thisProxy(thisIndex.x+1, thisIndex.y, thisIndex.z).receiveGhostsCurrentX(rightMsg);
330       //CkPrintf("Elem %d,%d,%d AFTER KERNEL sending to x+1 iter=%d\n", thisIndex.x, thisIndex.y, thisIndex.z, iterations);
331     }
332             
333     // Send my bottom face
334     if(thisIndex.y < num_chare_y-1){
335       ghostMsg *bottomMsg = new (blockDimX*blockDimZ) ghostMsg(TOP, blockDimX, blockDimZ);
336       CkSetRefNum(bottomMsg, iterations);
337       for(int i=0; i<blockDimX; ++i) 
338         for(int k=0; k<blockDimZ; ++k) {
339           bottomMsg->gh[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
340         }
341       thisProxy(thisIndex.x, thisIndex.y+1, thisIndex.z).receiveGhostsCurrentY(bottomMsg);
342       //CkPrintf("Elem %d,%d,%d AFTER KERNEL sending to y+1 iter=%d\n", thisIndex.x, thisIndex.y, thisIndex.z, iterations);
343     }
344       
345     // Send my back face
346     if(thisIndex.z < num_chare_z-1){
347       ghostMsg *backMsg = new (blockDimX*blockDimY) ghostMsg(FRONT, blockDimX, blockDimY);
348       CkSetRefNum(backMsg, iterations);
349       for(int i=0; i<blockDimX; ++i) 
350         for(int j=0; j<blockDimY; ++j) {
351           backMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
352         }
353       thisProxy(thisIndex.x, thisIndex.y, thisIndex.z+1).receiveGhostsCurrentZ(backMsg);
354       //CkPrintf("Elem %d,%d,%d AFTER KERNEL sending to z+1 iter=%d\n", thisIndex.x, thisIndex.y, thisIndex.z, iterations);
355     }
356
357
358     iterations++;
359
360     contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
361   }
362
363
364   void compute_kernel() {
365 #pragma unroll    
366     for(int i=1; i<blockDimX+1; ++i) {
367       for(int j=1; j<blockDimY+1; ++j) {
368         for(int k=1; k<blockDimZ+1; ++k) {
369           // update my value based on the surrounding values
370           temperature[index(i, j, k)] = (temperature[index(i-1, j, k)] 
371                                          +  temperature[index(i+1, j, k)]
372                                          +  temperature[index(i, j-1, k)]
373                                          +  temperature[index(i, j+1, k)]
374                                          +  temperature[index(i, j, k-1)]
375                                          +  temperature[index(i, j, k+1)]
376                                          +  temperature[index(i, j, k)] ) * DIVIDEBY7;
377         }
378       }
379     }
380   }
381
382   // Enforce some boundary conditions
383   void constrainBC() {
384     // Heat left, top and front faces of each chare's block
385     for(int i=1; i<blockDimX+1; ++i)
386       for(int k=1; k<blockDimZ+1; ++k)
387         temperature[index(i, 1, k)] = 255.0;
388     for(int j=1; j<blockDimY+1; ++j)
389       for(int k=1; k<blockDimZ+1; ++k)
390         temperature[index(1, j, k)] = 255.0;
391     for(int i=1; i<blockDimX+1; ++i)
392       for(int j=1; j<blockDimY+1; ++j)
393         temperature[index(i, j, 1)] = 255.0;
394   }
395
396 };
397
398
399 #include "gaussSeidel3d.def.h"