commit a lot of unwanted changes.
[charm.git] / examples / bigsim / sdag / jacobi-no-redn / parallelJacobi.C
1 #include "parallelJacobi.h"
2
3 #if CMK_BLUEGENE_CHARM
4 #define BgElapse BgElapse
5 #define BgGetTime   BgGetTime
6 #else
7 #define BgElapse(x)
8 #define BgGetTime   CmiWallTimer
9 #endif
10
11 #define ITER    3
12
13 #define DEBUG   0
14
15 #define indexof(i,j,ydim) ( ((i)*(ydim)) + (j))
16
17 CProxy_Main globalMainProxy;
18 CProxy_Chunk chunk_arr;
19
20 double startTime;
21 int numFinished;
22
23 Chunk::Chunk(int t, int x, int y){
24
25   int xdim,ydim;
26   __sdag_init();
27   xdim = x;
28   ydim = y;
29   total = t; 
30   iterations =0;
31   myMax = 99999.999;
32
33
34   // CkPrintf("[%d] x is %d, y is %d, t is %d %f\n",CkMyPe(),x,y,t,BgGetTime());  
35   myxdim = int(xdim/total);
36   counter=0;
37
38   if(thisIndex == total-1) 
39     myxdim = xdim - myxdim*(total-1);    
40
41   myydim = ydim;
42
43   if((thisIndex != 0)&&(thisIndex != total-1)){
44     A = new double[(myxdim+2)*myydim];
45     B = new double[(myxdim+2)*myydim];
46       //Initialize everything to zero
47     for (int i=0; i<myxdim+2; i++)
48       for (int j=0; j<myydim; j++) 
49         A[indexof(i,j,ydim)] = B[indexof(i,j,ydim)] = 0.0;    
50   }
51   else {
52     A = new double[(myxdim+1)*myydim];
53     B = new double[(myxdim+1)*myydim];
54     //Initialize everything to zero
55   for (int i=0; i<myxdim+1; i++)
56     for (int j=0; j<myydim; j++) 
57       A[indexof(i,j,ydim)] = B[indexof(i,j,ydim)] = 0.0;  
58   }
59
60   usesAtSync = false;
61   //LBDatabase *lbdb = getLBDB();
62   //lbdb->SetLBPeriod(50);
63
64   BgElapse(.5e-5);
65 }
66
67
68 Chunk::Chunk(CkMigrateMessage *m){
69 }
70
71
72 void Chunk::resetBoundary() {
73
74   if((thisIndex !=0))
75     if(thisIndex < (int)(total/2))
76       for(int i=1;i<myxdim+1;i++)
77         A[indexof(i,0,myydim)] = 1.0;
78
79   if(thisIndex ==0){
80     //if(thisIndex < (int)(total/2))
81       for(int i=0;i<myxdim;i++)
82         A[indexof(i,0,myydim)] = 1.0;
83     
84     for (int i = 0;2*i<myydim; i++) 
85         A[indexof(0,i,myydim)] = 1.0;
86   
87 }
88 }
89
90
91 void Chunk::print() {
92
93   if ((myxdim>100)||(myydim>100)) return;
94
95 #if 1
96   CkPrintf("thisIndex = %d,myxdim=%d,myydim=%d\n",thisIndex,myxdim,myydim);
97
98   if(thisIndex !=0)
99     for (int i=0; i<myydim; i++) {
100       for (int j=1; j<myxdim+1; j++) 
101         CkPrintf("%lf ", A[indexof(j,i,myydim)]) ;
102       CkPrintf("\n");
103     }
104   else
105     for (int i=0; i<myydim; i++) {
106       for (int j=0; j<myxdim; j++) 
107         CkPrintf("%lf ", A[indexof(j,i,myydim)]) ;
108       CkPrintf("\n");
109     }
110 #endif
111 }
112
113
114
115 void Chunk::testEnd(){
116
117   if(iterations == ITER)  {
118                          
119     if(CkMyPe() != 0)
120       return;
121
122     CkPrintf("Numfin=%d, total=%d, Pes = %d\n",numFinished,total,CkNumPes());
123     if((CkMyPe()==0)&&(numFinished != total/CkNumPes()))
124       return;
125
126     double elapt = BgGetTime()-startTime;
127     CkPrintf("Finished in %fs %fs/step and iters is %d\n", elapt, elapt/iterations,iterations);
128     // BgPrint("ENDED at: %f\n");
129     CkExit();
130     return;
131   }
132
133 }
134
135
136 void Chunk::startWork(){
137
138   double* temp = (double*) malloc(sizeof(double)*myydim);
139
140   if(iterations == ITER)  {
141     //       BgPrint("FINISHED iterations at :%f\n");
142     if(CkMyPe() != 0)
143      return;
144
145     CkPrintf("Numfin=%d, total=%d, Pes = %d\n",numFinished,total,CkNumPes());
146     if((CkMyPe()==0)&&(numFinished != total/CkNumPes()))
147       return;
148
149     double elapt = BgGetTime()-startTime;
150     CkPrintf("Finished in %fs %fs/step and iters is %d\n", elapt, elapt/iterations,iterations);
151     //BgPrint0("ENDED at: %f\n");
152     CkExit();
153     return;
154   }
155
156      //CkPrintf("[%d]print in startWork %f\n", CkMyPe(), BgGetTime());
157 #if DEBUG
158   print();
159   CkPrintf("\n\n\n");
160 #endif
161
162   if(thisIndex >0){
163     for(int i=0;i<myydim;i++)
164       temp[i] = A[indexof(1,i,myydim)];
165     chunk_arr[thisIndex-1].getStripfromright(new (myydim,0) Msg(myydim,temp));
166   } 
167   else{
168     //Send dummy if thisIndex==0 for dagger to work
169     chunk_arr[total-1].getStripfromright(new (myydim,0) Msg(myydim,temp));
170   }
171
172   BgElapse(0.25e-5); //Time for sending to getStripfromright
173
174   if(thisIndex < total-1){
175   
176       for(int i=0;i<myydim;i++)
177         temp[i] = A[indexof(myxdim,i,myydim)];
178       chunk_arr[thisIndex+1].getStripfromleft(new (myydim,0) Msg(myydim,temp));
179   }
180   else{
181     //Send dummy if thisIndex==total-1:For dagger to work
182     chunk_arr[0].getStripfromleft(new (myydim,0) Msg(myydim,temp));
183   }
184
185   BgElapse(0.25e-5); //Time for sending to getStripfromleft
186   
187 #if DEBUG
188     CkPrintf("\n\nA in end of startWork is \n");
189   print();
190 #endif
191 }
192
193
194 void Chunk::doWork(){
195
196   double maxChange = 0.0;
197   double * temp;
198
199   iterations++;
200   if((iterations == ITER)&&(CkMyPe()==0))
201     numFinished++;
202   if (thisIndex == 0)  
203     CkPrintf("Iteration: %d\n",iterations);
204
205   resetBoundary();
206
207     if((thisIndex !=0)&&(thisIndex != total-1))
208       for (int i=1; i<myxdim+1; i++)
209         for (int j=1; j<myydim-1; j++) {
210           B[indexof(i,j,myydim)] = 
211             (0.2)*(A[indexof(i,  j,  myydim)] +
212                    A[indexof(i,  j+1,myydim)] +
213                    A[indexof(i,  j-1,myydim)] +
214                    A[indexof(i+1,j,  myydim)] +
215                    A[indexof(i-1,j,  myydim)]);
216
217           double change =  B[indexof(i,j,myydim)] - A[indexof(i,j,myydim)];
218           if (change < 0) change = - change;
219           if (change > maxChange) maxChange = change;
220         }
221       
222     if(thisIndex == 0)
223       for (int i=1; i<myxdim; i++)
224         for (int j=1; j<myydim-1; j++) {
225           B[indexof(i,j,myydim)] = 
226             (0.2)*(A[indexof(i,  j,  myydim)] +
227                    A[indexof(i,  j+1,myydim)] +
228                    A[indexof(i,  j-1,myydim)] +
229                    A[indexof(i+1,j,  myydim)] +
230                    A[indexof(i-1,j,  myydim)]);
231
232           double change =  B[indexof(i,j,myydim)] - A[indexof(i,j,myydim)];
233           if (change < 0) change = - change;
234           if (change > maxChange) maxChange = change;
235         }
236       
237     if(thisIndex == total-1) {
238       for (int i=1; i<myxdim; i++)
239         for (int j=1; j<myydim-1; j++) {
240           B[indexof(i,j,myydim)] = 
241             (0.2)*(A[indexof(i,  j,  myydim)] +
242                    A[indexof(i,  j+1,myydim)] +
243                    A[indexof(i,  j-1,myydim)] +
244                    A[indexof(i+1,j,  myydim)] +
245                    A[indexof(i-1,j,  myydim)]);
246
247           double change =  B[indexof(i,j,myydim)] - A[indexof(i,j,myydim)];
248           if (change < 0) change = - change;
249           if (change > maxChange) maxChange = change;
250         }
251     }
252   
253   temp = A;
254   A =B; 
255   B=temp;  
256
257   BgElapse(20e-5);
258
259 }
260
261
262 void Chunk::processStripfromleft(Msg* aMessage){
263
264 //Do nothing if this is 0 Pe because the message will be a dummy for the 0 Pe.
265
266  BgElapse(2.5e-5);
267  //BgPrint("%f:TEST STRING LEFT\n");
268   if(thisIndex !=0) {
269                                 
270     for(int i=0;i<myydim;i++)
271       A[indexof(0,i,myydim)] = aMessage->strip[i];
272   }
273   
274 }
275
276 void Chunk::processStripfromright(Msg* aMessage){
277
278   //Do nothing if this is Pe number:(total-1) because this will be a dummy message for that Pe.
279
280   BgElapse(2.5e-5);
281  //BgPrint("%f:TEST STRING RIGHT\n");
282   if(thisIndex != total -1){
283     if(thisIndex != 0)
284       for(int i=0;i<myydim;i++)
285         A[indexof(myxdim+1,i,myydim)] = aMessage->strip[i];
286     else
287       for(int i=0;i<myydim;i++)
288         A[indexof(myxdim,i,myydim)] = aMessage->strip[i];
289   }
290 }
291
292
293 Main::Main(CkArgMsg *m)
294 {
295   int x,y,k;
296
297   if(m->argc != 4) CkAbort("Wrong parameters\n");
298         
299   x = atoi(m->argv[1]);
300   y = atoi(m->argv[2]);
301   k = atoi(m->argv[3]);
302
303   if(x < k) CkAbort("Xdim must be greater than k");
304   if (k < CkNumPes() || k%CkNumPes()) CkAbort("k must be multiple of numPes.");
305
306   chunk_arr = CProxy_Chunk::ckNew(k,x,y,k);
307   //chunk_arr.setReductionClient(workover, (void*)NULL);
308   //   CkCallback *cb = new CkCallback(CkIndex_Chunk::workover(NULL), CkArrayIndex1D(0), chunk_arr);
309   //chunk_arr.ckSetReductionClient(cb);
310
311
312   BgElapse(.25e-5);
313   //BgPrint0("STARTING FIRST STEP AT:%f\n");
314
315
316   //  chunk_arr.stepOver(new VoidMsg);
317   chunk_arr.singleStep(new VoidMsg());
318   startTime = BgGetTime();
319   numFinished = 0;
320 }
321
322
323
324
325 #include "parallelJacobi.def.h"
326
327
328