Merge branch 'charm' of charmgit:charm into charm
[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
36 #define MAX_ITER        26
37 #define WARM_ITER       5
38 #define TOP             1
39 #define BOTTOM          2
40
41 double startTime;
42 double endTime;
43
44 class Main : public CBase_Main
45 {
46   public:
47     CProxy_Jacobi array;
48     int iterations;
49
50     Main(CkArgMsg* m) {
51       if (m->argc != 4) {
52         CkPrintf("%s [array_size_X] [array_size_Y] [num_chares]\n", m->argv[0]);
53         CkAbort("Abort");
54       }
55
56       CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
57       arrayDimX = atoi(m->argv[1]);
58       arrayDimY = atoi(m->argv[2]);
59       num_chares = atoi(m->argv[3]);
60       if (arrayDimX < num_chares || arrayDimX % num_chares != 0)
61         CkAbort("array_size_X % num_chares != 0!");
62       blockDimX = arrayDimX / num_chares;
63
64       // store the main proxy
65       mainProxy = thisProxy;
66
67       // print info
68       CkPrintf("Running Jacobi on %d processors with (%d) elements\n", CkNumPes(), num_chares);
69       CkPrintf("Array Dimensions: %d %d\n", arrayDimX, arrayDimY);
70       CkPrintf("Block Dimensions: %d %d\n", blockDimX, arrayDimY);
71
72       // Create new array of worker chares
73       array = CProxy_Jacobi::ckNew(num_chares);
74
75       //Start the computation
76       iterations = 0;
77       array.begin_iteration();
78     }
79
80     // Each worker reports back to here when it completes an iteration
81     void report(CkReductionMsg *msg) {
82       iterations++;
83       if(iterations == WARM_ITER)
84         startTime = CmiWallTimer();
85       double error = *((double *)msg->getData());
86
87       if(iterations < MAX_ITER) {
88         CkPrintf("Start of iteration %d\n", iterations);
89         array.begin_iteration();
90       } else {
91         CkPrintf("Completed %d iterations\n", MAX_ITER-1);
92         endTime = CmiWallTimer();
93         CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(MAX_ITER-1-WARM_ITER));
94         CkExit();
95       }
96     }
97
98 };
99
100 class Jacobi: public CBase_Jacobi {
101   public:
102     int arrived_top;
103     int arrived_bottom;
104
105     double **temperature;
106     double **new_temperature;
107     void *sendLogs[4];
108     void *ackLogs[5];
109     int iterations;
110
111     // Constructor, initialize values
112     Jacobi() {
113       int i,j;
114       // allocate two dimensional arrays
115       temperature = new double*[blockDimX+2];
116       new_temperature = new double*[blockDimX+2];
117       for (i=0; i<blockDimX+2; i++) {
118         temperature[i] = new double[arrayDimY];
119         new_temperature[i] = new double[arrayDimY];
120       }
121       for(i=0;i<blockDimX+2; i++) {
122         for(j=0;j<arrayDimY; j++) {
123           temperature[i][j] = 0.5;
124           new_temperature[i][j] = 0.5;
125         }
126       }
127
128       arrived_top = 0;
129       arrived_bottom = 0;
130       iterations = 0;
131       constrainBC();
132     }
133
134     Jacobi(CkMigrateMessage* m) {}
135
136     ~Jacobi() { 
137       for (int i=0; i<blockDimX+2; i++) {
138         delete [] temperature[i];
139         delete [] new_temperature[i];
140       }
141       delete [] temperature; 
142       delete [] new_temperature; 
143     }
144
145     // Perform one iteration of work
146     void begin_iteration(void) {
147       iterations++;
148         
149       // Send my top edge
150       thisProxy(wrap_y(thisIndex)).receiveGhosts(BOTTOM, arrayDimY, &temperature[1][1]);
151       // Send my bottom edge
152       thisProxy(wrap_y(thisIndex)).receiveGhosts(TOP, arrayDimY, &temperature[blockDimX][1]);
153     }
154
155     void receiveGhosts(int dir, int size, double gh[]) {
156       int i, j;
157
158       switch(dir) {
159         case TOP:
160           arrived_top++;
161           for(j=0; j<size; j++)
162             temperature[0][j+1] = gh[j];
163           break;
164         case BOTTOM:
165           arrived_bottom++;
166           for(j=0; j<size; j++)
167             temperature[blockDimX+1][j+1] = gh[j];
168           break;
169         default:
170           CkAbort("ERROR\n");
171       }
172       check_and_compute();
173     }
174
175     void check_and_compute() {
176       double error = 0.0, max_error = 0.0;
177
178       if (arrived_top >=1 && arrived_bottom >= 1) {
179         arrived_top--;
180         arrived_bottom--;
181
182         compute_kernel();       
183
184         for(int i=1; i<blockDimX+1; i++) {
185           for(int j=0; j<arrayDimY; j++) {
186             error = fabs(new_temperature[i][j] - temperature[i][j]);
187             if(error > max_error) {
188               max_error = error;
189             }
190           }
191         }
192
193         double **tmp;
194         tmp = temperature;
195         temperature = new_temperature;
196         new_temperature = tmp;
197
198         constrainBC();
199
200         contribute(sizeof(double), &max_error, CkReduction::max_double,
201               CkCallback(CkIndex_Main::report(NULL), mainProxy));
202       }
203     }
204
205     // Check to see if we have received all neighbor values yet
206     // If all neighbor values have been received, we update our values and proceed
207     void compute_kernel()
208     {
209       for(int i=1; i<blockDimX+1; i++) {
210         for(int j=0; j<arrayDimY; j++) {
211           // update my value based on the surrounding values
212           new_temperature[i][j] = (temperature[i-1][j]+temperature[i+1][j]+temperature[i][j-1]+temperature[i][j+1]+temperature[i][j]) * 0.2;
213         }
214       }
215     }
216
217     // Enforce some boundary conditions
218     void constrainBC()
219     {
220      if(thisIndex <= num_chares/2) {
221         for(int i=1; i<=blockDimX; i++)
222           temperature[i][1] = 1.0;
223       }
224
225       if(thisIndex == num_chares-1) {
226         for(int j=arrayDimY/2; j<arrayDimY; j++)
227           temperature[blockDimX][j] = 0.0;
228       }
229     }
230
231 };
232
233 #include "jacobi1d.def.h"