communication aware lb strategy
[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) ((a) + (b)*(blockDimX+2))
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         CkPrintf("Start of iteration %d\n", iterations);
91         array.begin_iteration();
92       } else {
93         CkPrintf("Completed %d iterations\n", MAX_ITER-1);
94         endTime = CmiWallTimer();
95         CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(MAX_ITER-1-WARM_ITER));
96         CkExit();
97       }
98     }
99
100 };
101
102 class Jacobi: public CBase_Jacobi {
103   public:
104     int arrived_top;
105     int arrived_bottom;
106
107     double *temperature;
108     double *new_temperature;
109     void *sendLogs[4];
110     void *ackLogs[5];
111     int iterations;
112     int work;
113
114     // Constructor, initialize values
115     Jacobi() {
116       usesAtSync=CmiTrue;
117       int i,j;
118       // allocate two dimensional arrays
119       temperature = new double[(blockDimX+2) * arrayDimY];
120       new_temperature = new double[(blockDimX+2) * arrayDimY];
121      // for (i=0; i<blockDimX+2; i++) {
122      //   temperature[i] = new double[arrayDimY];
123      //   new_temperature[i] = new double[arrayDimY];
124      // }
125       for(i=0;i<blockDimX+2; i++) {
126         for(j=0;j<arrayDimY; j++) {
127           temperature[index(i,j)] = 0.5;
128           new_temperature[index(i,j)] = 0.5;
129         }
130       }
131
132       arrived_top = 0;
133       arrived_bottom = 0;
134       iterations = 0;
135
136       //work = thisIndex;
137       work = 1;
138       constrainBC();
139     }
140
141     Jacobi(CkMigrateMessage* m) {}
142
143     ~Jacobi() { 
144      // for (int i=0; i<blockDimX+2; i++) {
145      //   delete [] temperature[i];
146      //   delete [] new_temperature[i];
147      // }
148       delete [] temperature; 
149       delete [] new_temperature; 
150     }
151
152     void pup(PUP::er &p) {
153       CBase_Jacobi::pup(p);
154       if (p.isUnpacking()) {
155         temperature = new double[(blockDimX+2) * arrayDimY];
156         new_temperature = new double[(blockDimX+2) * arrayDimY];
157       }
158       p(temperature, (blockDimX+2) * arrayDimY);
159       p(new_temperature, (blockDimX+2) * arrayDimY);
160
161       p|iterations;
162       p|work;
163     }
164
165     void ResumeFromSync() {
166       double max_error = 0.0;
167       contribute(sizeof(double), &max_error, CkReduction::max_double,
168           CkCallback(CkIndex_Main::report(NULL), mainProxy));
169     }
170
171     // Perform one iteration of work
172     void begin_iteration(void) {
173       if (iterations % LBPERIOD == 0) {
174         AtSync();
175         return;
176       }
177
178       // Send my top edge
179       thisProxy(wrap_y(thisIndex)).receiveGhosts(BOTTOM, arrayDimY, &temperature[index(1, 1)]);
180       // Send my bottom edge
181       thisProxy(wrap_y(thisIndex)).receiveGhosts(TOP, arrayDimY, &temperature[index(blockDimX,1)]);
182       iterations++;
183     }
184
185     void receiveGhosts(int dir, int size, double gh[]) {
186       int i, j;
187
188       switch(dir) {
189         case TOP:
190           arrived_top++;
191           for(j=0; j<size; j++)
192             temperature[index(0,j+1)] = gh[j];
193           break;
194         case BOTTOM:
195           arrived_bottom++;
196           for(j=0; j<size; j++)
197             temperature[index(blockDimX+1,j+1)] = gh[j];
198           break;
199         default:
200           CkAbort("ERROR\n");
201       }
202       check_and_compute();
203     }
204
205     void check_and_compute() {
206       double error = 0.0, max_error = 0.0;
207
208       if (arrived_top >=1 && arrived_bottom >= 1) {
209         arrived_top--;
210         arrived_bottom--;
211
212         compute_kernel();       
213
214         for(int k = 0; k< work; k++) {
215           for(int i=1; i<blockDimX+1; i++) {
216             for(int j=0; j<arrayDimY; j++) {
217               error = fabs(new_temperature[index(i,j)] - temperature[index(i,j)]);
218               if(error > max_error) {
219                 max_error = error;
220               }
221             }
222           }
223         }
224
225         double *tmp;
226         tmp = temperature;
227         temperature = new_temperature;
228         new_temperature = tmp;
229
230         constrainBC();
231
232         contribute(sizeof(double), &max_error, CkReduction::max_double,
233             CkCallback(CkIndex_Main::report(NULL), mainProxy));
234       }
235     }
236
237     // Check to see if we have received all neighbor values yet
238     // If all neighbor values have been received, we update our values and proceed
239     void compute_kernel()
240     {
241       for(int i=1; i<blockDimX+1; i++) {
242         for(int j=0; j<arrayDimY; j++) {
243           // update my value based on the surrounding values
244           new_temperature[index(i,j)] =
245           (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;
246         }
247       }
248     }
249
250     // Enforce some boundary conditions
251     void constrainBC()
252     {
253       if(thisIndex <= num_chares/2) {
254         for(int i=1; i<=blockDimX; i++)
255           temperature[index(i,1)] = 1.0;
256       }
257
258       if(thisIndex == num_chares-1) {
259         for(int j=arrayDimY/2; j<arrayDimY; j++)
260           temperature[index(blockDimX,j)] = 0.0;
261       }
262     }
263
264 };
265
266 #include "jacobi1d.def.h"