Merging
[charm.git] / examples / charm++ / jacobi1d / jacobi1d.C
1 /** \file jacobi1d.C
2  *  Author: Abhinav S Bhatele
3  *  Date Created: July 16th, 2009
4  *
5  *
6  *    ***********  ^
7  *    *         *  |
8  *    *         *  |
9  *    *         *  X
10  *    *         *  |
11  *    *         *  |
12  *    ***********  ~
13  *    <--- Y --->
14  *
15  *    X: blockDimX, arrayDimX --> wrap_x
16  *    Y: arrayDimY --> wrap_y
17  * 
18  *  This example does a 1D decomposition of a 2D data array
19  */
20
21 #include "jacobi1d.decl.h"
22 #include "TopoManager.h"
23
24 /*readonly*/ CProxy_Main mainProxy;
25 /*readonly*/ int blockDimX;
26 /*readonly*/ int arrayDimX;
27 /*readonly*/ int arrayDimY;
28
29 // specify the number of worker chares in each dimension
30 /*readonly*/ int num_chares;
31
32 // We want to wrap entries around, and because mod operator % 
33 // sometimes misbehaves on negative values. -1 maps to the highest value.
34 #define wrap_y(a)  (((a)+num_chares)%num_chares)
35 #define index(a,b) ((b) + (a)*(arrayDimY))
36
37 #define MAX_ITER        100
38 #define WARM_ITER       5
39 #define TOP             1
40 #define BOTTOM          2
41 #define LBPERIOD 5
42
43 double startTime;
44 double endTime;
45
46 class Main : public CBase_Main
47 {
48   public:
49     CProxy_Jacobi array;
50     int iterations;
51
52     Main(CkArgMsg* m) {
53       if (m->argc != 4) {
54         CkPrintf("%s [array_size_X] [array_size_Y] [num_chares]\n", m->argv[0]);
55         CkAbort("Abort");
56       }
57
58       CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
59       arrayDimX = atoi(m->argv[1]);
60       arrayDimY = atoi(m->argv[2]);
61       num_chares = atoi(m->argv[3]);
62       if (arrayDimX < num_chares || arrayDimX % num_chares != 0)
63         CkAbort("array_size_X % num_chares != 0!");
64       blockDimX = arrayDimX / num_chares;
65
66       // store the main proxy
67       mainProxy = thisProxy;
68
69       // print info
70       CkPrintf("Running Jacobi on %d processors with (%d) elements\n", CkNumPes(), num_chares);
71       CkPrintf("Array Dimensions: %d %d\n", arrayDimX, arrayDimY);
72       CkPrintf("Block Dimensions: %d %d\n", blockDimX, arrayDimY);
73
74       // Create new array of worker chares
75       array = CProxy_Jacobi::ckNew(num_chares);
76
77       //Start the computation
78       iterations = 0;
79       array.begin_iteration();
80     }
81
82     // Each worker reports back to here when it completes an iteration
83     void report(CkReductionMsg *msg) {
84       iterations++;
85       if(iterations == WARM_ITER)
86         startTime = CmiWallTimer();
87       double error = *((double *)msg->getData());
88
89       if(iterations < MAX_ITER) {
90         if (iterations % LBPERIOD == 0) {
91           array.pause_for_lb();
92         } else {
93           CkPrintf("Start of iteration %d\n", iterations);
94           array.begin_iteration();
95         }
96       } else {
97         CkPrintf("Completed %d iterations\n", MAX_ITER-1);
98         endTime = CmiWallTimer();
99         CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(MAX_ITER-1-WARM_ITER));
100         CkExit();
101       }
102     }
103
104 };
105
106 class Jacobi: public CBase_Jacobi {
107   public:
108     int arrived_top;
109     int arrived_bottom;
110
111     double *temperature;
112     double *new_temperature;
113     void *sendLogs[4];
114     void *ackLogs[5];
115     int work;
116
117     // Constructor, initialize values
118     Jacobi() {
119       usesAtSync=CmiTrue;
120       int i,j;
121       // allocate two dimensional arrays
122       temperature = new double[(blockDimX+2) * arrayDimY];
123       new_temperature = new double[(blockDimX+2) * arrayDimY];
124      // for (i=0; i<blockDimX+2; i++) {
125      //   temperature[i] = new double[arrayDimY];
126      //   new_temperature[i] = new double[arrayDimY];
127      // }
128       for(i=0;i<blockDimX+2; i++) {
129         for(j=0;j<arrayDimY; j++) {
130           temperature[index(i,j)] = 0.5;
131           new_temperature[index(i,j)] = 0.5;
132         }
133       }
134
135       arrived_top = 0;
136       arrived_bottom = 0;
137
138       work = thisIndex/num_chares * 10 + 1;
139       work = 1;
140       constrainBC();
141     }
142
143     Jacobi(CkMigrateMessage* m) {}
144
145     ~Jacobi() { 
146      // for (int i=0; i<blockDimX+2; i++) {
147      //   delete [] temperature[i];
148      //   delete [] new_temperature[i];
149      // }
150       delete [] temperature; 
151       delete [] new_temperature; 
152     }
153
154     void pup(PUP::er &p) {
155       CBase_Jacobi::pup(p);
156       if (p.isUnpacking()) {
157         temperature = new double[(blockDimX+2) * arrayDimY];
158         new_temperature = new double[(blockDimX+2) * arrayDimY];
159       }
160       p|arrived_top;
161       p|arrived_bottom;
162       p(temperature, (blockDimX+2) * arrayDimY);
163       p(new_temperature, (blockDimX+2) * arrayDimY);
164
165       p|work;
166     }
167
168     void ResumeFromSync() {
169       double max_error = 0.0;
170       contribute(sizeof(double), &max_error, CkReduction::max_double,
171           CkCallback(CkIndex_Main::report(NULL), mainProxy));
172     }
173
174     // Perform one iteration of work
175     void begin_iteration(void) {
176       // Send my top edge
177     //  thisProxy(wrap_y(thisIndex)).receiveGhosts(BOTTOM, arrayDimY, &temperature[index(1, 1)]);
178     //  // Send my bottom edge
179     //  thisProxy(wrap_y(thisIndex)).receiveGhosts(TOP, arrayDimY, &temperature[index(blockDimX,1)]);
180     check_and_compute();
181     }
182
183     void receiveGhosts(int dir, int size, double gh[]) {
184       int i, j;
185
186       switch(dir) {
187         case TOP:
188           arrived_top++;
189           for(j=0; j<(size-1); j++)
190             temperature[index(0,j+1)] = gh[j];
191           break;
192         case BOTTOM:
193           arrived_bottom++;
194           for(j=0; j<(size-1); j++)
195             temperature[index(blockDimX+1,j+1)] = gh[j];
196           break;
197         default:
198           CkAbort("ERROR\n");
199       }
200       check_and_compute();
201     }
202
203     void pause_for_lb() {
204       AtSync();
205     }
206
207     void check_and_compute() {
208       double error = 0.0, max_error = 0.0;
209       arrived_top = 1; arrived_bottom = 1;
210
211       if (arrived_top >=1 && arrived_bottom >= 1) {
212         arrived_top--;
213         arrived_bottom--;
214
215         compute_kernel();       
216
217         for(int k = 0; k< work; k++) {
218           for(int i=1; i<blockDimX+1; i++) {
219             for(int j=0; j<arrayDimY; j++) {
220               error = fabs(new_temperature[index(i,j)] - temperature[index(i,j)]);
221               if(error > max_error) {
222                 max_error = error;
223               }
224             }
225           }
226         }
227
228         double *tmp;
229         tmp = temperature;
230         temperature = new_temperature;
231         new_temperature = tmp;
232
233         constrainBC();
234
235         contribute(sizeof(double), &max_error, CkReduction::max_double,
236             CkCallback(CkIndex_Main::report(NULL), mainProxy));
237       }
238     } 
239     // Check to see if we have received all neighbor values yet
240     // If all neighbor values have been received, we update our values and proceed
241     void compute_kernel()
242     {
243       for(int i=1; i<blockDimX+1; i++) {
244         for(int j=0; j<arrayDimY; j++) {
245           // update my value based on the surrounding values
246           new_temperature[index(i,j)] =
247           (temperature[index(i-1,j)]+temperature[index(i+1,j)]+temperature[index(i,j-1)]+temperature[index(i,j+1)]+temperature[index(i,j)]) * 0.2;
248         }
249       }
250     }
251
252     // Enforce some boundary conditions
253     void constrainBC()
254     {
255       if(thisIndex <= num_chares/2) {
256         for(int i=1; i<=blockDimX; i++)
257           temperature[index(i,1)] = 1.0;
258       }
259
260       if(thisIndex == num_chares-1) {
261         for(int j=arrayDimY/2; j<arrayDimY; j++)
262           temperature[index(blockDimX,j)] = 0.0;
263       }
264     }
265
266 };
267
268 #include "jacobi1d.def.h"