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