Enable node-queue when using multicore build
[charm.git] / examples / charm++ / jacobi3d-sdag / jacobi3d.C
1 #define MAX_ITER                26
2 #define LEFT                    1
3 #define RIGHT                   2
4 #define TOP                     3
5 #define BOTTOM                  4
6 #define FRONT                   5
7 #define BACK                    6
8 #define DIVIDEBY7               0.14285714285714285714
9 #define DELTA                   0.01
10 #define LBPERIOD                50
11 #define CHECKPOINTPERIOD        500
12
13 #include "jacobi3d.decl.h"
14
15 /*readonly*/ CProxy_Main mainProxy;
16 /*readonly*/ int arrayDimX;
17 /*readonly*/ int arrayDimY;
18 /*readonly*/ int arrayDimZ;
19 /*readonly*/ int blockDimX;
20 /*readonly*/ int blockDimY;
21 /*readonly*/ int blockDimZ;
22
23 // specify the number of worker chares in each dimension
24 /*readonly*/ int num_chare_x;
25 /*readonly*/ int num_chare_y;
26 /*readonly*/ int num_chare_z;
27
28 #define wrapX(a)        (((a)+num_chare_x)%num_chare_x)
29 #define wrapY(a)        (((a)+num_chare_y)%num_chare_y)
30 #define wrapZ(a)        (((a)+num_chare_z)%num_chare_z)
31
32 #define index(a,b,c)    ((a)+(b)*(blockDimX+2)+(c)*(blockDimX+2)*(blockDimY+2))
33
34 double startTime;
35 double endTime;
36
37 /** \class Main
38  *
39  */
40 class Main : public CBase_Main {
41 public:
42   CProxy_Jacobi array;
43   int iterations;
44
45   Main(CkArgMsg* m) {
46     if ( (m->argc != 3) && (m->argc != 7) ) {
47       CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
48       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]);
49       CkAbort("Abort");
50     }
51
52     // set iteration counter to zero
53     iterations = 0;
54
55     // store the main proxy
56     mainProxy = thisProxy;
57         
58     if(m->argc == 3) {
59       arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
60       blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]); 
61     }
62     else if (m->argc == 7) {
63       arrayDimX = atoi(m->argv[1]);
64       arrayDimY = atoi(m->argv[2]);
65       arrayDimZ = atoi(m->argv[3]);
66       blockDimX = atoi(m->argv[4]); 
67       blockDimY = atoi(m->argv[5]); 
68       blockDimZ = atoi(m->argv[6]);
69     }
70
71     if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
72       CkAbort("array_size_X % block_size_X != 0!");
73     if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
74       CkAbort("array_size_Y % block_size_Y != 0!");
75     if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
76       CkAbort("array_size_Z % block_size_Z != 0!");
77
78     num_chare_x = arrayDimX / blockDimX;
79     num_chare_y = arrayDimY / blockDimY;
80     num_chare_z = arrayDimZ / blockDimZ;
81
82     // print info
83     CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
84     CkPrintf("Running Jacobi on %d processors with (%d, %d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y, num_chare_z);
85     CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
86     CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
87
88     // Create new array of worker chares
89     array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y, num_chare_z);
90
91     //Start the computation
92     array.run();
93     startTime = CkWallTimer();
94   }
95
96   void done(int iterations) {
97     CkPrintf("Completed %d iterations\n", iterations);
98     endTime = CkWallTimer();
99     CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime) / iterations);
100     CkExit();
101   }
102 };
103
104 /** \class Jacobi
105  *
106  */
107
108 class Jacobi: public CBase_Jacobi {
109   Jacobi_SDAG_CODE
110
111 public:
112   int iterations;
113   int remoteCount;
114
115   double *temperature;
116   double *new_temperature;
117   bool converged;
118
119   // Constructor, initialize values
120   Jacobi() {
121     usesAtSync = true;
122     converged = false;
123
124     // allocate a three dimensional array
125     temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
126     new_temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
127
128     for(int k=0; k<blockDimZ+2; ++k)
129       for(int j=0; j<blockDimY+2; ++j)
130         for(int i=0; i<blockDimX+2; ++i)
131           temperature[index(i, j, k)] = 0.0;
132
133     iterations = 0;
134     constrainBC();
135   }
136
137   void pup(PUP::er &p)
138   {
139     CBase_Jacobi::pup(p);
140     __sdag_pup(p);
141     p|iterations;
142     p|remoteCount;
143
144     size_t size = (blockDimX+2) * (blockDimY+2) * (blockDimZ+2);
145     if (p.isUnpacking()) {
146       temperature = new double[size];
147       new_temperature = new double[size];
148     }
149     p(temperature, size);
150     p(new_temperature, size);
151   }
152
153   Jacobi(CkMigrateMessage* m) { }
154
155   ~Jacobi() { 
156     delete [] temperature; 
157     delete [] new_temperature; 
158   }
159
160   // Send ghost faces to the six neighbors
161   void begin_iteration(void) {
162     // Copy different faces into messages
163     double *leftGhost =  new double[blockDimY*blockDimZ];
164     double *rightGhost =  new double[blockDimY*blockDimZ];
165     double *topGhost =  new double[blockDimX*blockDimZ];
166     double *bottomGhost =  new double[blockDimX*blockDimZ];
167     double *frontGhost =  new double[blockDimX*blockDimY];
168     double *backGhost =  new double[blockDimX*blockDimY];
169
170     for(int k=0; k<blockDimZ; ++k)
171       for(int j=0; j<blockDimY; ++j) {
172         leftGhost[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
173         rightGhost[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
174       }
175
176     for(int k=0; k<blockDimZ; ++k)
177       for(int i=0; i<blockDimX; ++i) {
178         topGhost[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
179         bottomGhost[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
180       }
181
182     for(int j=0; j<blockDimY; ++j)
183       for(int i=0; i<blockDimX; ++i) {
184         frontGhost[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
185         backGhost[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
186       }
187
188     int x = thisIndex.x, y = thisIndex.y, z = thisIndex.z;
189     thisProxy(wrapX(x-1),y,z).updateGhosts(iterations, RIGHT,  blockDimY, blockDimZ, rightGhost);
190     thisProxy(wrapX(x+1),y,z).updateGhosts(iterations, LEFT,   blockDimY, blockDimZ, leftGhost);
191     thisProxy(x,wrapY(y-1),z).updateGhosts(iterations, TOP,    blockDimX, blockDimZ, topGhost);
192     thisProxy(x,wrapY(y+1),z).updateGhosts(iterations, BOTTOM, blockDimX, blockDimZ, bottomGhost);
193     thisProxy(x,y,wrapZ(z-1)).updateGhosts(iterations, BACK,   blockDimX, blockDimY, backGhost);
194     thisProxy(x,y,wrapZ(z+1)).updateGhosts(iterations, FRONT,  blockDimX, blockDimY, frontGhost);
195
196     delete [] leftGhost;
197     delete [] rightGhost;
198     delete [] bottomGhost;
199     delete [] topGhost;
200     delete [] frontGhost;
201     delete [] backGhost;
202   }
203
204   void updateBoundary(int dir, int height, int width, double* gh) {
205     switch(dir) {
206     case LEFT:
207       for(int k=0; k<width; ++k)
208         for(int j=0; j<height; ++j) {
209           temperature[index(0, j+1, k+1)] = gh[k*height+j];
210         }
211       break;
212     case RIGHT:
213       for(int k=0; k<width; ++k)
214         for(int j=0; j<height; ++j) {
215           temperature[index(blockDimX+1, j+1, k+1)] = gh[k*height+j];
216         }
217       break;
218     case BOTTOM:
219       for(int k=0; k<width; ++k)
220         for(int i=0; i<height; ++i) {
221           temperature[index(i+1, 0, k+1)] = gh[k*height+i];
222         }
223       break;
224     case TOP:
225       for(int k=0; k<width; ++k)
226         for(int i=0; i<height; ++i) {
227           temperature[index(i+1, blockDimY+1, k+1)] = gh[k*height+i];
228         }
229       break;
230     case FRONT:
231       for(int j=0; j<width; ++j)
232         for(int i=0; i<height; ++i) {
233           temperature[index(i+1, j+1, 0)] = gh[j*height+i];
234         }
235       break;
236     case BACK:
237       for(int j=0; j<width; ++j)
238         for(int i=0; i<height; ++i) {
239           temperature[index(i+1, j+1, blockDimZ+1)] = gh[j*height+i];
240         }
241       break;
242     default:
243       CkAbort("ERROR\n");
244     }
245   }
246
247   // Check to see if we have received all neighbor values yet
248   // If all neighbor values have been received, we update our values and proceed
249   double computeKernel() {
250 #pragma unroll    
251     for(int k=1; k<blockDimZ+1; ++k)
252       for(int j=1; j<blockDimY+1; ++j)
253         for(int i=1; i<blockDimX+1; ++i) {
254           // update my value based on the surrounding values
255           new_temperature[index(i, j, k)] = (temperature[index(i-1, j, k)] 
256                                              +  temperature[index(i+1, j, k)]
257                                              +  temperature[index(i, j-1, k)]
258                                              +  temperature[index(i, j+1, k)]
259                                              +  temperature[index(i, j, k-1)]
260                                              +  temperature[index(i, j, k+1)]
261                                              +  temperature[index(i, j, k)] ) * DIVIDEBY7;
262         } // end for
263     
264     double error = 0.0, max_error = 0.0;
265
266     for(int k=1; k<blockDimZ+1; ++k)
267       for(int j=1; j<blockDimY+1; ++j)
268         for(int i=1; i<blockDimX+1; ++i) {
269           error = fabs(new_temperature[index(i,j,k)] - temperature[index(i,j,k)]);
270           if (error > max_error) {
271             max_error = error;
272           }
273         }
274
275     double *tmp;
276     tmp = temperature;
277     temperature = new_temperature;
278     new_temperature = tmp;
279
280     constrainBC();
281
282     return max_error;
283   }
284
285   // Enforce some boundary conditions
286   void constrainBC() {
287     // Heat left, top and front faces of each chare's block
288     for(int k=1; k<blockDimZ+1; ++k)
289       for(int i=1; i<blockDimX+1; ++i)
290         temperature[index(i, 1, k)] = 255.0;
291     for(int k=1; k<blockDimZ+1; ++k)
292       for(int j=1; j<blockDimY+1; ++j)
293         temperature[index(1, j, k)] = 255.0;
294     for(int j=1; j<blockDimY+1; ++j)
295       for(int i=1; i<blockDimX+1; ++i)
296         temperature[index(i, j, 1)] = 255.0;
297   }
298
299 };
300
301 #include "jacobi3d.def.h"