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