fix gnu compilation error on Blue Gene/Q
[charm.git] / tests / charm++ / jacobi3d-gausssiedel / 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: June 01st, 2009
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 "jacobi3d.h"
35 #include "TopoManager.h"
36
37 #include "CkLoopAPI.h"
38 CProxy_FuncCkLoop ckLoopProxy;
39 #ifdef JACOBI_OPENMP
40 #include <omp.h>
41 #endif
42
43 //#define PRINT_DEBUG 1
44 #define  LOW_VALUE 0 
45 #define  HIGH_VALUE 255 
46 //#define JACOBI  1 
47
48 int GAUSS_ITER=1; 
49 #define THRESHOLD 0.001 
50 // See README for documentation
51
52 /*readonly*/ CProxy_Main mainProxy;
53 /*readonly*/ int arrayDimX;
54 /*readonly*/ int arrayDimY;
55 /*readonly*/ int blockDimX;
56 /*readonly*/ int arrayDimZ;
57 /*readonly*/ int blockDimY;
58 /*readonly*/ int blockDimZ;
59
60 // specify the number of worker chares in each dimension
61 /*readonly*/ int num_chare_x;
62 /*readonly*/ int num_chare_y;
63 /*readonly*/ int num_chare_z;
64
65 /*readonly*/ int globalBarrier;
66
67 static unsigned long next = 1;
68
69 int myrand(int numpes) {
70   next = next * 1103515245 + 12345;
71   return((unsigned)(next/65536) % numpes);
72 }
73
74 // We want to wrap entries around, and because mod operator % 
75 // sometimes misbehaves on negative values. -1 maps to the highest value.
76 #define wrap_x(a)       (((a)+num_chare_x)%num_chare_x)
77 #define wrap_y(a)       (((a)+num_chare_y)%num_chare_y)
78 #define wrap_z(a)       (((a)+num_chare_z)%num_chare_z)
79
80 #define index(a, b, c)  ( (a)*(blockDimY+2)*(blockDimZ+2) + (b)*(blockDimZ+2) + (c) )
81
82 #define  START_ITER     10
83 #define   END_ITER      20
84 #define PRINT_FREQ     100 
85 #define CKP_FREQ                100
86 #define MAX_ITER         10     
87 #define WARM_ITER               5
88 #define LEFT                    1
89 #define RIGHT                   2
90 #define TOP                     3
91 #define BOTTOM                  4
92 #define FRONT                   5
93 #define BACK                    6
94 #define DIVIDEBY7               0.14285714285714285714
95
96 double startTime;
97 double endTime;
98
99 /** \class ghostMsg
100  *
101  */
102 class ghostMsg: public CMessage_ghostMsg {
103   public:
104     int dir;
105     int height;
106     int width;
107     double* gh;
108
109     ghostMsg(int _d, int _h, int _w) : dir(_d), height(_h), width(_w) {
110     }
111 };
112
113 /** \class Main
114  *
115  */
116 class Main : public CBase_Main {
117   public:
118     CProxy_Jacobi array;
119     CProxy_TraceControl _traceControl;
120     int iterations;
121     Main(CkArgMsg* m) {
122       if ( (m->argc != 3) && (m->argc < 7) ) {
123         CkPrintf("%s [array_size] [block_size] gauss-iteration\n", m->argv[0]);
124         CkPrintf("OR %s [array_size_X] [array_size_Y] [array_size_Z] [block_size_X] [block_size_Y] [block_size_Z] [gauss-iter]\n", m->argv[0]);
125         CkAbort("Abort");
126       }
127
128       // set iteration counter to zero
129       iterations = 0;
130
131       // store the main proxy
132       mainProxy = thisProxy;
133         
134       if(m->argc == 3 ||  m->argc == 4) {
135           arrayDimX = arrayDimY = arrayDimZ = atoi(m->argv[1]);
136           blockDimX = blockDimY = blockDimZ = atoi(m->argv[2]); 
137           if(m->argc == 4 )
138               GAUSS_ITER = atoi(m->argv[3]);
139       }
140       else if (m->argc >= 7) {
141           arrayDimX = atoi(m->argv[1]);
142           arrayDimY = atoi(m->argv[2]);
143           arrayDimZ = atoi(m->argv[3]);
144           blockDimX = atoi(m->argv[4]); 
145           blockDimY = atoi(m->argv[5]); 
146           blockDimZ = atoi(m->argv[6]);
147           if(m->argc > 7)
148               GAUSS_ITER = atoi(m->argv[7]);
149       } 
150
151       if (arrayDimX < blockDimX || arrayDimX % blockDimX != 0)
152         CkAbort("array_size_X % block_size_X != 0!");
153       if (arrayDimY < blockDimY || arrayDimY % blockDimY != 0)
154         CkAbort("array_size_Y % block_size_Y != 0!");
155       if (arrayDimZ < blockDimZ || arrayDimZ % blockDimZ != 0)
156         CkAbort("array_size_Z % block_size_Z != 0!");
157
158       num_chare_x = arrayDimX / blockDimX;
159       num_chare_y = arrayDimY / blockDimY;
160       num_chare_z = arrayDimZ / blockDimZ;
161
162       // print info
163       CkPrintf("\nSTENCIL COMPUTATION WITH NO BARRIERS\n");
164       CkPrintf("Running Jacobi on %d processors with (%d, %d, %d) chares\n", CkNumPes(), num_chare_x, num_chare_y, num_chare_z);
165       CkPrintf("Array Dimensions: %d %d %d\n", arrayDimX, arrayDimY, arrayDimZ);
166       CkPrintf("Block Dimensions: %d %d %d\n", blockDimX, blockDimY, blockDimZ);
167
168       // Create new array of worker chares
169 #if USE_TOPOMAP
170       CProxy_JacobiMap map = CProxy_JacobiMap::ckNew(num_chare_x, num_chare_y, num_chare_z);
171       CkPrintf("Topology Mapping is being done ... \n");
172       CkArrayOptions opts(num_chare_x, num_chare_y, num_chare_z);
173       opts.setMap(map);
174       array = CProxy_Jacobi::ckNew(opts);
175 #else
176 #if CKLOOP
177       CProxy_JacobiNodeMap map = CProxy_JacobiNodeMap::ckNew();
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 #endif
185
186       TopoManager tmgr;
187       CkArray *jarr = array.ckLocalBranch();
188       int jmap[num_chare_x][num_chare_y][num_chare_z];
189
190       int hops=0, p;
191       for(int i=0; i<num_chare_x; i++)
192         for(int j=0; j<num_chare_y; j++)
193           for(int k=0; k<num_chare_z; k++) {
194             jmap[i][j][k] = jarr->procNum(CkArrayIndex3D(i, j, k));
195           }
196
197       for(int i=0; i<num_chare_x; i++)
198         for(int j=0; j<num_chare_y; j++)
199           for(int k=0; k<num_chare_z; k++) {
200             p = jmap[i][j][k];
201             hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i+1)][j][k]);
202             hops += tmgr.getHopsBetweenRanks(p, jmap[wrap_x(i-1)][j][k]);
203             hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j+1)][k]);
204             hops += tmgr.getHopsBetweenRanks(p, jmap[i][wrap_y(j-1)][k]);
205             hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k+1)]);
206             hops += tmgr.getHopsBetweenRanks(p, jmap[i][j][wrap_z(k-1)]);
207           }
208       CkPrintf("Total Hops: %d\n", hops);
209
210       _traceControl = CProxy_TraceControl::ckNew();
211 #if CKLOOP
212       ckLoopProxy = CkLoop_Init(CkMyNodeSize());
213 #endif
214       start();
215     }
216
217     void start() {
218       array.doStep();
219     }
220
221     void converge(double maxDiff)
222     {
223         iterations++;
224         if(iterations == START_ITER)
225             _traceControl.startTrace();
226         if(iterations == END_ITER)
227             _traceControl.endTrace();
228
229         if(iterations == WARM_ITER)
230             startTime = CmiWallTimer();                
231         //if(maxDiff <= THRESHOLD || iterations > 20)
232         if(maxDiff <= THRESHOLD )
233         {
234             endTime = CmiWallTimer();
235             CkPrintf("Completed:<x,y,z>chares  msgbytes iterations, times, timeperstep(ms)\n benchmark[%d:%d] [ %d %d %d ] \t\t %d \t\t %d \t\t %d\t\t %.3f \t\t %.3f \n", CkNumPes(), CkNumNodes(), num_chare_x, num_chare_y, num_chare_z, num_chare_x *num_chare_y * num_chare_z, sizeof(double) * blockDimX*blockDimY, iterations*GAUSS_ITER, endTime - startTime, (endTime - startTime)/(iterations-WARM_ITER)/GAUSS_ITER*1000);
236             //array.print();
237             CkExit();
238         }else
239         {
240             if(iterations % PRINT_FREQ== 0)
241                 CkPrintf("iteration %d  diff=%e\n", iterations, maxDiff); 
242                         array.doStep();
243         }
244     }
245     // Each worker reports back to here when it completes an iteration
246         void report() {
247         iterations += CKP_FREQ;
248         if (iterations < MAX_ITER) {
249 #ifdef CMK_MEM_CHECKPOINT
250             CkCallback cb (CkIndex_Jacobi::doStep(), array);
251             CkStartMemCheckpoint(cb);
252 #else
253                         array.doStep();
254 #endif
255         } else {
256             CkPrintf("Completed %d iterations\n", MAX_ITER-1);
257             endTime = CmiWallTimer();
258             CkPrintf("Time elapsed per iteration: %f\n", (endTime - startTime)/(MAX_ITER-1-WARM_ITER));
259             CkExit();
260         }
261             CkExit();
262     }
263 };
264
265 /** \class Jacobi
266  *
267  */
268
269 extern "C" void doCalc(int first,int last, void *result, int paramNum, void * param)
270 {
271 #if JACOBI
272     int i;
273     double maxDiff=0;
274     double tmp, update=0;
275     maxDiff = 0;
276
277     Jacobi  *obj = (Jacobi*)param;
278     int x_c = arrayDimX/2/blockDimX;
279     int y_c = arrayDimY/2/blockDimY;
280     int z_c = arrayDimZ/2/blockDimZ;
281     int x_p = (arrayDimX/2)%blockDimX+1;
282     int y_p = (arrayDimY/2)%blockDimY+1;
283     int z_p = (arrayDimZ/2)%blockDimZ+1;
284     //CkPrintf(" ckloop on %d  first last %d %d \n", CkMyPe(), first, last);  
285     for(i=first; i<=last; ++i) {
286           for(int j=1; j<blockDimY+1; ++j) {
287               for(int k=1; k<blockDimZ+1; ++k) {
288                   if(obj->thisIndex_x == x_c && y_c == obj->thisIndex_y && z_c == obj->thisIndex_z && i==x_p && j==y_p && k==z_p)
289                       continue;
290                   update = 0;//temperature[index(i, j, k)];
291                   update += obj->temperature[index(i+1, j, k)];
292                   update += obj->temperature[index(i-1, j, k)];
293                   update += obj->temperature[index(i, j+1, k)];
294                   update += obj->temperature[index(i, j-1, k)];
295                   update += obj->temperature[index(i, j, k+1)];
296                   update += obj->temperature[index(i, j, k-1)];
297                   update /= 6;
298                  
299                   double diff = obj->temperature[index(i, j, k)] - update;
300                   if(diff<0) diff = -1*diff;
301                   if(diff>maxDiff) maxDiff = diff;
302                   obj->new_temperature[index(i, j, k)] = update;
303               }
304           }
305         }
306         *((double*)result) = maxDiff;
307 #else
308 #endif 
309 }
310 extern "C" void doUpdate(int first,int last, void *result, int paramNum, void *param){
311 #if JACOBI
312         int i;
313         Jacobi  *obj = (Jacobi*)param;
314         int x_c = arrayDimX/2/blockDimX;
315         int y_c = arrayDimY/2/blockDimY;
316         int z_c = arrayDimZ/2/blockDimZ;
317         int x_p = (arrayDimX/2)%blockDimX+1;
318         int y_p = (arrayDimY/2)%blockDimY+1;
319         int z_p = (arrayDimZ/2)%blockDimZ+1;
320         //CkPrintf("\n {%d:%d:%d  %d:%d:%d  %d:%d:%d }", obj->thisIndex_x, obj->thisIndex_y, obj->thisIndex_z, x_c, y_c, z_c, x_p, y_p, z_p);
321         for(i=first; i<=last; ++i) {
322             for(int j=1; j<blockDimY+1; ++j) {
323                 for(int k=1; k<blockDimZ+1; ++k) {
324                     if(!(obj->thisIndex_x == x_c && y_c == obj->thisIndex_y && z_c == obj->thisIndex_z && i==x_p && j==y_p && k==z_p))
325                         obj->temperature[index(i, j, k)] = obj->new_temperature[index(i, j, k)];
326                     //CkPrintf("[ %d:%d:%d  .%2f ]", i, j, k, obj->temperature[index(i, j, k)]); 
327                 }
328             }
329         }
330        //CkPrintf("\n\n");
331 #else
332 #endif
333 }
334
335
336 Jacobi::Jacobi() {
337       // This call is an anachronism - the need to call __sdag_init() has been
338       // removed. We still call it here to test backward compatibility.
339       __sdag_init();
340
341       int i, j, k;
342       // allocate a three dimensional array
343       temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
344 #if  JACOBI
345       new_temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
346 #endif
347       for(i=0; i<blockDimX+2; ++i) {
348           for(j=0; j<blockDimY+2; ++j) {
349               for(k=0; k<blockDimZ+2; ++k) {
350                   temperature[index(i, j, k)] = 1.0;
351 #if  JACOBI
352                   new_temperature[index(i, j, k)] = 1.0;
353 #endif
354               }
355           } 
356       }
357       thisIndex_x = thisIndex.x;
358       thisIndex_y = thisIndex.y;
359       thisIndex_z = thisIndex.z;
360       iterations = 0;
361       imsg = 0;
362       constrainBC();
363
364       usesAtSync = CmiTrue;
365       neighbors = 6;
366       if(thisIndex.x == 0) 
367           neighbors--;
368       if( thisIndex.x== num_chare_x-1)
369           neighbors--;
370       if(thisIndex.y == 0) 
371           neighbors--;
372       if( thisIndex.y== num_chare_y-1)
373           neighbors--;
374       if(thisIndex.z == 0) 
375           neighbors--;
376       if( thisIndex.z== num_chare_z-1)
377           neighbors--;
378        // CkPrintf("neighbor = %d \n", neighbors);
379     }
380 void Jacobi::pup(PUP::er &p){
381                 
382                 // calling parent's pup
383                 CBase_Jacobi::pup(p);
384                 __sdag_pup(p);
385                 // pupping properties of this class
386                 p | iterations;
387                 p | imsg;
388                 // if unpacking, allocate the memory space
389                 if(p.isUnpacking()){
390                         temperature = new double[(blockDimX+2) * (blockDimY+2) * (blockDimZ+2)];
391                 }
392                 // pupping the arrays
393                 p((char *)temperature, (blockDimX+2) * (blockDimY+2) * (blockDimZ+2) * sizeof(double));
394         }
395     // Send ghost faces to the six neighbors
396     void Jacobi::begin_iteration(void) {
397       iterations++;
398
399       // Copy different faces into messages
400        ghostMsg *leftMsg ;
401        ghostMsg *rightMsg; 
402        ghostMsg *topMsg ;
403        ghostMsg *bottomMsg;
404        ghostMsg *frontMsg ;
405        ghostMsg *backMsg ;
406      
407        
408       if(thisIndex.x-1 >= 0)
409       {
410           leftMsg = new (blockDimY*blockDimZ) ghostMsg(RIGHT, blockDimY, blockDimZ);
411           for(int j=0; j<blockDimY; ++j) 
412               for(int k=0; k<blockDimZ; ++k) {
413                   leftMsg->gh[k*blockDimY+j] = temperature[index(1, j+1, k+1)];
414               }
415           CkSetRefNum(leftMsg, iterations);
416           thisProxy(wrap_x(thisIndex.x-1), thisIndex.y, thisIndex.z).receiveGhosts(leftMsg);
417       }
418
419       if(thisIndex.x+1 <num_chare_x) 
420       {
421           rightMsg = new (blockDimY*blockDimZ) ghostMsg(LEFT, blockDimY, blockDimZ);
422           for(int j=0; j<blockDimY; ++j) 
423               for(int k=0; k<blockDimZ; ++k) {
424                   rightMsg->gh[k*blockDimY+j] = temperature[index(blockDimX, j+1, k+1)];}
425           CkSetRefNum(rightMsg, iterations);
426           thisProxy(wrap_x(thisIndex.x+1), thisIndex.y, thisIndex.z).receiveGhosts(rightMsg);
427       }
428
429       if(thisIndex.y-1 >= 0)
430       {
431           topMsg = new (blockDimX*blockDimZ) ghostMsg(BOTTOM, blockDimX, blockDimZ);
432           for(int i=0; i<blockDimX; ++i) 
433               for(int k=0; k<blockDimZ; ++k) {
434                   topMsg->gh[k*blockDimX+i] = temperature[index(i+1, 1, k+1)];
435               }
436
437           CkSetRefNum(topMsg, iterations);
438           thisProxy(thisIndex.x, wrap_y(thisIndex.y-1), thisIndex.z).receiveGhosts(topMsg);
439       }
440
441       if(thisIndex.y+1 <num_chare_y)
442       {
443           bottomMsg = new (blockDimX*blockDimZ) ghostMsg(TOP, blockDimX, blockDimZ);
444           for(int i=0; i<blockDimX; ++i) 
445           for(int k=0; k<blockDimZ; ++k) {
446               bottomMsg->gh[k*blockDimX+i] = temperature[index(i+1, blockDimY, k+1)];
447           }
448           CkSetRefNum(bottomMsg, iterations);
449           thisProxy(thisIndex.x, wrap_y(thisIndex.y+1), thisIndex.z).receiveGhosts(bottomMsg);
450       }
451
452       if(thisIndex.z-1 >= 0)
453       {
454           frontMsg = new (blockDimX*blockDimY) ghostMsg(BACK, blockDimX, blockDimY);
455           for(int i=0; i<blockDimX; ++i) 
456               for(int j=0; j<blockDimY; ++j) {
457                   frontMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, 1)];
458               }
459           CkSetRefNum(frontMsg, iterations);
460           thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z-1)).receiveGhosts(frontMsg);
461       }
462       
463       if(thisIndex.z+1 <num_chare_z)
464       {
465           backMsg = new (blockDimX*blockDimY) ghostMsg(FRONT, blockDimX, blockDimY);
466           for(int i=0; i<blockDimX; ++i) 
467               for(int j=0; j<blockDimY; ++j) {
468                   backMsg->gh[j*blockDimX+i] = temperature[index(i+1, j+1, blockDimZ)];
469               }
470           CkSetRefNum(backMsg, iterations);
471           thisProxy(thisIndex.x, thisIndex.y, wrap_z(thisIndex.z+1)).receiveGhosts(backMsg);
472       }
473     }
474
475     void Jacobi::processGhosts(ghostMsg *gmsg) {
476       int height = gmsg->height;
477       int width = gmsg->width;
478       switch(gmsg->dir) {
479       case LEFT:
480           for(int j=0; j<height; ++j) 
481               for(int k=0; k<width; ++k) {
482                   temperature[index(0, j+1, k+1)] = gmsg->gh[k*height+j];
483               }
484           break;
485       case RIGHT:
486           for(int j=0; j<height; ++j) 
487               for(int k=0; k<width; ++k) {
488                   temperature[index(blockDimX+1, j+1, k+1)] = gmsg->gh[k*height+j];
489               }
490           break;
491       case TOP:
492           for(int i=0; i<height; ++i) 
493               for(int k=0; k<width; ++k) {
494                   temperature[index(i+1, 0, k+1)] = gmsg->gh[k*height+i];
495               }
496           break;
497       case BOTTOM:
498           for(int i=0; i<height; ++i) 
499               for(int k=0; k<width; ++k) {
500                   temperature[index(i+1, blockDimY+1, k+1)] = gmsg->gh[k*height+i];
501               }
502           break;
503       case FRONT:
504           for(int i=0; i<height; ++i) 
505               for(int j=0; j<width; ++j) {
506                   temperature[index(i+1, j+1,0)] = gmsg->gh[j*height+i];
507               }
508           break;
509       case BACK:
510           for(int i=0; i<height; ++i) 
511               for(int j=0; j<width; ++j) {
512                   temperature[index(i+1, j+1, blockDimZ+1)] = gmsg->gh[j*height+i];
513               }
514           break;
515       default:
516           CkAbort("ERROR\n");
517       }
518
519       delete gmsg;
520     }
521
522
523         void Jacobi::check_and_compute() {
524         double maxDiff;
525 #if CKLOOP
526         CkLoop_Parallelize(doCalc, 1,  (void*)this, 4*CkMyNodeSize(), 1, blockDimX, 1, &maxDiff, CKLOOP_DOUBLE_MAX);
527         CkLoop_Parallelize(doUpdate, 1,  (void*)this, CkMyNodeSize(), 1, blockDimX, 1);
528 #else
529                 maxDiff = compute_kernel();
530 #endif
531                 // calculate error
532                 // not being done right now since we are doing a fixed no. of iterations
533         //CkPrintf(" iteration %d [%d,%d,%d] max=%e\n", iterations, thisIndex.x, thisIndex.y, thisIndex.z, maxDiff);
534                 //constrainBC();
535                 if (iterations % CKP_FREQ == 0 || iterations > MAX_ITER){
536 #ifdef CMK_MEM_CHECKPOINT
537                         contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
538 #elif CMK_MESSAGE_LOGGING
539                         if(iterations > MAX_ITER)
540                                 contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
541                         else
542                                 AtSync();
543 #else
544             CkCallback cb(CkReductionTarget(Main, converge), mainProxy);
545             contribute(sizeof(double), &maxDiff, CkReduction::max_double, cb);
546 #endif
547         } else {
548             CkCallback cb(CkReductionTarget(Main, converge), mainProxy);
549             contribute(sizeof(double), &maxDiff, CkReduction::max_double, cb);
550                         //doStep();
551         }
552     }
553
554         void Jacobi::ResumeFromSync(){
555                 doStep();
556         }
557
558
559 #if JACOBI
560     double Jacobi::compute_kernel() {     //Gauss-Siedal compute
561         int i;
562         double maxDiff=0;
563         double tmp, update=0;
564         maxDiff = 0;
565
566         int x_c = arrayDimX/2/blockDimX;
567         int y_c = arrayDimY/2/blockDimY;
568         int z_c = arrayDimZ/2/blockDimZ;
569         int x_p = (arrayDimX/2)%blockDimX+1;
570         int y_p = (arrayDimY/2)%blockDimY+1;
571         int z_p = (arrayDimZ/2)%blockDimZ+1;
572   #pragma omp parallel for schedule(dynamic) 
573         for(i=1; i<blockDimX+1; ++i) {
574 #ifdef JACOBI_OPENMP
575           printf("[%d] did  %d iteration out of %d \n", omp_get_thread_num(), i, blockDimX+1); 
576 #endif
577           for(int j=1; j<blockDimY+1; ++j) {
578               for(int k=1; k<blockDimZ+1; ++k) {
579         
580                   if(thisIndex.x == x_c && y_c == thisIndex.y && z_c == thisIndex.z && i==x_p && j==y_p && k==z_p)
581                       continue;
582                   update = 0;//temperature[index(i, j, k)];
583                   update += temperature[index(i+1, j, k)];
584                   update += temperature[index(i-1, j, k)];
585                   update += temperature[index(i, j+1, k)];
586                   update += temperature[index(i, j-1, k)];
587                   update += temperature[index(i, j, k+1)];
588                   update += temperature[index(i, j, k-1)];
589                   update /= 6;
590                  
591                   double diff = temperature[index(i, j, k)] - update;
592                   if(diff<0) diff = -1*diff;
593                   if(diff>maxDiff) maxDiff = diff;
594                   new_temperature[index(i, j, k)] = update;
595               }
596           }
597         }
598   #pragma omp parallel for schedule(dynamic) 
599         for(i=1; i<blockDimX+1; ++i) {
600             for(int j=1; j<blockDimY+1; ++j) {
601                 for(int k=1; k<blockDimZ+1; ++k) {
602                   if(thisIndex.x == x_c && y_c == thisIndex.y && z_c == thisIndex.z && i==x_p && j==y_p && k==z_p)
603                       continue;
604                     temperature[index(i, j, k)] = new_temperature[index(i, j, k)];
605                 }
606             }
607         }
608 #if     PRINT_DEBUG
609         CkPrintf("\n {%d:%d:%d  %d:%d:%d  %d:%d:%d %f}", thisIndex_x, thisIndex_y, thisIndex_z, x_c, y_c, z_c, x_p, y_p, z_p, maxDiff);
610         for(i=1; i<blockDimX+1; ++i) {
611           for(int j=1; j<blockDimY+1; ++j) {
612               for(int k=1; k<blockDimZ+1; ++k) {
613                   CkPrintf("([%d:%d:%d]: %.2f )", i, j, k, temperature[index(i, j, k)]); 
614               }
615           }
616         }
617         CkPrintf("\n\n");
618 #endif
619         return maxDiff;
620     }
621
622 #else
623
624     // Check to see if we have received all neighbor values yet
625     // If all neighbor values have been received, we update our values and proceed
626     double Jacobi::compute_kernel() {     //Gauss-Siedal compute
627         int i=1;
628         double maxDiff=0;
629         double tmp, update=0;
630         int count=1;
631         int gg;
632
633         int x_c = arrayDimX/2/blockDimX;
634         int y_c = arrayDimY/2/blockDimY;
635         int z_c = arrayDimZ/2/blockDimZ;
636         int x_p = (arrayDimX/2)%blockDimX+1;
637         int y_p = (arrayDimY/2)%blockDimY+1;
638         int z_p = (arrayDimZ/2)%blockDimZ+1;
639             //CkPrintf("\n\n START %d [%d:%d:%d] \n", iterations, thisIndex.x, thisIndex.y, thisIndex.z);
640   #pragma omp parallel for schedule(dynamic) 
641         for(gg=0; gg<GAUSS_ITER; gg++){
642             maxDiff = 0;
643             for(i=1; i<blockDimX+1; ++i) {
644           //printf("[%d] did  %d iteration out of %d \n", omp_get_thread_num(), i, blockDimX+1); 
645           for(int j=1; j<blockDimY+1; ++j) {
646               for(int k=1; k<blockDimZ+1; ++k) {
647                   if(thisIndex.x == x_c && y_c == thisIndex.y && z_c == thisIndex.z && i==x_p && j==y_p && k==z_p)
648                       continue;
649                   update = 0;//temperature[index(i, j, k)];
650                   update += temperature[index(i+1, j, k)];
651                   update += temperature[index(i-1, j, k)];
652                   update += temperature[index(i, j+1, k)];
653                   update += temperature[index(i, j-1, k)];
654                   update += temperature[index(i, j, k+1)];
655                   update += temperature[index(i, j, k-1)];
656                   update /= 6;
657                  
658                   double diff = temperature[index(i, j, k)] - update;
659                   if(diff<0) diff = -1*diff;
660                   if(diff>maxDiff) maxDiff = diff;
661                   temperature[index(i, j, k)] = update;
662               }
663           }
664         }
665         }
666 #if     PRINT_DEBUG
667         CkPrintf("[%d:%d:%d] iteration=%d max:%.3f ==== ", thisIndex.x, thisIndex.y, thisIndex.z, iterations, maxDiff); 
668         for(i=0; i<blockDimX+2; ++i) {
669           for(int j=0; j<blockDimY+2; ++j) {
670               for(int k=0; k<blockDimZ+2; ++k) {
671                   CkPrintf("([%d:%d:%d]: %.2f )", i, j, k, temperature[index(i, j, k)]); 
672               }
673           }
674         }
675         CkPrintf("\n\n");
676 #endif
677         return maxDiff;
678     }
679 #endif
680     // Enforce some boundary conditions
681     void Jacobi::constrainBC() {
682       
683         // Heat right, left
684         if(thisIndex.x == 0 || thisIndex.x == num_chare_x-1)
685             for(int j=0; j<blockDimY+2; ++j)
686                 for(int k=0; k<blockDimZ+2; ++k)
687                 {   
688                     if(thisIndex.x == 0)
689                         temperature[index(0, j, k)] = LOW_VALUE;
690                     else
691                         temperature[index(blockDimX+1, j, k)] = LOW_VALUE;
692                 }
693         if(thisIndex.y == 0 || thisIndex.y == num_chare_y-1)
694             for(int j=0; j<blockDimX+2; ++j)
695                 for(int k=0; k<blockDimZ+2; ++k)
696                 {   
697                     if(thisIndex.y == 0)
698                         temperature[index(j,0, k)] = LOW_VALUE;
699                     else
700                         temperature[index(j, blockDimY+1, k)] = LOW_VALUE;
701                 }
702         if(thisIndex.z == 0 || thisIndex.z == num_chare_z-1)
703             for(int j=0; j<blockDimX+2; ++j)
704                 for(int k=0; k<blockDimY+2; ++k)
705                 {   
706                     if(thisIndex.z == 0)
707                         temperature[index(j, k, 0)] = LOW_VALUE;
708                     else
709                         temperature[index(j, k, blockDimZ+1)] = LOW_VALUE;
710                 }
711         int x_c = arrayDimX/2/blockDimX;
712         int y_c = arrayDimY/2/blockDimY;
713         int z_c = arrayDimZ/2/blockDimZ;
714         int x_p = (arrayDimX/2)%blockDimX+1;
715         int y_p = (arrayDimY/2)%blockDimY+1;
716         int z_p = (arrayDimZ/2)%blockDimZ+1;
717         if(thisIndex.x == x_c && y_c == thisIndex.y && z_c == thisIndex.z)
718         {
719             temperature[index(x_p, y_p, z_p)] = HIGH_VALUE;
720         }
721
722     }
723
724     void Jacobi::print()
725     {
726         CkPrintf(" print %d:%d:%d ", thisIndex.x, thisIndex.y, thisIndex.z); 
727           //printf("[%d] did  %d iteration out of %d \n", omp_get_thread_num(), i, blockDimX+1); 
728           for(int i=1; i<blockDimX+1;++i){CkPrintf("==\n");
729               for(int j=1;j<blockDimY+1;++j)
730               for(int k=1; k<blockDimZ+1; ++k) {
731                   CkPrintf(" %e ",  temperature[index(i, j, k)]) ;
732               }
733               CkPrintf("------\n-----\n");
734           }
735               contribute(0, 0, CkReduction::concat, CkCallback(CkIndex_Main::report(), mainProxy));
736     }
737 /** \class JacobiMap
738  *
739  */
740
741 class JacobiNodeMap : public CkArrayMap {
742
743 public:
744     JacobiNodeMap() {}
745     ~JacobiNodeMap() { 
746     }
747
748     int procNum(int, const CkArrayIndex &idx) {
749       int *index = (int *)idx.data();
750       return (CkMyNodeSize() * (index[0]*num_chare_x*num_chare_y + index[1]*num_chare_y + index[2]))%CkNumPes(); 
751     }
752
753 };
754 class JacobiMap : public CkArrayMap {
755   public:
756     int X, Y, Z;
757     int *mapping;
758
759     JacobiMap(int x, int y, int z) {
760       X = x; Y = y; Z = z;
761       mapping = new int[X*Y*Z];
762
763       // we are assuming that the no. of chares in each dimension is a 
764       // multiple of the torus dimension
765
766       TopoManager tmgr;
767       int dimNX, dimNY, dimNZ, dimNT;
768
769       dimNX = tmgr.getDimNX();
770       dimNY = tmgr.getDimNY();
771       dimNZ = tmgr.getDimNZ();
772       dimNT = tmgr.getDimNT();
773
774       // we are assuming that the no. of chares in each dimension is a 
775       // multiple of the torus dimension
776       int numCharesPerPe = X*Y*Z/CkNumPes();
777
778       int numCharesPerPeX = X / dimNX;
779       int numCharesPerPeY = Y / dimNY;
780       int numCharesPerPeZ = Z / dimNZ;
781       int pe = 0, pes = CkNumPes();
782
783 #if USE_BLOCK_RNDMAP
784       int used[pes];
785       for(int i=0; i<pes; i++)
786         used[i] = 0;
787 #endif
788
789       if(dimNT < 2) {   // one core per node
790         if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ); 
791         for(int i=0; i<dimNX; i++)
792           for(int j=0; j<dimNY; j++)
793             for(int k=0; k<dimNZ; k++)
794             {
795 #if USE_BLOCK_RNDMAP
796               pe = myrand(pes); 
797               while(used[pe]!=0) {
798                 pe = myrand(pes); 
799               }
800               used[pe] = 1;
801 #endif
802
803               for(int ci=i*numCharesPerPeX; ci<(i+1)*numCharesPerPeX; ci++)
804                 for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
805                   for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
806 #if USE_TOPOMAP
807                     mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k);
808 #elif USE_BLOCK_RNDMAP
809                     mapping[ci*Y*Z + cj*Z + ck] = pe;
810 #endif
811                   }
812             }
813       } else {          // multiple cores per node
814         // In this case, we split the chares in the X dimension among the
815         // cores on the same node. The strange thing I figured out is that
816         // doing this in the Z dimension is not as good.
817         numCharesPerPeX /= dimNT;
818         if(CkMyPe()==0) CkPrintf("%d %d %d %d : %d %d %d \n", dimNX, dimNY, dimNZ, dimNT, numCharesPerPeX, numCharesPerPeY, numCharesPerPeZ);
819
820         for(int i=0; i<dimNX; i++)
821           for(int j=0; j<dimNY; j++)
822             for(int k=0; k<dimNZ; k++)
823               for(int l=0; l<dimNT; l++)
824                 for(int ci=(dimNT*i+l)*numCharesPerPeX; ci<(dimNT*i+l+1)*numCharesPerPeX; ci++)
825                   for(int cj=j*numCharesPerPeY; cj<(j+1)*numCharesPerPeY; cj++)
826                     for(int ck=k*numCharesPerPeZ; ck<(k+1)*numCharesPerPeZ; ck++) {
827                       mapping[ci*Y*Z + cj*Z + ck] = tmgr.coordinatesToRank(i, j, k, l);
828                     }
829       } // end of if
830
831       if(CkMyPe() == 0) CkPrintf("Map generated ... \n");
832     }
833
834     ~JacobiMap() { 
835       delete [] mapping;
836     }
837
838     int procNum(int, const CkArrayIndex &idx) {
839       int *index = (int *)idx.data();
840       return mapping[index[0]*Y*Z + index[1]*Z + index[2]]; 
841     }
842 };
843
844 class OmpInitializer : public CBase_OmpInitializer{
845 public:
846         OmpInitializer(int numThreads) {
847 #ifdef JACOBI_OPENMP
848           if (numThreads < 1) {
849             numThreads = 1; 
850           }
851           //omp_set_num_threads(numThreads);
852           if(CkMyPe() == 0)
853           {
854 #pragma omp parallel
855               { if(omp_get_thread_num() == 0 )
856                   CkPrintf("Computation loop will be parallelized"
857                " with %d OpenMP threads\n", omp_get_num_threads());
858               }
859           }
860 #endif
861           mainProxy.start();
862         }
863 };
864
865 class TraceControl : public Group 
866 {
867 public:
868     TraceControl() {}
869
870     void startTrace() { traceBegin(); }
871     
872     void endTrace() { traceEnd(); }
873 };
874
875 #include "jacobi3d.def.h"