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