jacobi3d in sdag
[charm.git] / examples / bigsim / sdag / jacobi3d / jacobi3d.C
1 /*****************************************************************************
2  * $Source$
3  * $Author$
4  * $Date$
5  * $Revision$
6  *****************************************************************************/
7
8 /** \file jacobi3d.C
9  *  Author: Abhinav S Bhatele
10  *  Date Created: October 24th, 2007
11  *
12  *  This does a topological placement for a 3d jacobi.
13  *
14  *      
15  *            *****************
16  *         *               *  *
17  *   ^  *****************     *
18  *   |  *               *     *
19  *   |  *               *     *
20  *   |  *               *     *
21  *   Y  *               *     *
22  *   |  *               *     *
23  *   |  *               *     *
24  *   |  *               *  * 
25  *   ~  *****************    Z
26  *      <------ X ------> 
27  *
28  *   X: left, right --> wrap_x
29  *   Y: top, bottom --> wrap_y
30  *   Z: front, back --> wrap_z
31  */
32
33 #include "jacobi3d.decl.h"
34 #include "TopoManager.h"
35
36 // See README for documentation
37
38 /*readonly*/ CProxy_Main mainProxy;
39 /*readonly*/ int arrayDimX;
40 /*readonly*/ int arrayDimY;
41 /*readonly*/ int arrayDimZ;
42 /*readonly*/ int blockDimX;
43 /*readonly*/ int blockDimY;
44 /*readonly*/ int blockDimZ;
45 /*readonly*/ int torusDimX;
46 /*readonly*/ int torusDimY;
47 /*readonly*/ int torusDimZ;
48
49 // specify the number of worker chares in each dimension
50 /*readonly*/ int num_chare_x;
51 /*readonly*/ int num_chare_y;
52 /*readonly*/ int num_chare_z;
53
54 /*readonly*/ int globalBarrier;
55 /*readonly*/ int localBarrier;
56
57 static unsigned long next = 1;
58
59 int myrand(int numpes) {
60   next = next * 1103515245 + 12345;
61   return((unsigned)(next/65536) % numpes);
62 }
63
64 #define SMPWAYX                 1
65 #define SMPWAYY                 2
66 #define SMPWAYZ                 2
67 #define USE_TOPOMAP             0
68 #define USE_RRMAP               0
69 #define USE_BLOCKMAP            1
70 #define USE_BLOCK_RNDMAP        0
71 #define USE_BLOCK_RRMAP         1
72 #define USE_SMPMAP              0
73 #define USE_3D_ARRAYS           0
74
75 // We want to wrap entries around, and because mod operator % 
76 // sometimes misbehaves on negative values. -1 maps to the highest value.
77 #define wrap_x(a)       (((a)+num_chare_x)%num_chare_x)
78 #define wrap_y(a)       (((a)+num_chare_y)%num_chare_y)
79 #define wrap_z(a)       (((a)+num_chare_z)%num_chare_z)
80
81 #if USE_3D_ARRAYS
82 #define index(a, b, c)  a][b][c 
83 #else
84 #define index(a, b, c)  (a*(blockDimY+2)*(blockDimZ+2) + b*(blockDimZ+2) + c)
85 #endif
86
87 #define MAX_ITER                21
88 #define WARM_ITER               5
89 #define LEFT                    1
90 #define RIGHT                   2
91 #define TOP                     3
92 #define BOTTOM                  4
93 #define FRONT                   5
94 #define BACK                    6
95 #define DIVIDEBY7               0.14285714285714285714
96
97 double startTime;
98 double endTime;
99
100 /** \class ghostMsg
101  *
102  */
103 class ghostMsg: public CMessage_ghostMsg {
104   public:
105     int dir;
106     int height;
107     int width;
108     double* gh;
109
110     ghostMsg(int _d, int _h, int _w) : dir(_d), height(_h), width(_w) {
111     }
112 };
113
114 /** \class Main
115  *
116  */
117 class Main : public CBase_Main {
118   public:
119     CProxy_Jacobi array;
120     int iterations;
121
122     Main(CkArgMsg* m) {
123       if ( (m->argc != 3) && (m->argc != 7) && (m->argc != 10) ) {
124         CkPrintf("%s [array_size] [block_size]\n", m->argv[0]);
125         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]);
126         CkAbort("Abort");
127       }
128
129       // set iteration counter to zero
130       iterations = 0;
131
132       // store the main proxy
133       mainProxy = thisProxy;
134         
135       if(m->argc == 3) {
136         arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
137         blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]); 
138       }
139       else if (m->argc == 7) {
140         arrayDimX = atoi(m->argv[1]);
141         arrayDimY = atoi(m->argv[2]);
142         arrayDimZ = atoi(m->argv[3]);
143         blockDimX = atoi(m->argv[4]); 
144         blockDimY = atoi(m->argv[5]); 
145         blockDimZ = atoi(m->argv[6]);
146       } else {
147         arrayDimX = atoi(m->argv[1]);
148         arrayDimY = atoi(m->argv[2]);
149         arrayDimZ = atoi(m->argv[3]);
150         blockDimX = atoi(m->argv[4]); 
151         blockDimY = atoi(m->argv[5]); 
152         blockDimZ = atoi(m->argv[6]);
153         torusDimX = atoi(m->argv[7]);
154         torusDimY = atoi(m->argv[8]);
155         torusDimZ = atoi(m->argv[9]);
156       }
157
158       if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
159         CkAbort("array_size_X % block_size_X != 0!");
160       if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
161         CkAbort("array_size_Y % block_size_Y != 0!");
162       if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
163         CkAbort("array_size_Z % block_size_Z != 0!");
164
165       num_chare_x = arrayDimX / blockDimX;
166       num_chare_y = arrayDimY / blockDimY;
167       num_chare_z = arrayDimZ / blockDimZ;
168
169       // print info
170       CkPrintf("Running Jacobi on %d processors with (%d, %d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y, num_chare_z);
171       CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
172       CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
173
174       // Create new array of worker chares
175 #if USE_TOPOMAP || USE_RRMAP || USE_BLOCKMAP || USE_BLOCK_RRMAP || USE_BLOCK_RNDMAP || USE_SMPMAP
176        CkPrintf("Topology Mapping is being done ... %d %d %d %d\n", USE_TOPOMAP, USE_RRMAP, USE_BLOCKMAP, USE_SMPMAP);
177       CProxy_JacobiMap map = CProxy_JacobiMap::ckNew(num_chare_x, num_chare_y, num_chare_z, torusDimX, torusDimY, torusDimZ);
178       CkArrayOptions opts(num_chare_x, num_chare_y, num_chare_z);
179       opts.setMap(map);
180       array = CProxy_Jacobi::ckNew(opts);
181 #else
182       array = CProxy_Jacobi::ckNew(num_chare_x, num_chare_y, num_chare_z);
183 #endif
184
185       TopoManager tmgr;
186       CkArray *jarr = array.ckLocalBranch();
187       int jmap[num_chare_x][num_chare_y][num_chare_z];
188
189       int hops=0, p;
190       for(int i=0; i<num_chare_x; i++)
191         for(int j=0; j<num_chare_y; j++)
192           for(int k=0; k<num_chare_z; k++) {
193             jmap[i][j][k] = jarr->procNum(CkArrayIndex3D(i, j, k));
194           }
195
196       for(int i=0; i<num_chare_x; i++)
197         for(int j=0; j<num_chare_y; j++)
198           for(int k=0; k<num_chare_z; k++) {
199             p = jmap[i][j][k];
200             hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i+1)][j][k]);
201             hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i-1)][j][k]);
202             hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j+1)][k]);
203             hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j-1)][k]);
204             hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k+1)]);
205             hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k-1)]);
206           }
207       CkPrintf("Total Hops: %d\n", hops);
208
209       //Start the computation
210       array.doStep();
211     }
212
213     // Each worker reports back to here when it completes an iteration
214     void report() {
215       if (iterations == 0) {
216         iterations++;
217         startTime = CmiWallTimer();
218         array.doStep();
219       }
220       else {
221         endTime = CmiWallTimer();
222         CkPrintf("TIME : %f\n", (endTime - startTime)/(MAX_ITER-WARM_ITER-1));
223         CkExit();
224       }
225     }
226 };
227
228 /** \class Jacobi
229  *
230  */
231
232 class Jacobi: public CBase_Jacobi {
233   Jacobi_SDAG_CODE
234
235   public:
236     int arrived_left;
237     int arrived_right;
238     int arrived_top;
239     int arrived_bottom;
240     int arrived_front;
241     int arrived_back;
242     int iterations;
243     int imsg;
244
245 #if USE_3D_ARRAYS
246     double ***temperature;
247     double ***new_temperature;
248 #else
249     double *temperature;
250     double *new_temperature;
251 #endif
252
253     // Constructor, initialize values
254     Jacobi() {
255       __sdag_init();
256
257       int i, j, k;
258       // allocate a three dimensional array
259 #if USE_3D_ARRAYS
260       temperature = new double**[blockDimX+2];
261       new_temperature = new double**[blockDimX+2];
262       for (i=0; i<blockDimX+2; i++) {
263         temperature[i] = new double*[blockDimY+2];
264         new_temperature[i] = new double*[blockDimY+2];
265         for(j=0; j<blockDimY+2; j++) {
266           temperature[i][j] = new double[blockDimZ+2];
267           new_temperature[i][j] = new double[blockDimZ+2];
268         }
269       }
270 #else
271       temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
272       new_temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
273 #endif
274
275       for(i=0; i<blockDimX+2; ++i) {
276         for(j=0; j<blockDimY+2; ++j) {
277           for(k=0; k<blockDimZ+2; ++k) {
278             temperature[index(i, j, k)] = 0.0;
279           }
280         } 
281       }
282
283       iterations = 0;
284       arrived_left = 0;
285       arrived_right = 0;
286       arrived_top = 0;
287       arrived_bottom = 0;
288       arrived_front = 0;
289       arrived_back = 0;
290       imsg = 0;
291       constrainBC();
292     }
293
294     // Enforce some boundary conditions
295     void constrainBC() {
296         // Heat left, top and front faces of each chare's block
297         for(int i=1; i<blockDimX+1; ++i)
298           for(int k=1; k<blockDimZ+1; ++k)
299             temperature[index(i, 1, k)] = 255.0;
300         for(int j=1; j<blockDimY+1; ++j)
301           for(int k=1; k<blockDimZ+1; ++k)
302             temperature[index(1, j, k)] = 255.0;
303         for(int i=1; i<blockDimX+1; ++i)
304           for(int j=1; j<blockDimY+1; ++j)
305             temperature[index(i, j, 1)] = 255.0;
306     }
307
308     // a necessary function which we ignore now
309     // if we were to use load balancing and migration
310     // this function might become useful
311     Jacobi(CkMigrateMessage* m) {}
312
313     ~Jacobi() { 
314 #if USE_3D_ARRAYS
315       for (int i=0; i<blockDimX+2; i++) {
316         for(int j=0; j<blockDimY+2; j++) {
317           delete [] temperature[i][j];
318           delete [] new_temperature[i][j];
319         }
320         delete [] temperature[i];
321         delete [] new_temperature[i];
322       }
323       delete [] temperature; 
324       delete [] new_temperature; 
325 #else
326       delete [] temperature; 
327       delete [] new_temperature; 
328 #endif
329     }
330
331     // Perform one iteration of work
332     // The first step is to send the local state to the neighbors
333     void begin_iteration(void) {
334       // Copy different faces into messages
335       ghostMsg *leftMsg = new (blockDimY*blockDimZ) ghostMsg(RIGHT, blockDimY, blockDimZ);
336       ghostMsg *rightMsg = new (blockDimY*blockDimZ) ghostMsg(LEFT, blockDimY, blockDimZ);
337       ghostMsg *topMsg = new (blockDimX*blockDimZ) ghostMsg(BOTTOM, blockDimX, blockDimZ);
338       ghostMsg *bottomMsg = new (blockDimX*blockDimZ) ghostMsg(TOP, blockDimX, blockDimZ);
339       ghostMsg *frontMsg = new (blockDimX*blockDimY) ghostMsg(BACK, blockDimX, blockDimY);
340       ghostMsg *backMsg = new (blockDimX*blockDimY) ghostMsg(FRONT, blockDimX, blockDimY);
341
342
343       for(int j=0; j<blockDimY; ++j) 
344         for(int k=0; k<blockDimZ; ++k) {
345           leftMsg->gh[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
346           rightMsg->gh[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];
347       }
348
349       for(int i=0; i<blockDimX; ++i) 
350         for(int k=0; k<blockDimZ; ++k) {
351           topMsg->gh[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
352           bottomMsg->gh[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
353       }
354
355       for(int i=0; i<blockDimX; ++i) 
356         for(int j=0; j<blockDimY; ++j) {
357           frontMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
358           backMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
359       }
360
361       // Send my left face
362       thisProxy(wrap_x(thisIndex.x-1), thisIndex.y, thisIndex.z).receiveGhosts(leftMsg);
363       // Send my right face
364       thisProxy(wrap_x(thisIndex.x+1), thisIndex.y, thisIndex.z).receiveGhosts(rightMsg);
365       // Send my top face
366       thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z).receiveGhosts(topMsg);
367       // Send my bottom face
368       thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z).receiveGhosts(bottomMsg);
369       // Send my front face
370       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z-1)).receiveGhosts(frontMsg);
371       // Send my back face
372       thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z-1)).receiveGhosts(backMsg);
373     }
374
375     void processGhosts(ghostMsg *gmsg) {
376       int height = gmsg->height;
377       int width = gmsg->width;
378
379       switch(gmsg->dir) {
380         case LEFT:
381           arrived_left++;
382           for(int j=0; j<height; ++j) 
383             for(int k=0; k<width; ++k) {
384               temperature[index(0, j+1, k+1)] = gmsg->gh[k*height+j];
385           }
386           break;
387         case RIGHT:
388           arrived_right++;
389           for(int j=0; j<height; ++j) 
390             for(int k=0; k<width; ++k) {
391               temperature[index(blockDimX+1, j+1, k+1)] = gmsg->gh[k*height+j];
392           }
393           break;
394         case TOP:
395           arrived_top++;
396           for(int i=0; i<height; ++i) 
397             for(int k=0; k<width; ++k) {
398               temperature[index(i+1, 0, k+1)] = gmsg->gh[k*height+i];
399           }
400           break;
401         case BOTTOM:
402           arrived_bottom++;
403           for(int i=0; i<height; ++i) 
404             for(int k=0; k<width; ++k) {
405               temperature[index(i+1, blockDimY+1, k+1)] = gmsg->gh[k*height+i];
406           }
407           break;
408         case FRONT:
409           arrived_front++;
410           for(int i=0; i<height; ++i) 
411             for(int j=0; j<width; ++j) {
412               temperature[index(i+1, j+1, blockDimZ+1)] = gmsg->gh[j*height+i];
413           }
414           break;
415         case BACK:
416           arrived_back++;
417           for(int i=0; i<height; ++i) 
418             for(int j=0; j<width; ++j) {
419               temperature[index(i+1, j+1, 0)] = gmsg->gh[j*height+i];
420           }
421           break;
422         default:
423           CkAbort("ERROR\n");
424       }
425     }
426
427
428     void check_and_compute() {
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_kernel();
437         iterations++;
438         if (iterations < WARM_ITER || iterations == MAX_ITER)
439           contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
440         else
441           doStep();
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_kernel() {
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       constrainBC();
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"