b43ce61757e0612b1b67521de6abef7e7ac18694
[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        26
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       constrainBC();
138     }
139
140     Jacobi(CkMigrateMessage* m) {}
141
142     ~Jacobi() { 
143      // for (int i=0; i<blockDimX+2; i++) {
144      //   delete [] temperature[i];
145      //   delete [] new_temperature[i];
146      // }
147       delete [] temperature; 
148       delete [] new_temperature; 
149     }
150
151     void pup(PUP::er &p) {
152       CBase_Jacobi::pup(p);
153       if (p.isUnpacking()) {
154         temperature = new double[(blockDimX+2) * arrayDimY];
155         new_temperature = new double[(blockDimX+2) * arrayDimY];
156       }
157       p(temperature, (blockDimX+2) * arrayDimY);
158       p(new_temperature, (blockDimX+2) * arrayDimY);
159
160       p|iterations;
161       p|work;
162     }
163
164     void ResumeFromSync() {
165       double max_error = 0.0;
166       contribute(sizeof(double), &max_error, CkReduction::max_double,
167           CkCallback(CkIndex_Main::report(NULL), mainProxy));
168     }
169
170     // Perform one iteration of work
171     void begin_iteration(void) {
172       if (iterations % LBPERIOD == 0) {
173         AtSync();
174         return;
175       }
176
177       // Send my top edge
178       thisProxy(wrap_y(thisIndex)).receiveGhosts(BOTTOM, arrayDimY, &temperature[index(1, 1)]);
179       // Send my bottom edge
180       thisProxy(wrap_y(thisIndex)).receiveGhosts(TOP, arrayDimY, &temperature[index(blockDimX,1)]);
181       iterations++;
182     }
183
184     void receiveGhosts(int dir, int size, double gh[]) {
185       int i, j;
186
187       switch(dir) {
188         case TOP:
189           arrived_top++;
190           for(j=0; j<size; j++)
191             temperature[index(0,j+1)] = gh[j];
192           break;
193         case BOTTOM:
194           arrived_bottom++;
195           for(j=0; j<size; j++)
196             temperature[index(blockDimX+1,j+1)] = gh[j];
197           break;
198         default:
199           CkAbort("ERROR\n");
200       }
201       check_and_compute();
202     }
203
204     void check_and_compute() {
205       double error = 0.0, max_error = 0.0;
206
207       if (arrived_top >=1 && arrived_bottom >= 1) {
208         arrived_top--;
209         arrived_bottom--;
210
211         compute_kernel();       
212
213         for(int k = 0; k< work; k++) {
214           for(int i=1; i<blockDimX+1; i++) {
215             for(int j=0; j<arrayDimY; j++) {
216               error = fabs(new_temperature[index(i,j)] - temperature[index(i,j)]);
217               if(error > max_error) {
218                 max_error = error;
219               }
220             }
221           }
222         }
223
224         double *tmp;
225         tmp = temperature;
226         temperature = new_temperature;
227         new_temperature = tmp;
228
229         constrainBC();
230
231         contribute(sizeof(double), &max_error, CkReduction::max_double,
232             CkCallback(CkIndex_Main::report(NULL), mainProxy));
233       }
234     }
235
236     // Check to see if we have received all neighbor values yet
237     // If all neighbor values have been received, we update our values and proceed
238     void compute_kernel()
239     {
240       for(int i=1; i<blockDimX+1; i++) {
241         for(int j=0; j<arrayDimY; j++) {
242           // update my value based on the surrounding values
243           new_temperature[index(i,j)] =
244           (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;
245         }
246       }
247     }
248
249     // Enforce some boundary conditions
250     void constrainBC()
251     {
252       if(thisIndex <= num_chares/2) {
253         for(int i=1; i<=blockDimX; i++)
254           temperature[index(i,1)] = 1.0;
255       }
256
257       if(thisIndex == num_chares-1) {
258         for(int j=arrayDimY/2; j<arrayDimY; j++)
259           temperature[index(blockDimX,j)] = 0.0;
260       }
261     }
262
263 };
264
265 #include "jacobi1d.def.h"