atof seems to be not thread-safe on darwin (hope.cs). make sure only rank 0 calls it.
[charm.git] / examples / bigsim / emulator / jacobi3D.C
1 //
2 // 3D Jacobi code for the Blue Gene emulator
3 // - Version 0.01 (27 Dec 2000) Chee Wai Lee
4
5 // partly rewritten by Gengbin Zheng on 2/20/2001 porting to my Converse based
6 // BlueGene emulator
7
8 #include <math.h>
9 #include <stdlib.h>
10 #include "blue.h"
11
12 #define MAX_ARRAY_SIZE 36 // the suns don't like a larger size
13 #define A_SIZE_DEFAULT 16
14 #define NUM_DBLMSG_COUNT 15
15
16 int reduceID;
17 int computeID;
18 int ghostID;
19 int exchangeID;
20 int outputID;
21
22 // source ghost region tags
23 #define LEFT            1
24 #define RIGHT           2
25 #define BELOW           3
26 #define ABOVE           4
27 #define FRONT           5
28 #define BACK            6
29
30 extern "C" void reduce(char *);
31 extern "C" void compute(char *);
32 extern "C" void ghostrecv(char *);
33 extern "C" void ghostexchange(char *);
34 extern "C" void outputData(char *);
35
36 void initArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
37                int, int, int);
38 void copyArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
39                double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
40                int, int, int,
41                int, int, int);
42 void copyXArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
43                 double [1][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
44                 int, int, int,
45                 int, int, int);
46 void copyYArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
47                 double [MAX_ARRAY_SIZE][1][MAX_ARRAY_SIZE],
48                 int, int, int,
49                 int, int, int);
50 void copyZArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
51                 double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][1],
52                 int, int, int,
53                 int, int, int);
54 void printArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
55                 int, int, int);
56
57 class reduceMsg
58 {
59 public:
60   char core[CmiBlueGeneMsgHeaderSizeBytes];
61   bool done;
62   // IMPORTANT! kernel messages are allocated using CmiAlloc
63   void *operator new(size_t s) { return CmiAlloc(s); }
64   void operator delete(void* ptr) { CmiFree(ptr); }
65 };
66
67 class computeMsg
68 {
69 public:
70   char core[CmiBlueGeneMsgHeaderSizeBytes];
71   void *operator new(size_t s) { return CmiAlloc(s); }
72   void operator delete(void* ptr) { CmiFree(ptr); }
73 };
74
75 class ghostMsg
76 {
77 public:
78   char core[CmiBlueGeneMsgHeaderSizeBytes];
79   // respect the 128-byte limit of a bluegene message
80   double data[NUM_DBLMSG_COUNT]; 
81   int datacount;
82   int source;
83   void *operator new(size_t s) { return CmiAlloc(s); }
84   void operator delete(void* ptr) { CmiFree(ptr); }
85 };
86
87 class exchangeMsg
88 {
89 public:
90   char core[CmiBlueGeneMsgHeaderSizeBytes];
91   void *operator new(size_t s) { return CmiAlloc(s); }
92   void operator delete(void* ptr) { CmiFree(ptr); }
93 };
94
95 class outputMsg
96 {
97 public:
98   char core[CmiBlueGeneMsgHeaderSizeBytes];
99   void *operator new(size_t s) { return CmiAlloc(s); }
100   void operator delete(void* ptr) { CmiFree(ptr); }
101 };
102
103 typedef struct {
104   double gdata[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
105   double maindata[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
106   double tempdata[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
107   double ghost_x1[1][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
108   double ghost_x2[1][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
109   double ghost_y1[MAX_ARRAY_SIZE][1][MAX_ARRAY_SIZE];
110   double ghost_y2[MAX_ARRAY_SIZE][1][MAX_ARRAY_SIZE];
111   double ghost_z1[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][1];
112   double ghost_z2[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][1];
113   int my_x_size;
114   int my_y_size;
115   int my_z_size;
116   int reduction_count;   // maintained only by PE[0][0][0] 
117   int ghost_x1_elements_total;
118   int ghost_x2_elements_total;
119   int ghost_y1_elements_total;
120   int ghost_y2_elements_total;
121   int ghost_z1_elements_total;
122   int ghost_z2_elements_total;
123   int ghost_x1_elements_current;
124   int ghost_x2_elements_current;
125   int ghost_y1_elements_current;
126   int ghost_y2_elements_current;
127   int ghost_z1_elements_current;
128   int ghost_z2_elements_current;
129   bool done;             // maintained only by PE[0][0][0]
130   int iteration_count;
131 } nodeData;
132
133 int g_x_blocksize;
134 int g_y_blocksize;
135 int g_z_blocksize;
136 double g_epsilon;
137 int g_iteration_count = 0;
138 int g_dataXsize = A_SIZE_DEFAULT;
139 int g_dataYsize = A_SIZE_DEFAULT;
140 int g_dataZsize = A_SIZE_DEFAULT;
141
142 void BgEmulatorInit(int argc, char **argv) 
143 {
144   if (argc < 7) {
145     CmiPrintf("Usage: jacobi3D <BG_x> <BG_y> <BG_z> <numCommTh> <numWorkTh>");
146     CmiAbort("<epsilon> [<x>] [<y>] [<z>]\n");
147   } 
148
149   BgSetSize(atoi(argv[1]), atoi(argv[2]), atoi(argv[3]));
150   BgSetNumCommThread(atoi(argv[4]));
151   BgSetNumWorkThread(atoi(argv[5]));
152
153   if (CmiMyRank() == 0)
154   g_epsilon = atof(argv[6]);
155
156     switch (argc) {
157       case 8: {
158         g_dataXsize = atoi(argv[7]);
159         break;
160       }
161       case 9: {
162         g_dataXsize = atoi(argv[7]);
163         g_dataYsize = atoi(argv[8]);
164         break;
165       }
166       case 10: {
167         g_dataXsize = atoi(argv[7]);
168         g_dataYsize = atoi(argv[8]);
169         g_dataZsize = atoi(argv[9]);
170         break;
171       }
172     }
173
174     int numBgX, numBgY, numBgZ;
175     BgGetSize(&numBgX, &numBgY, &numBgZ);
176     g_x_blocksize = g_dataXsize/numBgX;
177     g_y_blocksize = g_dataYsize/numBgY;
178     g_z_blocksize = g_dataZsize/numBgZ;
179
180     if (CmiMyPe() == 0) {
181     CmiPrintf("Bluegene size: %d %d %d\n",
182              numBgX, numBgY, numBgZ);
183     CmiPrintf("Parameters: %d %d %d, Epsilon: %e\n",
184              g_dataXsize, g_dataYsize, g_dataZsize, g_epsilon);
185     }
186 }
187
188 void BgNodeStart(int argc, char **argv) 
189 {
190   reduceID = BgRegisterHandler(reduce);
191   computeID = BgRegisterHandler(compute);
192   ghostID = BgRegisterHandler(ghostrecv);
193   exchangeID = BgRegisterHandler(ghostexchange);
194   outputID = BgRegisterHandler(outputData);
195
196   /*
197   CmiPrintf("Starting BgInit at Node %d %d %d\n", 
198            bgNode->thisIndex.x,
199            bgNode->thisIndex.y,
200            bgNode->thisIndex.z);
201   */
202   int x_start, y_start, z_start;
203   int x_end, y_end, z_end;
204
205   int x,y,z;
206   BgGetMyXYZ(&x, &y, &z);
207   // determine partition indices given BG node address by simple uniform
208   // partitioning.
209   x_start = x * g_x_blocksize;
210   y_start = y * g_y_blocksize;
211   z_start = z * g_z_blocksize;
212
213   int numBgX, numBgY, numBgZ;
214   BgGetSize(&numBgX, &numBgY, &numBgZ);
215   if (x == numBgX - 1) {
216     x_end = g_dataXsize - 1;
217   } else {
218     x_end = x_start + g_x_blocksize - 1;
219   }
220   if (y == numBgY - 1) {
221     y_end = g_dataYsize - 1;
222   } else {
223     y_end = y_start + g_y_blocksize - 1;
224   }
225   if (z == numBgZ - 1) {
226     z_end = g_dataZsize - 1;
227   } else {
228     z_end = z_start + g_z_blocksize - 1;
229   }
230
231   //create node private data
232   nodeData *nd = new nodeData;
233
234   // populate global data with default values
235   initArray(nd->gdata, g_dataXsize, g_dataYsize, g_dataZsize);
236
237   nd->my_x_size = x_end - x_start + 1;
238   nd->my_y_size = y_end - y_start + 1;
239   nd->my_z_size = z_end - z_start + 1;
240   nd->reduction_count = numBgX * numBgY * numBgZ;
241   if (x != 0) {
242     nd->ghost_x1_elements_total = (nd->my_y_size)*(nd->my_z_size);
243   } else {
244     nd->ghost_x1_elements_total = 0;
245   }
246   if (x != numBgX - 1) {
247     nd->ghost_x2_elements_total = (nd->my_y_size)*(nd->my_z_size);
248   } else {
249     nd->ghost_x2_elements_total = 0;
250   }
251   if (y != 0) {
252     nd->ghost_y1_elements_total = (nd->my_x_size)*(nd->my_z_size);
253   } else {
254     nd->ghost_y1_elements_total = 0;
255   }
256   if (y != numBgY - 1) {
257     nd->ghost_y2_elements_total = (nd->my_x_size)*(nd->my_z_size);
258   } else {
259     nd->ghost_y2_elements_total = 0;
260   }
261   if (z != 0) {
262     nd->ghost_z1_elements_total = (nd->my_x_size)*(nd->my_y_size);
263   } else {
264     nd->ghost_z1_elements_total = 0;
265   }
266   if (z != numBgZ - 1) {
267     nd->ghost_z2_elements_total = (nd->my_x_size)*(nd->my_y_size);
268   } else {
269     nd->ghost_z2_elements_total = 0;
270   }
271   nd->ghost_x1_elements_current = 0; 
272   nd->ghost_x2_elements_current = 0; 
273   nd->ghost_y1_elements_current = 0; 
274   nd->ghost_y2_elements_current = 0; 
275   nd->ghost_z1_elements_current = 0; 
276   nd->ghost_z2_elements_current = 0; 
277   nd->done = true; // assume computation complete
278   nd->iteration_count = 0;
279   /*
280   CmiPrintf("Node (%d,%d,%d) has block size %d %d %d\n", 
281            bgNode->thisIndex.x, bgNode->thisIndex.y, bgNode->thisIndex.z,
282            nd->my_x_size, nd->my_y_size, nd->my_z_size);
283   */
284   // create and populate main data from global data structure
285   copyArray(nd->gdata, nd->maindata, 
286             nd->my_x_size, nd->my_y_size, nd->my_z_size,
287             x_start, y_start, z_start);
288
289   //  printArray(nd->maindata, nd->my_x_size, nd->my_y_size, nd->my_z_size);
290
291   // create and populate the size ghost regions
292   // boundary conditions have to be checked
293   // "x-left" plane
294   if  (x != 0) {
295     copyXArray(nd->gdata, nd->ghost_x1,
296                1, nd->my_y_size, nd->my_z_size,
297                x_start - 1, y_start, z_start);
298   }
299   // "x-right" plane
300   if (x != numBgX) {
301     copyXArray(nd->gdata, nd->ghost_x2, 
302                1, nd->my_y_size, nd->my_z_size,
303                x_end + 1, y_start, z_start);
304   }
305   // "y-bottom" plane
306   if  (y != 0) {
307     copyYArray(nd->gdata, nd->ghost_y1, 
308                nd->my_x_size, 1, nd->my_z_size,
309                x_start, y_start - 1, z_start);
310   }
311   // "y-top" plane
312   if  (y != numBgY) {
313     copyYArray(nd->gdata, nd->ghost_y2, 
314                nd->my_x_size, 1, nd->my_z_size,
315                x_start, y_end + 1, z_start);
316   }
317   // "z-front" plane
318   if  (z != 0) {
319     copyZArray(nd->gdata, nd->ghost_z1, 
320                nd->my_x_size, nd->my_y_size, 1,
321                x_start, y_start, z_start - 1);
322   }
323   // "z-back" plane
324   if  (z != numBgZ) {
325     copyZArray(nd->gdata, nd->ghost_z2, 
326                nd->my_x_size, nd->my_y_size, 1,
327                x_start, y_start, z_end + 1);
328   }
329
330   computeMsg *msg = new computeMsg;
331   BgSendLocalPacket(ANYTHREAD,computeID, LARGE_WORK, sizeof(computeMsg), (char *)msg);
332
333   BgSetNodeData((char *)nd);
334 }
335
336
337 void compute(char *info) 
338 {
339   delete (computeMsg *)info;
340   bool done = true; // assume done before computation begins
341   nodeData *localdata = (nodeData *)BgGetNodeData();
342
343   int x,y,z;
344   int x_size, y_size, z_size;
345   BgGetMyXYZ(&x,&y,&z);
346   x_size = localdata->my_x_size;
347   y_size = localdata->my_y_size;
348   z_size = localdata->my_z_size;
349
350   // CmiPrintf("Node %d %d %d computing values\n",x,y,z);
351
352   int numBgX, numBgY, numBgZ;
353   BgGetSize(&numBgX, &numBgY, &numBgZ);
354
355   // one iteration of jacobi computation
356   for (int i=0; i < x_size; i++) {
357     for (int j=0; j < y_size; j++) {
358       for (int k=0; k < z_size; k++) {
359         int count = 0;
360         double sum = 0.0;
361
362         // decide if node is on x1 boundary of bluegene configuration
363         if (x != 0) {
364           // decide if element requires ghost region data
365           if (i == 0) {
366             sum += (localdata->ghost_x1)[0][j][k];
367           } else {
368             sum += (localdata->maindata)[i-1][j][k];
369           }
370           count++;
371         } else {
372           // no ghost regions to work on. just ignore.
373           if (i != 0) {
374             sum += (localdata->maindata)[i-1][j][k];
375             count++;
376           }
377         }
378         if (x != numBgX - 1) {
379           if (i == x_size - 1) {
380             sum += (localdata->ghost_x2)[0][j][k];
381           } else {
382             sum += (localdata->maindata)[i+1][j][k];
383           }
384           count++;
385         } else {
386           // no ghost regions to work on. just ignore.
387           if (i != x_size - 1) {
388             sum += (localdata->maindata)[i+1][j][k];
389             count++;
390           }
391         }
392         if (y != 0) {
393           if (j == 0) {
394             sum += (localdata->ghost_y1)[i][0][k];
395           } else {
396             sum += (localdata->maindata)[i][j-1][k];
397           }
398           count++;
399         } else {
400           // no ghost regions to work on. just ignore.
401           if (j != 0) {
402             sum += (localdata->maindata)[i][j-1][k];
403             count++;
404           }
405         }
406         if (y != numBgY - 1) {
407           if (j == y_size - 1) {
408             sum += (localdata->ghost_y2)[i][0][k];
409           } else {
410             sum += (localdata->maindata)[i][j+1][k];
411           }
412           count++;
413         } else {
414           // no ghost regions to work on. just ignore.
415           if (j != y_size - 1) {
416             sum += (localdata->maindata)[i][j+1][k];
417             count++;
418           }
419         }
420         if (z != 0) {
421           if (k == 0) {
422             sum += (localdata->ghost_z1)[i][j][0];
423           } else {
424             sum += (localdata->maindata)[i][j][k-1];
425           }
426           count++;
427         } else {
428           // no ghost regions to work on. just ignore.
429           if (k != 0) {
430             sum += (localdata->maindata)[i][j][k-1];
431             count++;
432           }
433         }
434         if (z != numBgZ - 1) {
435           if (k == z_size - 1) {
436             sum += (localdata->ghost_z2)[i][j][0];
437           } else {
438             sum += (localdata->maindata)[i][j][k+1];
439           }
440           count++;
441         } else {
442           // no ghost regions to work on. just ignore.
443           if (k != z_size - 1) {
444             sum += (localdata->maindata)[i][j][k+1];
445             count++;
446           }
447         }
448
449         (localdata->tempdata)[i][j][k] = sum / count;
450         /*
451         CmiPrintf("New: %f, Old: %f, Diff: %f, Epsilon: %f\n",
452                  (localdata->tempdata)[i][j][k],
453                  (localdata->maindata)[i][j][k],
454                  fabs((localdata->tempdata)[i][j][k] - 
455                       (localdata->maindata)[i][j][k]),
456                  g_epsilon);
457         */
458         if (fabs((localdata->maindata)[i][j][k] - 
459                  (localdata->tempdata)[i][j][k])
460             > g_epsilon) {
461           done = false;  // we're not finished yet!
462         }
463       }
464     }
465   }  // end of for loop in jacobi iteration
466   
467   copyArray(localdata->tempdata, localdata->maindata,
468             x_size, y_size, z_size,
469             0, 0, 0);
470   localdata->iteration_count++;
471   /*
472   CmiPrintf("Array computation of Node (%d,%d,%d) at iteration %d\n",
473            x, y, z,
474            localdata->iteration_count);
475   printArray(localdata->maindata,x_size,y_size,z_size);
476   */
477   // perform reduction to confirm if all regions are done. In this 
478   // simplified version, everyone sends to processor 0,0,0
479   reduceMsg *msg = new reduceMsg();
480   msg->done = done;
481
482   if ((x == 0) && (y == 0) && (z == 0)) {
483     BgSendLocalPacket(ANYTHREAD,reduceID, SMALL_WORK, sizeof(reduceMsg), (char *)msg);
484   } else {
485     BgSendPacket(0,0,0,ANYTHREAD,reduceID, SMALL_WORK, sizeof(reduceMsg), (char *)msg);
486   }
487   
488 }
489
490 void ghostexchange(char *info) {
491   delete (exchangeMsg *)info;
492   int x,y,z;
493   BgGetMyXYZ(&x,&y,&z);
494
495   /*
496   CmiPrintf("exchange procedure being activated on Node (%d,%d,%d)\n",
497            x, y, z);
498   */
499   nodeData *localdata = (nodeData *)BgGetNodeData();
500
501   double tempframe[NUM_DBLMSG_COUNT];
502
503   int numBgX, numBgY, numBgZ;
504   BgGetSize(&numBgX, &numBgY, &numBgZ);
505   // No exchange to be done! Immediately go to compute!
506   if ((numBgX == 1) &&
507       (numBgY == 1) &&
508       (numBgZ == 1)) {
509     computeMsg *msg = new computeMsg;
510     BgSendLocalPacket(ANYTHREAD,computeID, LARGE_WORK, sizeof(computeMsg), (char *)msg);
511     return;
512   }
513
514   // exchange computed ghost regions
515   // initialize message data for x1-planar ghost region
516   if (x != 0) {
517     int localcount = 0;
518     int x_size = 1;
519     int y_size = localdata->my_y_size;
520     int z_size = localdata->my_z_size;
521     for (int i=0; i < x_size; i++) {
522       for (int j=0; j < y_size; j++) {
523         for (int k=0; k < z_size; k++) {
524           tempframe[localcount] = (localdata->maindata)[i][j][k];
525           localcount++;
526           if ((localcount == NUM_DBLMSG_COUNT) ||
527               ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
528             ghostMsg *x1_msg = new ghostMsg;
529             for (int count=0; count < localcount; count++) {
530               (x1_msg->data)[count] = tempframe[count];
531             }
532             x1_msg->source = RIGHT; // sending message to left
533             x1_msg->datacount = localcount;
534             BgSendPacket(x-1,y,z,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)x1_msg);
535             localcount = 0;
536           }
537         }
538       }
539     }
540   }
541   
542   // initialize message data for x2-planar ghost region
543   if (x != numBgX - 1) {
544     int localcount = 0;
545     int x_size = 1;
546     int y_size = localdata->my_y_size;
547     int z_size = localdata->my_z_size;
548     for (int i=0; i < x_size; i++) {
549       for (int j=0; j < y_size; j++) {
550         for (int k=0; k < z_size; k++) {
551           tempframe[localcount] = 
552             (localdata->maindata)[localdata->my_x_size-1][j][k];
553           localcount++;
554           if ((localcount == NUM_DBLMSG_COUNT) ||
555               ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
556             ghostMsg *x2_msg = new ghostMsg;
557             for (int count=0;count < localcount; count++) {
558               (x2_msg->data)[count] = tempframe[count];
559             }
560             x2_msg->source = LEFT; // sending message to right
561             x2_msg->datacount = localcount;
562             BgSendPacket(x+1,y,z,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)x2_msg);
563             localcount = 0;
564           }
565         }
566       }
567     }
568   }
569
570   // initialize message data for y1-planar ghost region
571   if (y != 0) {
572     int localcount = 0;
573     int x_size = localdata->my_x_size;
574     int y_size = 1;
575     int z_size = localdata->my_z_size;
576     for (int i=0; i < x_size; i++) {
577       for (int j=0; j < y_size; j++) {
578         for (int k=0; k < z_size; k++) {
579           tempframe[localcount] = (localdata->maindata)[i][0][k];
580           localcount++;
581           if ((localcount == NUM_DBLMSG_COUNT) ||
582               ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
583             ghostMsg *y1_msg = new ghostMsg;
584             for (int count=0;count < localcount; count++) {
585               (y1_msg->data)[count] = tempframe[count];
586             }
587             y1_msg->source = ABOVE; // sending message below
588             y1_msg->datacount = localcount;
589             BgSendPacket(x,y-1,z,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)y1_msg);
590             localcount = 0;
591           }
592         }
593       }
594     }
595   }
596
597   // initialize message data for x2-planar ghost region
598   if (y != numBgY - 1) {
599     int localcount = 0;
600     int x_size = localdata->my_x_size;
601     int y_size = 1;
602     int z_size = localdata->my_z_size;
603     for (int i=0; i < x_size; i++) {
604       for (int j=0; j < y_size; j++) {
605         for (int k=0; k < z_size; k++) {
606           tempframe[localcount] = 
607             (localdata->maindata)[i][localdata->my_y_size-1][k];
608           localcount++;
609           if ((localcount == NUM_DBLMSG_COUNT) ||
610               ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
611             ghostMsg *y2_msg = new ghostMsg;
612             for (int count=0;count < localcount; count++) {
613               (y2_msg->data)[count] = tempframe[count];
614             }
615             y2_msg->source = BELOW; // sending message up
616             y2_msg->datacount = localcount;
617             BgSendPacket(x,y+1,z,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)y2_msg);
618             localcount = 0;
619           }
620         }
621       }
622     }
623   }
624   
625   // initialize message data for z1-planar ghost region
626   if (z != 0) {
627     int localcount = 0;
628     int x_size = localdata->my_x_size;
629     int y_size = localdata->my_y_size;
630     int z_size = 1;
631     for (int i=0; i < x_size; i++) {
632       for (int j=0; j < y_size; j++) {
633         for (int k=0; k < z_size; k++) {
634           tempframe[localcount] = (localdata->maindata)[i][j][0];
635           localcount++;
636           if ((localcount == NUM_DBLMSG_COUNT) ||
637               ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
638             ghostMsg *z1_msg = new ghostMsg;
639             for (int count=0; count < localcount; count++) {
640               (z1_msg->data)[count] = tempframe[count];
641             }
642             z1_msg->source = BACK; // sending message to the front
643             z1_msg->datacount = localcount;
644             BgSendPacket(x,y,z-1,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)z1_msg);
645             localcount = 0;
646           }
647         }
648       }
649     }
650   }
651
652   // initialize message data for z2-planar ghost region
653   if (z != numBgZ - 1) {
654     int localcount = 0;
655     int x_size = localdata->my_x_size;
656     int y_size = localdata->my_y_size;
657     int z_size = 1;
658     for (int i=0; i < x_size; i++) {
659       for (int j=0; j < y_size; j++) {
660         for (int k=0; k < z_size; k++) {
661           tempframe[localcount] = 
662             (localdata->maindata)[i][j][localdata->my_z_size-1];
663           localcount++;
664           if ((localcount == NUM_DBLMSG_COUNT) ||
665               ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
666             ghostMsg *z2_msg = new ghostMsg;
667             for (int count=0; count < localcount; count++) {
668               (z2_msg->data)[count] = tempframe[count];
669             }
670             z2_msg->source = FRONT; // sending message to the back
671             z2_msg->datacount = localcount;
672             BgSendPacket(x,y,z+1,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)z2_msg);
673             localcount = 0;
674           }
675         }
676       }
677     }
678   }
679 }
680
681 void reduce(char *info) 
682 {
683   // only received by processor 0,0,0
684   int x,y,z;
685   BgGetMyXYZ(&x,&y,&z);
686
687   /*
688   CmiPrintf("reduce procedure being activated on Node (%d,%d,%d)\n",
689            x, y, z);
690   */
691   nodeData *localdata = (nodeData *)BgGetNodeData();
692
693   localdata->done &= ((reduceMsg *)info)->done;
694   delete (reduceMsg *)info;
695   localdata->reduction_count--;
696
697   int numBgX, numBgY, numBgZ;
698   BgGetSize(&numBgX, &numBgY, &numBgZ);
699
700   // if that was the last message
701   if (localdata->reduction_count == 0) {
702     if (!localdata->done) {
703       // more work to be done, so ask every node to exchange ghost regions
704       for (int i=0; i<numBgX;i++) {
705         for (int j=0; j<numBgY;j++) {
706           for (int k=0; k<numBgZ;k++) {
707             exchangeMsg *msg = new exchangeMsg;
708             if ((i == 0) && (j == 0) && (k == 0)) {
709               BgSendLocalPacket(ANYTHREAD,exchangeID, LARGE_WORK, sizeof(exchangeMsg),(char *)msg);
710             } else {
711               BgSendPacket(i,j,k,ANYTHREAD,exchangeID,SMALL_WORK, sizeof(exchangeMsg),(char *)msg);
712             }
713           }
714         }
715       }
716       // resetting the value of the reduction count
717       localdata->reduction_count = numBgX * numBgY * numBgZ;
718       // resetting the completion status of the computation to assume true
719       localdata->done = true; 
720     } else {
721       // instruct all computing nodes to output their individual results
722       for (int i=0; i<numBgX;i++) {
723         for (int j=0; j<numBgY;j++) {
724           for (int k=0; k<numBgZ;k++) {
725             outputMsg *msg = new outputMsg;
726             if ((i == 0) && (j == 0) && (k == 0)) {
727               BgSendLocalPacket(ANYTHREAD,outputID, LARGE_WORK, sizeof(outputMsg), (char *)msg);
728             } else {
729               BgSendPacket(i,j,k,ANYTHREAD,outputID,SMALL_WORK, sizeof(outputMsg), (char *)msg);
730             }
731           }
732         }
733       }
734     }
735   } 
736 }
737
738 void ghostrecv(char *info)
739 {
740   int x,y,z;
741   BgGetMyXYZ(&x,&y,&z);
742
743   /*
744   CmiPrintf("Node (%d,%d,%d) receiving exchange message\n",
745            x, y, z);
746   */
747   nodeData *localdata = (nodeData *)BgGetNodeData();
748   
749   int x_size = localdata->my_x_size;
750   int y_size = localdata->my_y_size;
751   int z_size = localdata->my_z_size;
752
753   int x1_total = localdata->ghost_x1_elements_total;
754   int x2_total = localdata->ghost_x2_elements_total;
755   int y1_total = localdata->ghost_y1_elements_total;
756   int y2_total = localdata->ghost_y2_elements_total;
757   int z1_total = localdata->ghost_z1_elements_total;
758   int z2_total = localdata->ghost_z2_elements_total;
759
760   // determine the source of the ghost region and copy the data to the
761   // appropriate local ghost region.
762   if (((ghostMsg *)info)->source == LEFT) {
763     for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
764       int x_offset = 
765         ((localdata->ghost_x1_elements_current)/(y_size*z_size)) % 1;
766       int y_offset =
767         ((localdata->ghost_x1_elements_current)/z_size) % y_size;
768       int z_offset =
769         (localdata->ghost_x1_elements_current) % z_size;
770       (localdata->ghost_x1)[x_offset][y_offset][z_offset] =
771         ((ghostMsg*)info)->data[count];
772       localdata->ghost_x1_elements_current++;
773     }
774   } else if (((ghostMsg *)info)->source == RIGHT) {
775     for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
776       int x_offset = 
777         ((localdata->ghost_x2_elements_current)/(y_size*z_size)) % 1;
778       int y_offset =
779         ((localdata->ghost_x2_elements_current)/z_size) % y_size;
780       int z_offset =
781         (localdata->ghost_x2_elements_current) % z_size;
782       (localdata->ghost_x2)[x_offset][y_offset][z_offset] =
783         ((ghostMsg*)info)->data[count];
784       localdata->ghost_x2_elements_current++;
785     }
786   } else if (((ghostMsg *)info)->source == BELOW) {
787     for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
788       int x_offset = 
789         ((localdata->ghost_y1_elements_current)/(1*z_size)) % x_size;
790       int y_offset =
791         ((localdata->ghost_y1_elements_current)/z_size) % 1;
792       int z_offset =
793         (localdata->ghost_y1_elements_current) % z_size;
794       (localdata->ghost_y1)[x_offset][y_offset][z_offset] =
795         ((ghostMsg*)info)->data[count];
796       localdata->ghost_y1_elements_current++;
797     }
798   } else if (((ghostMsg *)info)->source == ABOVE) {
799     for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
800       int x_offset = 
801         ((localdata->ghost_y2_elements_current)/(1*z_size)) % x_size;
802       int y_offset =
803         ((localdata->ghost_y2_elements_current)/z_size) % 1;
804       int z_offset =
805         (localdata->ghost_y2_elements_current) % z_size;
806       (localdata->ghost_y2)[x_offset][y_offset][z_offset] =
807         ((ghostMsg*)info)->data[count];
808       localdata->ghost_y2_elements_current++;
809     }
810   } else if (((ghostMsg *)info)->source == FRONT) {
811     for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
812       int x_offset = 
813         ((localdata->ghost_z1_elements_current)/(y_size*1)) % x_size;
814       int y_offset =
815         ((localdata->ghost_z1_elements_current)/1) % y_size;
816       int z_offset =
817         (localdata->ghost_z1_elements_current) % 1;
818       (localdata->ghost_z1)[x_offset][y_offset][z_offset] =
819         ((ghostMsg*)info)->data[count];
820       localdata->ghost_z1_elements_current++;
821     }
822   } else if (((ghostMsg *)info)->source == BACK) {
823     for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
824       int x_offset = 
825         ((localdata->ghost_z2_elements_current)/(y_size*1)) % x_size;
826       int y_offset =
827         ((localdata->ghost_z2_elements_current)/1) % y_size;
828       int z_offset =
829         (localdata->ghost_z2_elements_current) % 1;
830       (localdata->ghost_z2)[x_offset][y_offset][z_offset] =
831         ((ghostMsg*)info)->data[count];
832       localdata->ghost_z2_elements_current++;
833     }
834   }
835
836   // if that was the last message
837   if ((x1_total == localdata->ghost_x1_elements_current) &&
838       (x2_total == localdata->ghost_x2_elements_current) &&
839       (y1_total == localdata->ghost_y1_elements_current) &&
840       (y2_total == localdata->ghost_y2_elements_current) &&
841       (z1_total == localdata->ghost_z1_elements_current) &&
842       (z2_total == localdata->ghost_z2_elements_current)) {
843     // reset exchange counts
844     localdata->ghost_x1_elements_current = 0;
845     localdata->ghost_x2_elements_current = 0;
846     localdata->ghost_y1_elements_current = 0;
847     localdata->ghost_y2_elements_current = 0;
848     localdata->ghost_z1_elements_current = 0;
849     localdata->ghost_z2_elements_current = 0;
850     
851     computeMsg *msg = new computeMsg;
852     BgSendLocalPacket(ANYTHREAD,computeID, LARGE_WORK, sizeof(computeMsg), (char *)msg); // get to work!
853   }
854   delete (ghostMsg*)info;
855 }
856
857 void outputData(char *info) {
858   delete (outputMsg *)info;
859   int x,y,z;
860   BgGetMyXYZ(&x,&y,&z);
861
862   /*
863   CmiPrintf("Node (%d,%d,%d) printing data.\n",
864            x, y, z);
865   */
866   nodeData *localdata = (nodeData *)BgGetNodeData();
867   
868   int x_size = localdata->my_x_size;
869   int y_size = localdata->my_y_size;
870   int z_size = localdata->my_z_size;
871
872   CmiPrintf("Final output at Node (%d,%d,%d) with iteration count = %d:\n",
873            x, y, z,
874            localdata->iteration_count);
875 //  printArray(localdata->maindata,x_size,y_size,z_size);
876
877   if ((x == 0) && (y == 0) && (z == 0)) {
878     BgShutdown();
879   }
880 }
881
882 void initArray(double target[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
883                int x, int y, int z) {
884   for (int i=0; i < x; i++) {
885     for (int j=0; j < y; j++) {
886       for (int k=0; k < z; k++) {
887         target[i][j][k] = (i+j+k)*1.0;
888       }
889     }
890   }
891 }
892
893 void copyArray(double source[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
894                double dest[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
895                int x, int y, int z,
896                int offset_x, int offset_y, int offset_z) {
897   for (int i=0; i < x; i++) {
898     for (int j=0; j < y; j++) {
899       for (int k=0; k < z; k++) {
900         dest[i][j][k] = source[i+offset_x][j+offset_y][k+offset_z];
901       }
902     }
903   }
904 }
905
906 void copyXArray(double source[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
907                 double dest[1][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
908                 int x, int y, int z,
909                 int offset_x, int offset_y, int offset_z) {
910   for (int i=0; i < x; i++) {
911     for (int j=0; j < y; j++) {
912       for (int k=0; k < z; k++) {
913         dest[i][j][k] = source[i+offset_x][j+offset_y][k+offset_z];
914       }
915     }
916   }
917 }
918
919 void copyYArray(double source[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
920                 double dest[MAX_ARRAY_SIZE][1][MAX_ARRAY_SIZE],
921                 int x, int y, int z,
922                 int offset_x, int offset_y, int offset_z) {
923   for (int i=0; i < x; i++) {
924     for (int j=0; j < y; j++) {
925       for (int k=0; k < z; k++) {
926         dest[i][j][k] = source[i+offset_x][j+offset_y][k+offset_z];
927       }
928     }
929   }
930 }
931
932 void copyZArray(double source[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
933                 double dest[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][1],
934                 int x, int y, int z,
935                 int offset_x, int offset_y, int offset_z) {
936   for (int i=0; i < x; i++) {
937     for (int j=0; j < y; j++) {
938       for (int k=0; k < z; k++) {
939         dest[i][j][k] = source[i+offset_x][j+offset_y][k+offset_z];
940       }
941     }
942   }
943 }
944
945 void printArray(double target[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
946                 int x, int y, int z) {
947   for (int i=0; i < x; i++) {
948     for (int j=0; j < y; j++) {
949       for (int k=0; k < z; k++) {
950         CmiPrintf("%f ",target[i][j][k]);
951       }
952       CmiPrintf("\n");
953     }
954   }
955 }
956