faa8544408efbf8e9bdbb438c62a0af29a5b1f3a
[charm.git] / examples / charm++ / topology / jacobi3d / jacobi3d.C
1 /** \file jacobi3d.C
2  *  Author: Abhinav S Bhatele
3  *  Date Created: October 24th, 2007
4  *
5  *  This does a topological placement for a 3d jacobi.
6  *
7  *      
8  *            *****************
9  *         *               *  *
10  *   ^  *****************     *
11  *   |  *               *     *
12  *   |  *               *     *
13  *   |  *               *     *
14  *   Y  *               *     *
15  *   |  *               *     *
16  *   |  *               *     *
17  *   |  *               *  * 
18  *   ~  *****************    Z
19  *      <------ X ------> 
20  *
21  *   X: left, right --> wrap_x
22  *   Y: top, bottom --> wrap_y
23  *   Z: front, back --> wrap_z
24  */
25
26 #include "jacobi3d.decl.h"
27 #include "TopoManager.h"
28
29 // See README for documentation
30
31 /*readonly*/ CProxy_Main mainProxy;
32 /*readonly*/ int arrayDimX;
33 /*readonly*/ int arrayDimY;
34 /*readonly*/ int arrayDimZ;
35 /*readonly*/ int blockDimX;
36 /*readonly*/ int blockDimY;
37 /*readonly*/ int blockDimZ;
38 /*readonly*/ int torusDimX;
39 /*readonly*/ int torusDimY;
40 /*readonly*/ int torusDimZ;
41
42 // specify the number of worker chares in each dimension
43 /*readonly*/ int num_chare_x;
44 /*readonly*/ int num_chare_y;
45 /*readonly*/ int num_chare_z;
46
47 static unsigned long next = 1;
48
49 int myrand(int numpes) {
50   next = next * 1103515245 + 12345;
51   return((unsigned)(next/65536) % numpes);
52 }
53
54 #define SMPWAYX                 1
55 #define SMPWAYY                 2
56 #define SMPWAYZ                 2
57 #define USE_TOPOMAP             0
58 #define USE_RRMAP               0
59 #define USE_BLOCKMAP            1
60 #define USE_BLOCK_RNDMAP        0
61 #define USE_BLOCK_RRMAP         1
62 #define USE_SMPMAP              0
63 #define USE_3D_ARRAYS           0
64
65 // We want to wrap entries around, and because mod operator % 
66 // sometimes misbehaves on negative values. -1 maps to the highest value.
67 #define wrap_x(a)       (((a)+num_chare_x)%num_chare_x)
68 #define wrap_y(a)       (((a)+num_chare_y)%num_chare_y)
69 #define wrap_z(a)       (((a)+num_chare_z)%num_chare_z)
70
71 #if USE_3D_ARRAYS
72 #define index(a, b, c)  a][b][c 
73 #else
74 #define index(a, b, c)  ( (a)*(blockDimY+2)*(blockDimZ+2) + (b)*(blockDimZ+2) + (c) )
75 #endif
76
77 #define MAX_ITER                21
78 #define LEFT                    1
79 #define RIGHT                   2
80 #define TOP                     3
81 #define BOTTOM                  4
82 #define FRONT                   5
83 #define BACK                    6
84 #define DIVIDEBY7               0.14285714285714285714
85
86 double startTime;
87 double endTime;
88
89 /** \class Main
90  *
91  */
92
93 class Main : public CBase_Main {
94   public:
95     CProxy_Jacobi array;
96     int iterations;
97
98     Main(CkArgMsg* m) {
99       if ( (m->argc != 3) && (m->argc != 7) && (m->argc != 10) ) {
100         CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
101         CkPrintf("OR %s [array_size_X] [array_size_Y] [array_size_Z] [block_size_X] [block_size_Y] [block_size_Z]\n", m->argv[0]);
102         CkAbort("Abort");
103       }
104
105       // set iteration counter to zero
106       iterations = 0;
107
108       // store the main proxy
109       mainProxy = thisProxy;
110         
111       if(m->argc == 3) {
112         arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
113         blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]); 
114       }
115       else if (m->argc == 7) {
116         arrayDimX = atoi(m->argv[1]);
117         arrayDimY = atoi(m->argv[2]);
118         arrayDimZ = atoi(m->argv[3]);
119         blockDimX = atoi(m->argv[4]); 
120         blockDimY = atoi(m->argv[5]); 
121         blockDimZ = atoi(m->argv[6]);
122       } else {
123         arrayDimX = atoi(m->argv[1]);
124         arrayDimY = atoi(m->argv[2]);
125         arrayDimZ = atoi(m->argv[3]);
126         blockDimX = atoi(m->argv[4]); 
127         blockDimY = atoi(m->argv[5]); 
128         blockDimZ = atoi(m->argv[6]);
129         torusDimX = atoi(m->argv[7]);
130         torusDimY = atoi(m->argv[8]);
131         torusDimZ = atoi(m->argv[9]);
132       }
133
134       if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
135         CkAbort("array_size_X % block_size_X != 0!");
136       if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
137         CkAbort("array_size_Y % block_size_Y != 0!");
138       if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
139         CkAbort("array_size_Z % block_size_Z != 0!");
140
141       num_chare_x = arrayDimX / blockDimX;
142       num_chare_y = arrayDimY / blockDimY;
143       num_chare_z = arrayDimZ / blockDimZ;
144
145       // print info
146       CkPrintf("Running Jacobi on %d processors with (%d, %d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y, num_chare_z);
147       CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
148       CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
149
150       // Create new array of worker chares
151 #if USE_TOPOMAP || USE_RRMAP || USE_BLOCKMAP || USE_BLOCK_RRMAP || USE_BLOCK_RNDMAP || USE_SMPMAP
152        CkPrintf("Topology Mapping is being done ... %d %d %d %d\n", USE_TOPOMAP, USE_RRMAP, USE_BLOCKMAP, USE_SMPMAP);
153       CProxy_JacobiMap map = CProxy_JacobiMap::ckNew(num_chare_x, num_chare_y, num_chare_z, torusDimX, torusDimY, torusDimZ);
154       CkArrayOptions opts(num_chare_x, num_chare_y, num_chare_z);
155       opts.setMap(map);
156       array = CProxy_Jacobi::ckNew(opts);
157 #else
158       array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y, num_chare_z);
159 #endif
160
161       TopoManager tmgr;
162       CkArray *jarr = array.ckLocalBranch();
163       int jmap[num_chare_x][num_chare_y][num_chare_z];
164
165       int hops=0, p;
166       for(int i=0; i<num_chare_x; i++)
167         for(int j=0; j<num_chare_y; j++)
168           for(int k=0; k<num_chare_z; k++) {
169             jmap[i][j][k] = jarr->procNum(CkArrayIndex3D(i, j, k));
170           }
171
172       for(int i=0; i<num_chare_x; i++)
173         for(int j=0; j<num_chare_y; j++)
174           for(int k=0; k<num_chare_z; k++) {
175             p = jmap[i][j][k];
176             hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i+1)][j][k]);
177             hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i-1)][j][k]);
178             hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j+1)][k]);
179             hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j-1)][k]);
180             hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k+1)]);
181             hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k-1)]);
182           }
183       CkPrintf("Total Hops: %d\n", hops);
184
185       //Start the computation
186       array.begin_iteration();
187     }
188
189     // Each worker reports back to here when it completes an iteration
190     void report() {
191       if (iterations == 0) {
192         iterations++;
193         startTime = CkWallTimer();
194         array.begin_iteration();
195       }
196       else {
197         endTime = CkWallTimer();
198         CkPrintf("TIME : %f\n", (endTime - startTime)/(MAX_ITER-1));
199         CkExit();
200       }
201     }
202 };
203
204 /** \class Jacobi
205  *
206  */
207
208 class Jacobi: public CBase_Jacobi {
209   public:
210     int arrived_left;
211     int arrived_right;
212     int arrived_top;
213     int arrived_bottom;
214     int arrived_front;
215     int arrived_back;
216     int iterations;
217     int msgs;
218
219 #if USE_3D_ARRAYS
220     double ***temperature;
221     double ***new_temperature;
222 #else
223     double *temperature;
224     double *new_temperature;
225 #endif
226
227     // Constructor, initialize values
228     Jacobi() {
229         int i, j, k;
230         // allocate a three dimensional array
231 #if USE_3D_ARRAYS
232         temperature = new double**[blockDimX+2];
233         new_temperature = new double**[blockDimX+2];
234         for (i=0; i<blockDimX+2; i++) {
235           temperature[i] = new double*[blockDimY+2];
236           new_temperature[i] = new double*[blockDimY+2];
237           for(j=0; j<blockDimY+2; j++) {
238             temperature[i][j] = new double[blockDimZ+2];
239             new_temperature[i][j] = new double[blockDimZ+2];
240           }
241         }
242 #else
243         temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
244         new_temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
245 #endif
246
247         for(i=0; i<blockDimX+2; ++i) {
248           for(j=0; j<blockDimY+2; ++j) {
249             for(k=0; k<blockDimZ+2; ++k) {
250               temperature[index(i, j, k)] = 0.0;
251             }
252           } 
253         }
254
255         iterations = 0;
256         arrived_left = 0;
257         arrived_right = 0;
258         arrived_top = 0;
259         arrived_bottom = 0;
260         arrived_front = 0;
261         arrived_back = 0;
262         msgs = 0;
263         BC();
264     }
265
266     // Enforce some boundary conditions
267     void BC() {
268         // Heat left, top and front faces of each chare's block
269         for(int i=1; i<blockDimX+1; ++i)
270           for(int k=1; k<blockDimZ+1; ++k)
271             temperature[index(i, 1, k)] = 255.0;
272         for(int j=1; j<blockDimY+1; ++j)
273           for(int k=1; k<blockDimZ+1; ++k)
274             temperature[index(1, j, k)] = 255.0;
275         for(int i=1; i<blockDimX+1; ++i)
276           for(int j=1; j<blockDimY+1; ++j)
277             temperature[index(i, j, 1)] = 255.0;
278     }
279
280     // a necessary function which we ignore now
281     // if we were to use load balancing and migration
282     // this function might become useful
283     Jacobi(CkMigrateMessage* m) {}
284
285     ~Jacobi() { 
286 #if USE_3D_ARRAYS
287       for (int i=0; i<blockDimX+2; i++) {
288         for(int j=0; j<blockDimY+2; j++) {
289           delete [] temperature[i][j];
290           delete [] new_temperature[i][j];
291         }
292         delete [] temperature[i];
293         delete [] new_temperature[i];
294       }
295       delete [] temperature; 
296       delete [] new_temperature; 
297 #else
298       delete [] temperature; 
299       delete [] new_temperature; 
300 #endif
301     }
302
303     // Perform one iteration of work
304     // The first step is to send the local state to the neighbors
305     void begin_iteration(void) {
306
307         // Copy different faces into temporary arrays
308         double *left_face = new double[blockDimY*blockDimZ];
309         double *right_face = new double[blockDimY*blockDimZ];
310         double *top_face = new double[blockDimX*blockDimZ];
311         double *bottom_face = new double[blockDimX*blockDimZ];
312         double *front_face = new double[blockDimX*blockDimY];
313         double *back_face = new double[blockDimX*blockDimY];
314
315         for(int j=0; j<blockDimY; ++j) 
316           for(int k=0; k<blockDimZ; ++k) {
317             left_face[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
318             right_face[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
319         }
320
321         for(int i=0; i<blockDimX; ++i) 
322           for(int k=0; k<blockDimZ; ++k) {
323             top_face[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
324             bottom_face[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
325         }
326
327         for(int i=0; i<blockDimX; ++i) 
328           for(int j=0; j<blockDimY; ++j) {
329             front_face[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
330             back_face[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
331         }
332
333         // Send my left face
334         thisProxy(wrap_x(thisIndex.x-1), thisIndex.y, thisIndex.z).ghostsFromRight(blockDimY, blockDimZ, left_face);
335         // Send my right face
336         thisProxy(wrap_x(thisIndex.x+1), thisIndex.y, thisIndex.z).ghostsFromLeft(blockDimY, blockDimZ, right_face);
337         // Send my top face
338         thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z).ghostsFromBottom(blockDimX, blockDimZ, top_face);
339         // Send my bottom face
340         thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z).ghostsFromTop(blockDimX, blockDimZ, bottom_face);
341         // Send my front face
342         thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z-1)).ghostsFromFront(blockDimX, blockDimY, back_face);
343         // Send my back face
344         thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z+1)).ghostsFromBack(blockDimX, blockDimY, front_face);
345
346         delete [] back_face;
347         delete [] front_face;
348         delete [] bottom_face;
349         delete [] top_face; 
350         delete [] right_face;
351         delete [] left_face;
352     }
353
354     void ghostsFromRight(int height, int width, double ghost_values[]) {
355         for(int j=0; j<height; ++j) 
356           for(int k=0; k<width; ++k) {
357             temperature[index(blockDimX+1, j+1, k+1)] = ghost_values[k*height+j];
358         }
359         check_and_compute(RIGHT);
360     }
361
362     void ghostsFromLeft(int height, int width, double ghost_values[]) {
363         for(int j=0; j<height; ++j) 
364           for(int k=0; k<width; ++k) {
365             temperature[index(0, j+1, k+1)] = ghost_values[k*height+j];
366         }
367         check_and_compute(LEFT);
368     }
369
370     void ghostsFromBottom(int height, int width, double ghost_values[]) {
371         for(int i=0; i<height; ++i) 
372           for(int k=0; k<width; ++k) {
373             temperature[index(i+1, blockDimY+1, k+1)] = ghost_values[k*height+i];
374         }
375         check_and_compute(BOTTOM);
376     }
377
378     void ghostsFromTop(int height, int width, double ghost_values[]) {
379         for(int i=0; i<height; ++i) 
380           for(int k=0; k<width; ++k) {
381             temperature[index(i+1, 0, k+1)] = ghost_values[k*height+i];
382         }
383         check_and_compute(TOP);
384     }
385
386     void ghostsFromFront(int height, int width, double ghost_values[]) {
387         for(int i=0; i<height; ++i) 
388           for(int j=0; j<width; ++j) {
389             temperature[index(i+1, j+1, blockDimZ+1)] = ghost_values[j*height+i];
390         }
391         check_and_compute(FRONT);
392     }
393
394     void ghostsFromBack(int height, int width, double ghost_values[]) {
395         for(int i=0; i<height; ++i) 
396           for(int j=0; j<width; ++j) {
397             temperature[index(i+1, j+1, 0)] = ghost_values[j*height+i];
398         }
399         check_and_compute(BACK);
400     }
401
402     void check_and_compute(int direction) {
403         switch(direction) {
404           case LEFT:
405             arrived_left++;
406             msgs++;
407             break;
408           case RIGHT:
409             arrived_right++;
410             msgs++;
411             break;
412           case TOP:
413             arrived_top++;
414             msgs++;
415             break;
416           case BOTTOM:
417             arrived_bottom++;
418             msgs++;
419             break;
420           case FRONT:
421             arrived_front++;
422             msgs++;
423             break;
424           case BACK:
425             arrived_back++;
426             msgs++;
427             break;
428         }
429         if (arrived_left >=1 && arrived_right >=1 && arrived_top >=1 && arrived_bottom >=1 && arrived_front >= 1 && arrived_back >= 1) {
430           arrived_left--;
431           arrived_right--;
432           arrived_top--;
433           arrived_bottom--;
434           arrived_front--;
435           arrived_back--;
436           compute();
437           iterations++;
438           if (iterations == 1 || iterations == MAX_ITER)
439             contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
440           else
441             begin_iteration();
442         }
443     }
444
445     // Check to see if we have received all neighbor values yet
446     // If all neighbor values have been received, we update our values and proceed
447     void compute() {
448         // We must create a new array for these values because we don't want to 
449         // update any of the values in temperature[][] array until using them first.
450         // Other schemes could be used to accomplish this same problem. We just put
451         // the new values in a temporary array and write them to temperature[][] 
452         // after all of the new values are computed.
453
454 #pragma unroll    
455         for(int i=1; i<blockDimX+1; ++i) {
456           for(int j=1; j<blockDimY+1; ++j) {
457             for(int k=1; k<blockDimZ+1; ++k) {
458               // update my value based on the surrounding values
459               new_temperature[index(i, j, k)] = (temperature[index(i-1, j, k)] 
460                                               +  temperature[index(i+1, j, k)]
461                                               +  temperature[index(i, j-1, k)]
462                                               +  temperature[index(i, j+1, k)]
463                                               +  temperature[index(i, j, k-1)]
464                                               +  temperature[index(i, j, k+1)]
465                                               +  temperature[index(i, j, k)] ) * DIVIDEBY7;
466             }
467           }
468         }
469
470 #pragma unroll
471         for(int i=0;i<blockDimX+2;++i)
472           for(int j=0;j<blockDimY+2;++j)
473             for(int k=1; k<blockDimZ+1; ++k)
474                 temperature[index(i, j, k)] = new_temperature[index(i, j, k)];
475
476         // Enforce the boundary conditions again
477         BC();
478     }
479
480 };
481
482 /** \class JacobiMap
483  *
484  */
485
486 class JacobiMap : public CkArrayMap {
487   public:
488     int X, Y, Z;
489     int *mapping;
490
491     JacobiMap(int x, int y, int z, int tx, int ty, int tz) {
492       X = x; Y = y; Z = z;
493       mapping = new int[X*Y*Z];
494
495       // we are assuming that the no. of chares in each dimension is a 
496       // multiple of the torus dimension
497
498       TopoManager tmgr;
499       int dimX, dimY, dimZ, dimT;
500
501 #if USE_TOPOMAP
502       dimX = tmgr.getDimNX();
503       dimY = tmgr.getDimNY();
504       dimZ = tmgr.getDimNZ();
505       dimT = tmgr.getDimNT();
506 #elif USE_BLOCKMAP
507       dimX = tx;
508       dimY = ty;
509       dimZ = tz;
510       dimT = 1;
511 #endif
512
513 #if USE_RRMAP
514       if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d\n", x, y, z, numCharesPerPe, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ); 
515       int pe = 0;
516       for(int i=0; i<x; i++)
517         for(int j=0; j<y; j++)
518           for(int k=0; k<z; k++) {
519             if(pe == CkNumPes()) {
520               pe = 0;
521             }
522             mapping[i*Y*Z + j*Z + k] = pe;
523             pe++;
524           }
525 #elif USE_SMPMAP
526       if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d\n", x, y, z, numCharesPerPe, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ); 
527       int pe = -1;
528       x /= (numCharesPerPeX*SMPWAYX);
529       y /= (numCharesPerPeY*SMPWAYY);
530       z /= (numCharesPerPeZ*SMPWAYZ);
531
532       for(int i=0; i<x; i++)
533         for(int j=0; j<y; j++)
534           for(int k=0; k<z; k++)
535             for(int bi=i*SMPWAYX; bi<(i+1)*SMPWAYX; bi++)
536               for(int bj=j*SMPWAYY; bj<(j+1)*SMPWAYY; bj++)
537                 for(int bk=k*SMPWAYZ; bk<(k+1)*SMPWAYZ; bk++) {
538                   pe++;
539                   for(int ci=bi*numCharesPerPeX; ci<(bi+1)*numCharesPerPeX; ci++)
540                     for(int cj=bj*numCharesPerPeY; cj<(bj+1)*numCharesPerPeY; cj++)
541                       for(int ck=bk*numCharesPerPeZ; ck<(bk+1)*numCharesPerPeZ; ck++) {
542                         //if(CkMyPe()==0) CkPrintf("%d %d %d %d %d %d [%d]\n", i, j, k, ci, cj, ck, pe);
543                         mapping[ci*Y*Z + cj*Z + ck] = pe;
544                       }
545                 }
546 #else
547
548       // we are assuming that the no. of chares in each dimension is a 
549       // multiple of the torus dimension
550       int numCharesPerPe = X*Y*Z/CkNumPes();
551
552       int numCharesPerPeX = X / dimX;
553       int numCharesPerPeY = Y / dimY;
554       int numCharesPerPeZ = Z / dimZ;
555       int pe = 0, pes = CkNumPes();
556
557 #if USE_BLOCK_RNDMAP
558       int used[pes];
559       for(int i=0; i<pes; i++)
560         used[i] = 0;
561 #endif
562
563       if(dimT < 2) {    // one core per node
564         if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimX, dimY, dimZ, dimT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ); 
565         for(int i=0; i<dimX; i++)
566           for(int j=0; j<dimY; j++)
567             for(int k=0; k<dimZ; k++)
568             {
569 #if USE_BLOCK_RNDMAP
570               pe = myrand(pes); 
571               while(used[pe]!=0) {
572                 pe = myrand(pes); 
573               }
574               used[pe] = 1;
575 #endif
576
577               for(int ci=i*numCharesPerPeX; ci<(i+1)*numCharesPerPeX; ci++)
578                 for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
579                   for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
580 #if USE_TOPOMAP
581                     mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k);
582 #elif USE_BLOCKMAP
583                     mapping[ci*Y*Z + cj*Z + ck] = i + j*dimX + k*dimX*dimY;
584   #if USE_BLOCK_RNDMAP
585                     mapping[ci*Y*Z + cj*Z + ck] = pe;
586   #elif USE_BLOCK_RRMAP
587                     mapping[ci*Y*Z + cj*Z + ck] = pe;
588   #endif
589 #endif
590                   }
591 #if USE_BLOCK_RRMAP
592               pe++;
593 #endif
594             }
595       } else {          // multiple cores per node
596       // In this case, we split the chares in the X dimension among the
597       // cores on the same node. The strange thing I figured out is that
598       // doing this in the Z dimension is not as good.
599       numCharesPerPeX /= dimT;
600       if(CkMyPe()==0) CkPrintf("%d %d %d : %d %d %d %d : %d %d %d \n", x, y, z, dimX, dimY, dimZ, dimT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ); 
601       for(int i=0; i<dimX; i++)
602         for(int j=0; j<dimY; j++)
603           for(int k=0; k<dimZ; k++)
604             for(int l=0; l<dimT; l++)
605               for(int ci=(dimT*i+l)*numCharesPerPeX; ci<(dimT*i+l+1)*numCharesPerPeX; ci++)
606                 for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
607                   for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
608                     mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k, l);
609                   }
610       } // end of if
611 #endif
612
613       if(CkMyPe() == 0) CkPrintf("Map generated ... \n");
614     }
615
616     ~JacobiMap() { 
617       delete [] mapping;
618     }
619
620     int procNum(int, const CkArrayIndex &idx) {
621       int *index = (int *)idx.data();
622       return mapping[index[0]*Y*Z + index[1]*Z + index[2]]; 
623     }
624 };
625
626 #include "jacobi3d.def.h"