7dbe3e9c660671f4872053156e83099a87397715
[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   // CkPrintf("[%d] x is %d, y is %d, t is %d %f\n",CkMyPe(),x,y,t,BgGetTime());  
34   myxdim = int(xdim/total);
35   counter=0;
36
37   if(thisIndex == total-1) 
38     myxdim = xdim - myxdim*(total-1);    
39
40   myydim = ydim;
41
42   if((thisIndex != 0)&&(thisIndex != total-1)){
43     A = new double[(myxdim+2)*myydim];
44     B = new double[(myxdim+2)*myydim];
45       //Initialize everything to zero
46     for (int i=0; i<myxdim+2; i++)
47       for (int j=0; j<myydim; j++) 
48         A[indexof(i,j,ydim)] = B[indexof(i,j,ydim)] = 0.0;    
49   }
50   else {
51     A = new double[(myxdim+1)*myydim];
52     B = new double[(myxdim+1)*myydim];
53     //Initialize everything to zero
54   for (int i=0; i<myxdim+1; i++)
55     for (int j=0; j<myydim; j++) 
56       A[indexof(i,j,ydim)] = B[indexof(i,j,ydim)] = 0.0;  
57   }
58
59   usesAtSync = false;
60   //LBDatabase *lbdb = getLBDB();
61   //lbdb->SetLBPeriod(50);
62
63   BgElapse(.5e-5);
64 }
65
66
67 Chunk::Chunk(CkMigrateMessage *m){
68 }
69
70
71 void Chunk::resetBoundary() {
72
73   if((thisIndex !=0))
74     if(thisIndex < (int)(total/2))
75       for(int i=1;i<myxdim+1;i++)
76         A[indexof(i,0,myydim)] = 1.0;
77
78   if(thisIndex ==0){
79     //if(thisIndex < (int)(total/2))
80       for(int i=0;i<myxdim;i++)
81         A[indexof(i,0,myydim)] = 1.0;
82     
83     for (int i = 0;2*i<myydim; i++) 
84         A[indexof(0,i,myydim)] = 1.0;
85   
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 int xx[1024*4+100];
192 for (int i=0;i<1024*4+100;i++) xx[i]=i;
193 }
194
195
196 void Chunk::doWork(){
197
198   double maxChange = 0.0;
199   double * temp;
200
201   iterations++;
202   if((iterations == ITER)&&(CkMyPe()==0))
203     numFinished++;
204   if (thisIndex == 0)  
205     CkPrintf("Iteration: %d\n",iterations);
206
207   resetBoundary();
208
209     if((thisIndex !=0)&&(thisIndex != total-1))
210       for (int i=1; i<myxdim+1; i++)
211         for (int j=1; j<myydim-1; j++) {
212           B[indexof(i,j,myydim)] = 
213             (0.2)*(A[indexof(i,  j,  myydim)] +
214                    A[indexof(i,  j+1,myydim)] +
215                    A[indexof(i,  j-1,myydim)] +
216                    A[indexof(i+1,j,  myydim)] +
217                    A[indexof(i-1,j,  myydim)]);
218
219           double change =  B[indexof(i,j,myydim)] - A[indexof(i,j,myydim)];
220           if (change < 0) change = - change;
221           if (change > maxChange) maxChange = change;
222         }
223       
224     if(thisIndex == 0)
225       for (int i=1; i<myxdim; i++)
226         for (int j=1; j<myydim-1; j++) {
227           B[indexof(i,j,myydim)] = 
228             (0.2)*(A[indexof(i,  j,  myydim)] +
229                    A[indexof(i,  j+1,myydim)] +
230                    A[indexof(i,  j-1,myydim)] +
231                    A[indexof(i+1,j,  myydim)] +
232                    A[indexof(i-1,j,  myydim)]);
233
234           double change =  B[indexof(i,j,myydim)] - A[indexof(i,j,myydim)];
235           if (change < 0) change = - change;
236           if (change > maxChange) maxChange = change;
237         }
238       
239     if(thisIndex == total-1) {
240       for (int i=1; i<myxdim; i++)
241         for (int j=1; j<myydim-1; j++) {
242           B[indexof(i,j,myydim)] = 
243             (0.2)*(A[indexof(i,  j,  myydim)] +
244                    A[indexof(i,  j+1,myydim)] +
245                    A[indexof(i,  j-1,myydim)] +
246                    A[indexof(i+1,j,  myydim)] +
247                    A[indexof(i-1,j,  myydim)]);
248
249           double change =  B[indexof(i,j,myydim)] - A[indexof(i,j,myydim)];
250           if (change < 0) change = - change;
251           if (change > maxChange) maxChange = change;
252         }
253     }
254   
255   temp = A;
256   A =B; 
257   B=temp;  
258
259   BgElapse(20e-5);
260
261 }
262
263
264 void Chunk::processStripfromleft(Msg* aMessage){
265
266 //Do nothing if this is 0 Pe because the message will be a dummy for the 0 Pe.
267
268  BgElapse(2.5e-5);
269  //BgPrint("%f:TEST STRING LEFT\n");
270   if(thisIndex !=0) {
271                                 
272     for(int i=0;i<myydim;i++)
273       A[indexof(0,i,myydim)] = aMessage->strip[i];
274   }
275   
276 }
277
278 void Chunk::processStripfromright(Msg* aMessage){
279
280   //Do nothing if this is Pe number:(total-1) because this will be a dummy message for that Pe.
281
282   BgElapse(2.5e-5);
283  //BgPrint("%f:TEST STRING RIGHT\n");
284   if(thisIndex != total -1){
285     if(thisIndex != 0)
286       for(int i=0;i<myydim;i++)
287         A[indexof(myxdim+1,i,myydim)] = aMessage->strip[i];
288     else
289       for(int i=0;i<myydim;i++)
290         A[indexof(myxdim,i,myydim)] = aMessage->strip[i];
291   }
292 }
293
294
295 Main::Main(CkArgMsg *m)
296 {
297   int x,y,k;
298
299   if(m->argc != 4) CkAbort("Wrong parameters\n");
300         
301   x = atoi(m->argv[1]);
302   y = atoi(m->argv[2]);
303   k = atoi(m->argv[3]);
304
305   if(x < k) CkAbort("Xdim must be greater than k");
306   if (k < CkNumPes() || k%CkNumPes()) CkAbort("k must be multiple of numPes.");
307
308   chunk_arr = CProxy_Chunk::ckNew(k,x,y,k);
309   //chunk_arr.setReductionClient(workover, (void*)NULL);
310   //   CkCallback *cb = new CkCallback(CkIndex_Chunk::workover(NULL), CkArrayIndex1D(0), chunk_arr);
311   //chunk_arr.ckSetReductionClient(cb);
312
313
314   BgElapse(.25e-5);
315   //BgPrint0("STARTING FIRST STEP AT:%f\n");
316
317
318   //  chunk_arr.stepOver(new VoidMsg);
319   chunk_arr.singleStep(new VoidMsg());
320   startTime = BgGetTime();
321   numFinished = 0;
322 }
323
324
325
326
327 #include "parallelJacobi.def.h"
328
329
330