Roll back to 12/04 to correct my accidental checkin of test code.
[charm.git] / examples / ampi / Cjacobi3D / jacobi.C
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include "mpi.h"
4
5 #if CMK_BLUEGENE_CHARM
6 extern void BgPrintf(char *);
7 #define BGPRINTF(x)  if (thisIndex == 0) BgPrintf(x);
8 #else
9 #define BGPRINTF(x)
10 #endif
11
12 #define DIMX 100
13 #define DIMY 100
14 #define DIMZ 100
15
16 int NX, NY, NZ;
17
18 class chunk {
19   public:
20     double t[DIMX+2][DIMY+2][DIMZ+2];
21     int xidx, yidx, zidx;
22     int xm, xp, ym, yp, zm, zp;
23     double sbxm[DIMY*DIMZ];
24     double sbxp[DIMY*DIMZ];
25     double sbym[DIMX*DIMZ];
26     double sbyp[DIMX*DIMZ];
27     double sbzm[DIMX*DIMY];
28     double sbzp[DIMX*DIMY];
29     double rbxm[DIMY*DIMZ];
30     double rbxp[DIMY*DIMZ];
31     double rbym[DIMX*DIMZ];
32     double rbyp[DIMX*DIMZ];
33     double rbzm[DIMX*DIMY];
34     double rbzp[DIMX*DIMY];
35 };
36
37 #ifdef AMPI
38 void chunk_pup(pup_er p, void *d)
39 {
40   chunk **cpp = (chunk **) d;
41   if(pup_isUnpacking(p))
42     *cpp = new chunk;
43   chunk *cp = *cpp;
44   pup_doubles(p, &cp->t[0][0][0], (DIMX+2)*(DIMY+2)*(DIMZ+2));
45   pup_int(p, &cp->xidx);
46   pup_int(p, &cp->yidx);
47   pup_int(p, &cp->zidx);
48   pup_int(p, &cp->xp);
49   pup_int(p, &cp->xm);
50   pup_int(p, &cp->yp);
51   pup_int(p, &cp->ym);
52   pup_int(p, &cp->zp);
53   pup_int(p, &cp->zm);
54   pup_doubles(p, cp->sbxm, (DIMY*DIMZ));
55   pup_doubles(p, cp->sbxp, (DIMY*DIMZ));
56   pup_doubles(p, cp->rbxm, (DIMY*DIMZ));
57   pup_doubles(p, cp->rbxp, (DIMY*DIMZ));
58   pup_doubles(p, cp->sbym, (DIMX*DIMZ));
59   pup_doubles(p, cp->sbyp, (DIMX*DIMZ));
60   pup_doubles(p, cp->rbym, (DIMX*DIMZ));
61   pup_doubles(p, cp->rbyp, (DIMX*DIMZ));
62   pup_doubles(p, cp->sbzm, (DIMX*DIMY));
63   pup_doubles(p, cp->sbzp, (DIMX*DIMY));
64   pup_doubles(p, cp->rbzm, (DIMX*DIMY));
65   pup_doubles(p, cp->rbzp, (DIMX*DIMY));
66   if(pup_isDeleting(p))
67     delete cp;
68 }
69 #endif
70
71 #define abs(x) ((x)<0.0 ? -(x) : (x))
72
73 int index1d(int ix, int iy, int iz)
74 {
75   return NY*NZ*ix + NZ*iy + iz;
76 }
77
78 void index3d(int index, int& ix, int& iy, int& iz)
79 {
80   ix = index/(NY*NZ);
81   iy = (index%(NY*NZ))/NZ;
82   iz = index%NZ;
83 }
84
85 static void copyout(double *d, double t[DIMX+2][DIMY+2][DIMZ+2],
86                     int sx, int ex, int sy, int ey, int sz, int ez)
87 {
88   int i, j, k;
89   int l = 0;
90   for(i=sx; i<=ex; i++)
91     for(j=sy; j<=ey; j++)
92       for(k=sz; k<=ez; k++, l++)
93         d[l] = t[i][j][k];
94 }
95
96 static void copyin(double *d, double t[DIMX+2][DIMY+2][DIMZ+2],
97                     int sx, int ex, int sy, int ey, int sz, int ez)
98 {
99   int i, j, k;
100   int l = 0;
101   for(i=sx; i<=ex; i++)
102     for(j=sy; j<=ey; j++)
103       for(k=sz; k<=ez; k++, l++)
104         t[i][j][k] = d[l];
105 }
106
107 int main(int ac, char** av)
108 {
109   int i,j,k,m,cidx;
110   int iter, niter;
111   MPI_Status status;
112   double error, tval, maxerr, tmpmaxerr, starttime, endtime, itertime;
113   chunk *cp;
114   int thisIndex, ierr, nblocks;
115
116   MPI_Init(&ac, &av);
117   MPI_Comm_rank(MPI_COMM_WORLD, &thisIndex);
118   MPI_Comm_size(MPI_COMM_WORLD, &nblocks);
119
120   if (ac < 4) {
121     if (thisIndex == 0)
122       printf("Usage: jacobi X Y Z [nIter].\n");
123     MPI_Finalize();
124   }
125   NX = atoi(av[1]);
126   NY = atoi(av[2]);
127   NZ = atoi(av[3]);
128   if (NX*NY*NZ != nblocks) {
129     if (thisIndex == 0) 
130       printf("%d x %d x %d != %d\n", NX,NY,NZ, nblocks);
131     MPI_Finalize();
132   }
133   if (ac == 5)
134     niter = atoi(av[4]);
135   else
136     niter = 20;
137
138 /*
139   if(thisIndex == 0)
140     niter = 20;
141 */
142
143   MPI_Bcast(&niter, 1, MPI_INT, 0, MPI_COMM_WORLD);
144
145   cp = new chunk;
146 #ifdef AMPI
147   MPI_Register((void*)&cp, (MPI_PupFn) chunk_pup);
148 #endif
149
150   index3d(thisIndex, cp->xidx, cp->yidx, cp->zidx);
151   cp->xp = index1d((cp->xidx+1)%NX,cp->yidx,cp->zidx);
152   cp->xm = index1d((cp->xidx+NX-1)%NX,cp->yidx,cp->zidx);
153   cp->yp = index1d(cp->xidx,(cp->yidx+1)%NY,cp->zidx);
154   cp->ym = index1d(cp->xidx,(cp->yidx+NY-1)%NY,cp->zidx);
155   cp->zp = index1d(cp->xidx,cp->yidx,(cp->zidx+1)%NZ);
156   cp->zm = index1d(cp->xidx,cp->yidx,(cp->zidx+NZ-1)%NZ);
157   for(i=1; i<=DIMZ; i++)
158     for(j=1; j<=DIMY; j++)
159       for(k=1; k<=DIMX; k++)
160         cp->t[k][j][i] = DIMY*DIMX*(i-1) + DIMX*(j-2) + (k-1);
161
162   MPI_Barrier(MPI_COMM_WORLD);
163   starttime = MPI_Wtime();
164
165   maxerr = 0.0;
166   for(iter=1; iter<=niter; iter++) {
167     BGPRINTF("interation starts at %f\n");
168     maxerr = 0.0;
169     copyout(cp->sbxm, cp->t, 1, 1, 1, DIMY, 1, DIMZ);
170     copyout(cp->sbxp, cp->t, DIMX, DIMX, 1, DIMY, 1, DIMZ);
171     copyout(cp->sbym, cp->t, 1, DIMX, 1, 1, 1, DIMZ);
172     copyout(cp->sbyp, cp->t, 1, DIMX, DIMY, DIMY, 1, DIMZ);
173     copyout(cp->sbzm, cp->t, 1, DIMX, 1, DIMY, 1, 1);
174     copyout(cp->sbzp, cp->t, 1, DIMX, 1, DIMY, DIMZ, DIMZ);
175
176     MPI_Request rreq[6];
177     MPI_Status rsts[6];
178
179     MPI_Irecv(cp->rbxp, DIMY*DIMZ, MPI_DOUBLE, cp->xp, 0, MPI_COMM_WORLD, &rreq[0]);
180     MPI_Irecv(cp->rbxm, DIMY*DIMZ, MPI_DOUBLE, cp->xm, 1, MPI_COMM_WORLD, &rreq[1]);
181     MPI_Irecv(cp->rbyp, DIMX*DIMZ, MPI_DOUBLE, cp->yp, 2, MPI_COMM_WORLD, &rreq[2]);
182     MPI_Irecv(cp->rbym, DIMX*DIMZ, MPI_DOUBLE, cp->ym, 3, MPI_COMM_WORLD, &rreq[3]);
183     MPI_Irecv(cp->rbzm, DIMX*DIMY, MPI_DOUBLE, cp->zm, 5, MPI_COMM_WORLD, &rreq[4]);
184     MPI_Irecv(cp->rbzp, DIMX*DIMY, MPI_DOUBLE, cp->zp, 4, MPI_COMM_WORLD, &rreq[5]);
185
186     MPI_Send(cp->sbxm, DIMY*DIMZ, MPI_DOUBLE, cp->xm, 0, MPI_COMM_WORLD);
187     MPI_Send(cp->sbxp, DIMY*DIMZ, MPI_DOUBLE, cp->xp, 1, MPI_COMM_WORLD);
188     MPI_Send(cp->sbym, DIMX*DIMZ, MPI_DOUBLE, cp->ym, 2, MPI_COMM_WORLD);
189     MPI_Send(cp->sbyp, DIMX*DIMZ, MPI_DOUBLE, cp->yp, 3, MPI_COMM_WORLD);
190     MPI_Send(cp->sbzm, DIMX*DIMY, MPI_DOUBLE, cp->zm, 4, MPI_COMM_WORLD);
191     MPI_Send(cp->sbzp, DIMX*DIMY, MPI_DOUBLE, cp->zp, 5, MPI_COMM_WORLD);
192
193     MPI_Waitall(6, rreq, rsts);
194
195     copyin(cp->sbxm, cp->t, 0, 0, 1, DIMY, 1, DIMZ);
196     copyin(cp->sbxp, cp->t, DIMX+1, DIMX+1, 1, DIMY, 1, DIMZ);
197     copyin(cp->sbym, cp->t, 1, DIMX, 0, 0, 1, DIMZ);
198     copyin(cp->sbyp, cp->t, 1, DIMX, DIMY+1, DIMY+1, 1, DIMZ);
199     copyin(cp->sbzm, cp->t, 1, DIMX, 1, DIMY, 0, 0);
200     copyin(cp->sbzp, cp->t, 1, DIMX, 1, DIMY, DIMZ+1, DIMZ+1);
201     if(iter > 25 &&  iter < 85 && thisIndex == 35)
202       m = 9;
203     else
204       m = 1;
205     for(; m>0; m--)
206       for(i=1; i<=DIMZ; i++)
207         for(j=1; j<=DIMY; j++)
208           for(k=1; k<=DIMX; k++) {
209             tval = (cp->t[k][j][i] + cp->t[k][j][i+1] +
210                  cp->t[k][j][i-1] + cp->t[k][j+1][i]+ 
211                  cp->t[k][j-1][i] + cp->t[k+1][j][i] + cp->t[k-1][j][i])/7.0;
212             error = abs(tval-cp->t[k][j][i]);
213             cp->t[k][j][i] = tval;
214             if (error > maxerr) maxerr = error;
215           }
216     MPI_Allreduce(&maxerr, &tmpmaxerr, 1, MPI_DOUBLE, MPI_MAX, 
217                    MPI_COMM_WORLD);
218     maxerr = tmpmaxerr;
219     endtime = MPI_Wtime();
220     itertime = endtime - starttime;
221     double  it;
222     MPI_Allreduce(&itertime, &it, 1, MPI_DOUBLE, MPI_SUM,
223                    MPI_COMM_WORLD);
224     itertime = it/nblocks;
225     if (thisIndex == 0)
226       printf("iter %d time: %lf maxerr: %lf\n", iter, itertime, maxerr);
227     starttime = MPI_Wtime();
228 #ifdef AMPI
229     if(iter%20 == 10) {
230       MPI_Migrate();
231     }
232 #endif
233   }
234   MPI_Finalize();
235   return 0;
236 }