bug fix
[charm.git] / examples / ampi / Cjacobi3D / jacobi-cpp.C
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include "mpi.h"
4 #ifdef AMPI
5 #include "charm++.h"
6 #endif
7
8 #if CMK_BLUEGENE_CHARM
9 extern void BgPrintf(char *);
10 #define BGPRINTF(x)  if (thisIndex == 0) BgPrintf(x);
11 #else
12 #define BGPRINTF(x)
13 #endif
14
15 int DIM, DIMX, DIMY, DIMZ, NX, NY, NZ;
16
17 class chunk {
18  public:
19   int dimx, dimy, dimz;
20   int xidx, yidx, zidx;
21   int xm, xp, ym, yp, zm, zp;
22   double ***t; 
23   double *sbxm;
24   double *sbxp;
25   double *sbym;
26   double *sbyp;
27   double *sbzm;
28   double *sbzp;
29   double *rbxm;
30   double *rbxp;
31   double *rbym;
32   double *rbyp;
33   double *rbzm;
34   double *rbzp;
35
36  private:
37   void create(){
38     t = new double** [dimx+2];
39     for(int i=0;i<dimx+2;i++){
40       t[i] = new double* [dimy+2];
41       for(int j=0;j<dimy+2;j++){
42         t[i][j] = new double [dimz+2];
43       }
44     }
45     sbxm = new double [dimy*dimz];
46     sbxp = new double [dimy*dimz];
47     sbym = new double [dimx*dimz];
48     sbyp = new double [dimx*dimz];
49     sbzm = new double [dimx*dimy];
50     sbzp = new double [dimx*dimy];
51     rbxm = new double [dimy*dimz];
52     rbxp = new double [dimy*dimz];
53     rbym = new double [dimx*dimz];
54     rbyp = new double [dimx*dimz];
55     rbzm = new double [dimx*dimy];
56     rbzp = new double [dimx*dimy];
57   }
58   void destroy(){
59     for(int i=0;i<dimx+2;i++){
60       for(int j=0;j<dimy+2;j++){
61         delete [] t[i][j];
62       }
63       delete [] t[i];
64     }
65     delete [] t;
66     delete [] sbxm;
67     delete [] sbxp;
68     delete [] sbym;
69     delete [] sbyp;
70     delete [] sbzm;
71     delete [] sbzp;
72     delete [] rbxm;
73     delete [] rbxp;
74     delete [] rbym;
75     delete [] rbyp;
76     delete [] rbzm;
77     delete [] rbzp;
78   }
79   inline int index1d(int ix, int iy, int iz){
80     return NY*NZ*ix + NZ*iy + iz;
81   }
82   
83   inline void index3d(int index, int& ix, int& iy, int& iz){
84     ix = index/(NY*NZ);
85     iy = (index%(NY*NZ))/NZ;
86     iz = index%NZ;
87   }
88   
89  public:
90   chunk(int x, int y, int z, int rank){
91     int i,j,k;
92     dimx = x; dimy = y; dimz = z;
93     create();
94     indexing(rank);
95     for(i=1; i<=DIMX; i++)
96       for(j=1; j<=DIMY; j++)
97         for(k=1; k<=DIMZ; k++)
98           t[i][j][k] = DIMY*DIMZ*(i-1) + DIMZ*(j-1) + (k-1);
99   }
100   ~chunk(){
101     destroy();
102   }
103 #ifdef AMPI
104   void pup(PUP::er& p){
105     p|dimx;p|dimy;p|dimz;
106     p|xidx;p|yidx;p|zidx;
107     p|xp;p|xm;p|yp;p|ym;p|zp;p|zm;
108     if(p.isUnpacking())
109       create();
110     for(int i=0;i<dimx+2;i++){
111       for(int j=0;j<dimy+2;j++){
112         p(t[i][j],dimz+2);
113       }
114     }
115     p(sbxm,dimy*dimz);
116     p(sbxp,dimy*dimz);
117     p(sbym,dimx*dimz);
118     p(sbyp,dimx*dimz);
119     p(sbzm,dimx*dimy);
120     p(sbzp,dimx*dimy);
121     p(rbxm,dimy*dimz);
122     p(rbxp,dimy*dimz);
123     p(rbym,dimx*dimz);
124     p(rbyp,dimx*dimz);
125     p(rbzm,dimx*dimy);
126     p(rbzp,dimx*dimy);
127     if(p.isDeleting())
128       destroy();
129   }
130 #endif
131   void copyin(){
132     int i, j;
133     int l = 0;
134     for(i=1;i<=dimy;i++)
135       for(j=1;j<=dimz;j++,l++){
136         t[0][i][j] = sbxm[l];
137         t[dimx+1][i][j] = sbxp[l];
138       }
139     l = 0;
140     for(i=1;i<=dimx;i++)
141       for(j=1;j<=dimz;j++,l++){
142         t[i][0][j] = sbym[l];
143         t[i][dimy+1][j] = sbyp[l];
144       }
145     l = 0;
146     for(i=1;i<=dimx;i++)
147       for(j=1;j<=dimy;j++,l++){
148         t[i][j][0] = sbzm[l];
149         t[i][j][dimz+1] = sbzp[l];
150       }
151   }
152
153   void copyout(){
154     int i, j;
155     int l = 0;
156     for(i=1;i<=dimy;i++)
157       for(j=1;j<=dimz;j++,l++){
158         sbxm[l] = t[1][i][j];
159         sbxp[l] = t[dimx][i][j];
160       }
161     l = 0;
162     for(i=1;i<=dimx;i++)
163       for(j=1;j<=dimz;j++,l++){
164         sbym[l] = t[i][1][j];
165         sbyp[l] = t[i][dimy][j];
166       }
167     l = 0;
168     for(i=1;i<=dimx;i++)
169       for(j=1;j<=dimy;j++,l++){
170         sbzm[l] = t[i][j][1];
171         sbzp[l] = t[i][j][dimz];
172       }
173   }
174
175   void indexing(int rank){
176     index3d(rank,xidx,yidx,zidx);
177     xp = index1d((xidx+1)%NX,yidx,zidx);
178     xm = index1d((xidx+NX-1)%NX,yidx,zidx);
179     yp = index1d(xidx,(yidx+1)%NY,zidx);
180     ym = index1d(xidx,(yidx+NY-1)%NY,zidx);
181     zp = index1d(xidx,yidx,(zidx+1)%NZ);
182     zm = index1d(xidx,yidx,(zidx+NZ-1)%NZ);
183   }
184
185   double calc(){
186     double error = 0.0, maxerr = 0.0, tval;
187     int i,j,k;
188     for(k=1; k<=DIMX; k++) 
189       for(j=1; j<=DIMY; j++)
190         for(i=1; i<=DIMZ; i++){
191           tval = (t[k][j][i] + t[k][j][i+1] + t[k][j][i-1] + t[k][j+1][i]+ 
192                   t[k][j-1][i] + t[k+1][j][i] + t[k-1][j][i])/7.0;
193           error = abs(tval-t[k][j][i]);
194           t[k][j][i] = tval;
195           if (error > maxerr) maxerr = error;
196         }
197     return maxerr;
198   }
199 };
200
201 #define abs(x) ((x)<0.0 ? -(x) : (x))
202
203 int main(int ac, char** av)
204 {
205   int i,j,k,m,cidx;
206   int iter, niter;
207   MPI_Status status;
208   double error, tval, maxerr, tmpmaxerr, starttime, endtime, itertime;
209   chunk *cp;
210   int thisIndex, ierr, nblocks;
211
212   MPI_Init(&ac, &av);
213   MPI_Comm_rank(MPI_COMM_WORLD, &thisIndex);
214   MPI_Comm_size(MPI_COMM_WORLD, &nblocks);
215
216   if (ac < 5) {
217     if (thisIndex == 0)
218       printf("Usage: jacobi DIM X Y Z [nIter].\n");
219     MPI_Finalize();
220   }
221   DIM = atoi(av[1]);
222   NX = atoi(av[2]);
223   NY = atoi(av[3]);
224   NZ = atoi(av[4]);
225   if (NX*NY*NZ != nblocks) {
226     if (thisIndex == 0) 
227       printf("%d x %d x %d != %d\n", NX,NY,NZ, nblocks);
228     MPI_Finalize();
229   }
230   if (ac == 6)
231     niter = atoi(av[5]);
232   else
233     niter = 10;
234
235   DIMX = DIM/NX;
236   DIMY = DIM/NY;
237   DIMZ = DIM/NZ;
238
239   MPI_Bcast(&niter, 1, MPI_INT, 0, MPI_COMM_WORLD);
240
241   cp = new chunk(DIMX,DIMY,DIMZ,thisIndex);
242
243   MPI_Barrier(MPI_COMM_WORLD);
244   starttime = MPI_Wtime();
245
246   for(iter=1; iter<=niter; iter++) {
247     BGPRINTF("interation starts at %f\n");
248     maxerr = 0.0;
249
250     cp->copyout();
251
252     MPI_Request rreq[6];
253     MPI_Status rsts[6];
254
255     MPI_Irecv(cp->rbxp, DIMY*DIMZ, MPI_DOUBLE, cp->xp, 0, MPI_COMM_WORLD, &rreq[0]);
256     MPI_Irecv(cp->rbxm, DIMY*DIMZ, MPI_DOUBLE, cp->xm, 1, MPI_COMM_WORLD, &rreq[1]);
257     MPI_Irecv(cp->rbyp, DIMX*DIMZ, MPI_DOUBLE, cp->yp, 2, MPI_COMM_WORLD, &rreq[2]);
258     MPI_Irecv(cp->rbym, DIMX*DIMZ, MPI_DOUBLE, cp->ym, 3, MPI_COMM_WORLD, &rreq[3]);
259     MPI_Irecv(cp->rbzp, DIMX*DIMY, MPI_DOUBLE, cp->zp, 4, MPI_COMM_WORLD, &rreq[4]);
260     MPI_Irecv(cp->rbzm, DIMX*DIMY, MPI_DOUBLE, cp->zm, 5, MPI_COMM_WORLD, &rreq[5]);
261
262     MPI_Send(cp->sbxm, DIMY*DIMZ, MPI_DOUBLE, cp->xm, 0, MPI_COMM_WORLD);
263     MPI_Send(cp->sbxp, DIMY*DIMZ, MPI_DOUBLE, cp->xp, 1, MPI_COMM_WORLD);
264     MPI_Send(cp->sbym, DIMX*DIMZ, MPI_DOUBLE, cp->ym, 2, MPI_COMM_WORLD);
265     MPI_Send(cp->sbyp, DIMX*DIMZ, MPI_DOUBLE, cp->yp, 3, MPI_COMM_WORLD);
266     MPI_Send(cp->sbzm, DIMX*DIMY, MPI_DOUBLE, cp->zm, 4, MPI_COMM_WORLD);
267     MPI_Send(cp->sbzp, DIMX*DIMY, MPI_DOUBLE, cp->zp, 5, MPI_COMM_WORLD);
268
269     MPI_Waitall(6, rreq, rsts);
270
271     cp->copyin();
272
273     maxerr = cp->calc(); 
274
275     MPI_Allreduce(&maxerr, &tmpmaxerr, 1, MPI_DOUBLE, MPI_MAX, 
276                    MPI_COMM_WORLD);
277     maxerr = tmpmaxerr;
278     endtime = MPI_Wtime();
279     itertime = endtime - starttime;
280     if (thisIndex == 0)
281       printf("iter %d time: %lf maxerr: %lf\n", iter, itertime, maxerr);
282     starttime = MPI_Wtime();
283
284 #ifdef AMPI
285     if(iter%20 == 10) {
286       MPI_Migrate();
287     }
288 #endif
289   }
290   MPI_Finalize();
291   return 0;
292 }