fix bug and clean up jacobi examples
[charm.git] / examples / charm++ / jacobi2d-1d-decomposition / jacobi1d.C
1 /** \file jacobi1d.C
2  *  Author: Harshitha Menon, Yanhua Sun
3  *
4  *    ***********  ^
5  *    *         *  |
6  *    *         *  |
7  *    *         *  X
8  *    *         *  |
9  *    *         *  |
10  *    ***********  ~
11  *    <--- Y --->
12  *
13  *  This example does a 1D decomposition of a 2D data array
14  */
15
16 #include "jacobi1d.decl.h"
17 #include "TopoManager.h"
18
19 /*readonly*/ CProxy_Main mainProxy;
20 /*readonly*/ int blockDimX;
21 /*readonly*/ int arrayDimX;
22 /*readonly*/ int arrayDimY;
23
24 // specify the number of worker chares 
25 /*readonly*/ int numChares;
26 /*readonly*/ int maxiterations;
27
28 #define MAX_ITER                100
29 #define TOP             1
30 #define BOTTOM          2
31 const double THRESHOLD =  0.004;
32
33 class Main : public CBase_Main
34 {
35   double startTime;
36   double endTime;
37
38 public:
39   CProxy_Jacobi array;
40
41   Main(CkArgMsg* m) {
42     if (m->argc!=4 && m->argc!=5 ) {
43       CkPrintf("%s [array_size_X] [array_size_Y] [numChares]\n", m->argv[0]);
44       CkPrintf("OR %s [array_size_X] [array_size_Y] [numChares] maxiterations\n", m->argv[0]);
45       CkAbort("Abort");
46     }
47
48     CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
49     arrayDimX = atoi(m->argv[1]);
50     arrayDimY = atoi(m->argv[2]);
51     numChares = atoi(m->argv[3]);
52     if (arrayDimX < numChares || arrayDimX % numChares != 0)
53       CkAbort("array_size_X % numChares != 0!");
54     blockDimX = arrayDimX / numChares;
55
56     maxiterations = MAX_ITER;
57     if(m->argc==5)
58       maxiterations = atoi(m->argv[4]); 
59     // store the main proxy
60     mainProxy = thisProxy;
61
62     // print info
63     CkPrintf("Running Jacobi on %d processors with (%d) elements\n", CkNumPes(), numChares);
64     CkPrintf("Array Dimensions: %d %d\n", arrayDimX, arrayDimY);
65     CkPrintf("Block Dimensions: %d %d\n", arrayDimX, blockDimX);
66
67     // Create new array of worker chares
68     array = CProxy_Jacobi::ckNew(numChares);
69
70     //Start the computation
71     startTime = CkWallTimer();
72     array.run();
73   }
74
75   void done(int totalIter) {
76     if(totalIter >= maxiterations)
77       CkPrintf("Finish due to max iterations %d, total time %.3f seconds. \n", totalIter, CkWallTimer()-startTime); 
78     else
79       CkPrintf("Finish due to convergence, iterations %d, total time %.3f seconds. \n", totalIter, CkWallTimer()-startTime); 
80     CkExit();
81   }
82 };
83
84 class Jacobi: public CBase_Jacobi {
85   Jacobi_SDAG_CODE
86
87 public:
88
89     double **temperature;
90     double **new_temperature;
91     int imsg;
92     int iterations;
93     int neighbors;
94     double max_error;
95     int converged;
96     bool topBound, bottomBound;
97     int istart, ifinish;
98
99     // Constructor, initialize values
100     Jacobi() {
101       int i,j;
102       temperature = new double*[blockDimX+2];
103       new_temperature = new double*[blockDimX+2];
104       for (i=0; i<blockDimX+2; i++) {
105         temperature[i] = new double[arrayDimY];
106         new_temperature[i] = new double[arrayDimY];
107       }
108       for (i=0; i<blockDimX+2; i++) {
109         for(j=0;j<arrayDimY; j++) {
110           temperature[i][j] = 0;
111           new_temperature[i][j] = 0;
112         }
113       }
114       converged = 0;
115       topBound = bottomBound = false;
116       neighbors = 0;
117       if(thisIndex == 0) {
118         topBound = true;
119         istart = 2;
120       }
121       else {
122         neighbors++;
123         istart = 1;
124       }
125
126       if(thisIndex == numChares -1){
127         bottomBound = true;
128         ifinish = blockDimX;
129       }
130       else {
131         neighbors++;
132         ifinish = blockDimX+1;
133       }
134
135       iterations = 0;
136       constrainBC();
137     }
138
139     Jacobi(CkMigrateMessage* m) {}
140
141     ~Jacobi() { 
142       for (int i=0; i<blockDimX+2; i++) {
143         delete [] temperature[i];
144         delete [] new_temperature[i];
145       }
146       delete [] temperature; 
147       delete [] new_temperature; 
148     }
149
150
151     void pup(PUP::er &p)
152     {
153       CBase_Jacobi::pup(p);
154       __sdag_pup(p);
155       p|imsg;
156       p|iterations;
157       p|neighbors;
158       p|max_error;
159       p|converged;
160       p|topBound;
161       p|bottomBound;
162       p|istart;
163       p|ifinish;
164
165       size_t size = (blockDimX+2) * (arrayDimY);
166       if (p.isUnpacking()) {
167         temperature = new double*[blockDimX+2];
168         new_temperature = new double*[blockDimX+2];
169
170         for (int i=0; i<blockDimX+2; i++) {
171           temperature[i] = new double[arrayDimY];
172           new_temperature[i] = new double[arrayDimY];
173         }
174       }
175
176       for(int i=0;i<blockDimX+2; i++) {
177         for(int j=0;j<arrayDimY; j++) {
178           p|temperature[i][j];
179           p|new_temperature[i][j];
180         }
181       }
182     }
183
184     void check_and_compute() {
185       double error = 0.0;
186       double **tmp;
187
188       max_error = 0.0;
189       for(int i=istart; i<ifinish; i++) {
190         for(int j=1; j<arrayDimY-1; j++) {
191           // update my value based on the surrounding values
192           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;
193           error = fabs(new_temperature[i][j] - temperature[i][j]);
194           if(error > max_error) {
195             max_error = error;
196           }
197         }
198       }
199
200       tmp = temperature;
201       temperature = new_temperature;
202       new_temperature = tmp;
203     }
204
205     // Enforce some boundary conditions
206     void constrainBC() {
207       int i;
208       //top boundary 
209       if(topBound) {
210         for(i=0;i<arrayDimY; i++) {
211           temperature[1][i] = 1.0;
212           new_temperature[1][i] = 1.0;
213         }
214       }
215       //left, right boundary
216       for(i=0;i<blockDimX+2; i++) 
217       {
218         temperature[i][0] = 1.0;
219         new_temperature[i][0] = 1.0;
220         temperature[i][arrayDimY-1] = 1.0;
221         new_temperature[i][arrayDimY-1] = 1.0;
222       }
223       //bottom boundary
224       if(bottomBound) {
225         for(i=0;i<arrayDimY; i++) {
226           temperature[blockDimX][i] = 1.0;
227           new_temperature[blockDimX][i] = 1.0;
228         }
229       }
230     } 
231     // for debugging
232     void dumpMatrix(double **matrix)
233     {
234       CkPrintf("\n\n[%d]\n",thisIndex);
235       for(int i=0; i<blockDimX+2;++i)
236       {
237         for(int j=0; j<arrayDimY;++j)
238         {
239           CkPrintf("%0.6lf ",matrix[i][j]);
240         }
241         CkPrintf("\n");
242       }
243     }
244
245 };
246
247 #include "jacobi1d.def.h"