emulator examples
authorGengbin Zheng <gzheng@illinois.edu>
Thu, 7 Oct 2004 06:40:49 +0000 (06:40 +0000)
committerGengbin Zheng <gzheng@illinois.edu>
Thu, 7 Oct 2004 06:40:49 +0000 (06:40 +0000)
15 files changed:
examples/bigsim/emulator/Makefile [new file with mode: 0644]
examples/bigsim/emulator/jacobi3D.C [new file with mode: 0644]
examples/bigsim/emulator/line.C [new file with mode: 0644]
examples/bigsim/emulator/littleMD/Atom.h [new file with mode: 0644]
examples/bigsim/emulator/littleMD/Cell.h [new file with mode: 0644]
examples/bigsim/emulator/littleMD/Handlers.C [new file with mode: 0644]
examples/bigsim/emulator/littleMD/Makefile [new file with mode: 0644]
examples/bigsim/emulator/littleMD/SimParams.h [new file with mode: 0644]
examples/bigsim/emulator/littleMD/bgMD.C [new file with mode: 0644]
examples/bigsim/emulator/littleMD/bgMD.h [new file with mode: 0644]
examples/bigsim/emulator/maxReduce.C [new file with mode: 0644]
examples/bigsim/emulator/octo.C [new file with mode: 0644]
examples/bigsim/emulator/ping.C [new file with mode: 0644]
examples/bigsim/emulator/prime.C [new file with mode: 0644]
examples/bigsim/emulator/ring.C [new file with mode: 0644]

diff --git a/examples/bigsim/emulator/Makefile b/examples/bigsim/emulator/Makefile
new file mode 100644 (file)
index 0000000..3083f71
--- /dev/null
@@ -0,0 +1,41 @@
+CHARMC=../../../bin/charmc $(OPTS)
+
+all: jacobi3D maxReduce prime ping octo line ring
+       cd littleMD; make
+
+jacobi3D: jacobi3D.C
+       $(CHARMC) -language bluegene -o jacobi3D jacobi3D.C #-memory paranoid
+
+maxReduce: maxReduce.C
+       $(CHARMC) -language bluegene -o maxReduce maxReduce.C
+
+prime: prime.C
+       $(CHARMC) -language bluegene -o prime prime.C
+
+ping: ping.C
+       $(CHARMC) -language bluegene -o ping ping.C
+
+octo: octo.C
+       $(CHARMC) -language bluegene -o octo octo.C
+
+line: line.C
+       $(CHARMC) -language bluegene -o line line.C
+
+ring: ring.C
+       $(CHARMC) -language bluegene -o ring ring.C
+
+test: all
+       ./charmrun ./maxReduce +p4 +cth3 +wth10
+       ./charmrun +p4 ./octo 4 5 6 4 5
+       ./charmrun +p4 ./line 6 6 6 5 6 10
+       ./charmrun +p4 ./jacobi3D 3 3 3 2 10 0.1
+       ./charmrun +p4 ./prime 3 3 3 3 4 1000
+       ./charmrun +p4 ./ring 2 2 2 1 1
+       ./charmrun +p4 littleMD/bgMD 6 6 6 5 6
+
+clean:
+       rm -f core *.cpm.h
+       rm -f TAGS *.o
+       rm -f jacobi3D maxReduce prime ping octo line
+       rm -f conv-host charmrun
+       cd littleMD; make clean
diff --git a/examples/bigsim/emulator/jacobi3D.C b/examples/bigsim/emulator/jacobi3D.C
new file mode 100644 (file)
index 0000000..9c82928
--- /dev/null
@@ -0,0 +1,953 @@
+//
+// 3D Jacobi code for the Blue Gene emulator
+// - Version 0.01 (27 Dec 2000) Chee Wai Lee
+
+// partly rewritten by Gengbin Zheng on 2/20/2001 porting to my Converse based
+// BlueGene emulator
+
+#include <math.h>
+#include <stdlib.h>
+#include "blue.h"
+
+#define MAX_ARRAY_SIZE 36 // the suns don't like a larger size
+#define A_SIZE_DEFAULT 16
+#define NUM_DBLMSG_COUNT 15
+
+int reduceID;
+int computeID;
+int ghostID;
+int exchangeID;
+int outputID;
+
+// source ghost region tags
+#define LEFT            1
+#define RIGHT           2
+#define BELOW           3
+#define ABOVE           4
+#define FRONT           5
+#define BACK            6
+
+extern "C" void reduce(char *);
+extern "C" void compute(char *);
+extern "C" void ghostrecv(char *);
+extern "C" void ghostexchange(char *);
+extern "C" void outputData(char *);
+
+void initArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+              int, int, int);
+void copyArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+              double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+              int, int, int,
+              int, int, int);
+void copyXArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+               double [1][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+               int, int, int,
+               int, int, int);
+void copyYArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+               double [MAX_ARRAY_SIZE][1][MAX_ARRAY_SIZE],
+               int, int, int,
+               int, int, int);
+void copyZArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+               double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][1],
+               int, int, int,
+               int, int, int);
+void printArray(double [MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+               int, int, int);
+
+class reduceMsg
+{
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  bool done;
+  // IMPORTANT! kernel messages are allocated using CmiAlloc
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+class computeMsg
+{
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+class ghostMsg
+{
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  // respect the 128-byte limit of a bluegene message
+  double data[NUM_DBLMSG_COUNT]; 
+  int datacount;
+  int source;
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+class exchangeMsg
+{
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+class outputMsg
+{
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+typedef struct {
+  double maindata[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
+  double tempdata[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
+  double ghost_x1[1][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
+  double ghost_x2[1][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
+  double ghost_y1[MAX_ARRAY_SIZE][1][MAX_ARRAY_SIZE];
+  double ghost_y2[MAX_ARRAY_SIZE][1][MAX_ARRAY_SIZE];
+  double ghost_z1[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][1];
+  double ghost_z2[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][1];
+  int my_x_size;
+  int my_y_size;
+  int my_z_size;
+  int reduction_count;   // maintained only by PE[0][0][0] 
+  int ghost_x1_elements_total;
+  int ghost_x2_elements_total;
+  int ghost_y1_elements_total;
+  int ghost_y2_elements_total;
+  int ghost_z1_elements_total;
+  int ghost_z2_elements_total;
+  int ghost_x1_elements_current;
+  int ghost_x2_elements_current;
+  int ghost_y1_elements_current;
+  int ghost_y2_elements_current;
+  int ghost_z1_elements_current;
+  int ghost_z2_elements_current;
+  bool done;             // maintained only by PE[0][0][0]
+  int iteration_count;
+} nodeData;
+
+double gdata[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
+int g_x_blocksize;
+int g_y_blocksize;
+int g_z_blocksize;
+double g_epsilon;
+int g_iteration_count = 0;
+int g_dataXsize = A_SIZE_DEFAULT;
+int g_dataYsize = A_SIZE_DEFAULT;
+int g_dataZsize = A_SIZE_DEFAULT;
+
+void BgEmulatorInit(int argc, char **argv) 
+{
+  if (argc < 7) {
+    CmiPrintf("Usage: jacobi3D <BG_x> <BG_y> <BG_z> <numCommTh> <numWorkTh>");
+    CmiAbort("<epsilon> [<x>] [<y>] [<z>]\n");
+  } 
+
+  BgSetSize(atoi(argv[1]), atoi(argv[2]), atoi(argv[3]));
+  BgSetNumCommThread(atoi(argv[4]));
+  BgSetNumWorkThread(atoi(argv[5]));
+
+  g_epsilon = atof(argv[6]);
+
+    switch (argc) {
+      case 8: {
+       g_dataXsize = atoi(argv[7]);
+       break;
+      }
+      case 9: {
+       g_dataXsize = atoi(argv[7]);
+       g_dataYsize = atoi(argv[8]);
+       break;
+      }
+      case 10: {
+       g_dataXsize = atoi(argv[7]);
+       g_dataYsize = atoi(argv[8]);
+       g_dataZsize = atoi(argv[9]);
+       break;
+      }
+    }
+
+    int numBgX, numBgY, numBgZ;
+    BgGetSize(&numBgX, &numBgY, &numBgZ);
+    g_x_blocksize = g_dataXsize/numBgX;
+    g_y_blocksize = g_dataYsize/numBgY;
+    g_z_blocksize = g_dataZsize/numBgZ;
+
+    // populate global data with default values
+    initArray(gdata, g_dataXsize, g_dataYsize, g_dataZsize);
+
+    CmiPrintf("Bluegene size: %d %d %d\n",
+            numBgX, numBgY, numBgZ);
+    CmiPrintf("Parameters: %d %d %d, Epsilon: %f\n",
+            g_dataXsize, g_dataYsize, g_dataZsize, g_epsilon);
+
+}
+
+void BgNodeStart(int argc, char **argv) 
+{
+  reduceID = BgRegisterHandler(reduce);
+  computeID = BgRegisterHandler(compute);
+  ghostID = BgRegisterHandler(ghostrecv);
+  exchangeID = BgRegisterHandler(ghostexchange);
+  outputID = BgRegisterHandler(outputData);
+
+  /*
+  CmiPrintf("Starting BgInit at Node %d %d %d\n", 
+          bgNode->thisIndex.x,
+          bgNode->thisIndex.y,
+          bgNode->thisIndex.z);
+  */
+  int x_start, y_start, z_start;
+  int x_end, y_end, z_end;
+
+  int x,y,z;
+  BgGetMyXYZ(&x, &y, &z);
+  // determine partition indices given BG node address by simple uniform
+  // partitioning.
+  x_start = x * g_x_blocksize;
+  y_start = y * g_y_blocksize;
+  z_start = z * g_z_blocksize;
+
+  int numBgX, numBgY, numBgZ;
+  BgGetSize(&numBgX, &numBgY, &numBgZ);
+  if (x == numBgX - 1) {
+    x_end = g_dataXsize - 1;
+  } else {
+    x_end = x_start + g_x_blocksize - 1;
+  }
+  if (y == numBgY - 1) {
+    y_end = g_dataYsize - 1;
+  } else {
+    y_end = y_start + g_y_blocksize - 1;
+  }
+  if (z == numBgZ - 1) {
+    z_end = g_dataZsize - 1;
+  } else {
+    z_end = z_start + g_z_blocksize - 1;
+  }
+
+  //create node private data
+  nodeData *nd = new nodeData;
+  nd->my_x_size = x_end - x_start + 1;
+  nd->my_y_size = y_end - y_start + 1;
+  nd->my_z_size = z_end - z_start + 1;
+  nd->reduction_count = numBgX * numBgY * numBgZ;
+  if (x != 0) {
+    nd->ghost_x1_elements_total = (nd->my_y_size)*(nd->my_z_size);
+  } else {
+    nd->ghost_x1_elements_total = 0;
+  }
+  if (x != numBgX - 1) {
+    nd->ghost_x2_elements_total = (nd->my_y_size)*(nd->my_z_size);
+  } else {
+    nd->ghost_x2_elements_total = 0;
+  }
+  if (y != 0) {
+    nd->ghost_y1_elements_total = (nd->my_x_size)*(nd->my_z_size);
+  } else {
+    nd->ghost_y1_elements_total = 0;
+  }
+  if (y != numBgY - 1) {
+    nd->ghost_y2_elements_total = (nd->my_x_size)*(nd->my_z_size);
+  } else {
+    nd->ghost_y2_elements_total = 0;
+  }
+  if (z != 0) {
+    nd->ghost_z1_elements_total = (nd->my_x_size)*(nd->my_y_size);
+  } else {
+    nd->ghost_z1_elements_total = 0;
+  }
+  if (z != numBgZ - 1) {
+    nd->ghost_z2_elements_total = (nd->my_x_size)*(nd->my_y_size);
+  } else {
+    nd->ghost_z2_elements_total = 0;
+  }
+  nd->ghost_x1_elements_current = 0; 
+  nd->ghost_x2_elements_current = 0; 
+  nd->ghost_y1_elements_current = 0; 
+  nd->ghost_y2_elements_current = 0; 
+  nd->ghost_z1_elements_current = 0; 
+  nd->ghost_z2_elements_current = 0; 
+  nd->done = true; // assume computation complete
+  nd->iteration_count = 0;
+  /*
+  CmiPrintf("Node (%d,%d,%d) has block size %d %d %d\n", 
+          bgNode->thisIndex.x, bgNode->thisIndex.y, bgNode->thisIndex.z,
+          nd->my_x_size, nd->my_y_size, nd->my_z_size);
+  */
+  // create and populate main data from global data structure
+  copyArray(gdata, nd->maindata, 
+           nd->my_x_size, nd->my_y_size, nd->my_z_size,
+           x_start, y_start, z_start);
+
+  //  printArray(nd->maindata, nd->my_x_size, nd->my_y_size, nd->my_z_size);
+
+  // create and populate the size ghost regions
+  // boundary conditions have to be checked
+  // "x-left" plane
+  if  (x != 0) {
+    copyXArray(gdata, nd->ghost_x1,
+              1, nd->my_y_size, nd->my_z_size,
+              x_start - 1, y_start, z_start);
+  }
+  // "x-right" plane
+  if (x != numBgX) {
+    copyXArray(gdata, nd->ghost_x2, 
+              1, nd->my_y_size, nd->my_z_size,
+              x_end + 1, y_start, z_start);
+  }
+  // "y-bottom" plane
+  if  (y != 0) {
+    copyYArray(gdata, nd->ghost_y1, 
+              nd->my_x_size, 1, nd->my_z_size,
+              x_start, y_start - 1, z_start);
+  }
+  // "y-top" plane
+  if  (y != numBgY) {
+    copyYArray(gdata, nd->ghost_y2, 
+              nd->my_x_size, 1, nd->my_z_size,
+              x_start, y_end + 1, z_start);
+  }
+  // "z-front" plane
+  if  (z != 0) {
+    copyZArray(gdata, nd->ghost_z1, 
+              nd->my_x_size, nd->my_y_size, 1,
+              x_start, y_start, z_start - 1);
+  }
+  // "z-back" plane
+  if  (z != numBgZ) {
+    copyZArray(gdata, nd->ghost_z2, 
+              nd->my_x_size, nd->my_y_size, 1,
+              x_start, y_start, z_end + 1);
+  }
+
+  computeMsg *msg = new computeMsg;
+  BgSendLocalPacket(ANYTHREAD,computeID, LARGE_WORK, sizeof(computeMsg), (char *)msg);
+
+  BgSetNodeData((char *)nd);
+}
+
+
+void compute(char *info) 
+{
+  delete (computeMsg *)info;
+  bool done = true; // assume done before computation begins
+  nodeData *localdata = (nodeData *)BgGetNodeData();
+
+  int x,y,z;
+  int x_size, y_size, z_size;
+  BgGetMyXYZ(&x,&y,&z);
+  x_size = localdata->my_x_size;
+  y_size = localdata->my_y_size;
+  z_size = localdata->my_z_size;
+
+  // CmiPrintf("Node %d %d %d computing values\n",x,y,z);
+
+  int numBgX, numBgY, numBgZ;
+  BgGetSize(&numBgX, &numBgY, &numBgZ);
+
+  // one iteration of jacobi computation
+  for (int i=0; i < x_size; i++) {
+    for (int j=0; j < y_size; j++) {
+      for (int k=0; k < z_size; k++) {
+       int count = 0;
+       double sum = 0.0;
+
+       // decide if node is on x1 boundary of bluegene configuration
+       if (x != 0) {
+         // decide if element requires ghost region data
+         if (i == 0) {
+           sum += (localdata->ghost_x1)[0][j][k];
+         } else {
+           sum += (localdata->maindata)[i-1][j][k];
+         }
+         count++;
+       } else {
+         // no ghost regions to work on. just ignore.
+         if (i != 0) {
+           sum += (localdata->maindata)[i-1][j][k];
+           count++;
+         }
+       }
+       if (x != numBgX - 1) {
+         if (i == x_size - 1) {
+           sum += (localdata->ghost_x2)[0][j][k];
+         } else {
+           sum += (localdata->maindata)[i+1][j][k];
+         }
+         count++;
+       } else {
+         // no ghost regions to work on. just ignore.
+         if (i != x_size - 1) {
+           sum += (localdata->maindata)[i+1][j][k];
+           count++;
+         }
+       }
+       if (y != 0) {
+         if (j == 0) {
+           sum += (localdata->ghost_y1)[i][0][k];
+         } else {
+           sum += (localdata->maindata)[i][j-1][k];
+         }
+         count++;
+       } else {
+         // no ghost regions to work on. just ignore.
+         if (j != 0) {
+           sum += (localdata->maindata)[i][j-1][k];
+           count++;
+         }
+       }
+       if (y != numBgY - 1) {
+         if (j == y_size - 1) {
+           sum += (localdata->ghost_y2)[i][0][k];
+         } else {
+           sum += (localdata->maindata)[i][j+1][k];
+         }
+         count++;
+       } else {
+         // no ghost regions to work on. just ignore.
+         if (j != y_size - 1) {
+           sum += (localdata->maindata)[i][j+1][k];
+           count++;
+         }
+       }
+       if (z != 0) {
+         if (k == 0) {
+           sum += (localdata->ghost_z1)[i][j][0];
+         } else {
+           sum += (localdata->maindata)[i][j][k-1];
+         }
+         count++;
+       } else {
+         // no ghost regions to work on. just ignore.
+         if (k != 0) {
+           sum += (localdata->maindata)[i][j][k-1];
+           count++;
+         }
+       }
+       if (z != numBgZ - 1) {
+         if (k == z_size - 1) {
+           sum += (localdata->ghost_z2)[i][j][0];
+         } else {
+           sum += (localdata->maindata)[i][j][k+1];
+         }
+         count++;
+       } else {
+         // no ghost regions to work on. just ignore.
+         if (k != z_size - 1) {
+           sum += (localdata->maindata)[i][j][k+1];
+           count++;
+         }
+       }
+
+       (localdata->tempdata)[i][j][k] = sum / count;
+       /*
+       CmiPrintf("New: %f, Old: %f, Diff: %f, Epsilon: %f\n",
+                (localdata->tempdata)[i][j][k],
+                (localdata->maindata)[i][j][k],
+                fabs((localdata->tempdata)[i][j][k] - 
+                     (localdata->maindata)[i][j][k]),
+                g_epsilon);
+       */
+       if (fabs((localdata->maindata)[i][j][k] - 
+                (localdata->tempdata)[i][j][k])
+           > g_epsilon) {
+         done = false;  // we're not finished yet!
+       }
+      }
+    }
+  }  // end of for loop in jacobi iteration
+  
+  copyArray(localdata->tempdata, localdata->maindata,
+           x_size, y_size, z_size,
+           0, 0, 0);
+  localdata->iteration_count++;
+  /*
+  CmiPrintf("Array computation of Node (%d,%d,%d) at iteration %d\n",
+          x, y, z,
+          localdata->iteration_count);
+  printArray(localdata->maindata,x_size,y_size,z_size);
+  */
+  // perform reduction to confirm if all regions are done. In this 
+  // simplified version, everyone sends to processor 0,0,0
+  reduceMsg *msg = new reduceMsg();
+  msg->done = done;
+
+  if ((x == 0) && (y == 0) && (z == 0)) {
+    BgSendLocalPacket(ANYTHREAD,reduceID, SMALL_WORK, sizeof(reduceMsg), (char *)msg);
+  } else {
+    BgSendPacket(0,0,0,ANYTHREAD,reduceID, SMALL_WORK, sizeof(reduceMsg), (char *)msg);
+  }
+  
+}
+
+void ghostexchange(char *info) {
+  delete (exchangeMsg *)info;
+  int x,y,z;
+  BgGetMyXYZ(&x,&y,&z);
+
+  /*
+  CmiPrintf("exchange procedure being activated on Node (%d,%d,%d)\n",
+          x, y, z);
+  */
+  nodeData *localdata = (nodeData *)BgGetNodeData();
+
+  double tempframe[NUM_DBLMSG_COUNT];
+
+  int numBgX, numBgY, numBgZ;
+  BgGetSize(&numBgX, &numBgY, &numBgZ);
+  // No exchange to be done! Immediately go to compute!
+  if ((numBgX == 1) &&
+      (numBgY == 1) &&
+      (numBgZ == 1)) {
+    computeMsg *msg = new computeMsg;
+    BgSendLocalPacket(ANYTHREAD,computeID, LARGE_WORK, sizeof(computeMsg), (char *)msg);
+    return;
+  }
+
+  // exchange computed ghost regions
+  // initialize message data for x1-planar ghost region
+  if (x != 0) {
+    int localcount = 0;
+    int x_size = 1;
+    int y_size = localdata->my_y_size;
+    int z_size = localdata->my_z_size;
+    for (int i=0; i < x_size; i++) {
+      for (int j=0; j < y_size; j++) {
+       for (int k=0; k < z_size; k++) {
+         tempframe[localcount] = (localdata->maindata)[i][j][k];
+         localcount++;
+         if ((localcount == NUM_DBLMSG_COUNT) ||
+             ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
+           ghostMsg *x1_msg = new ghostMsg;
+           for (int count=0; count < localcount; count++) {
+             (x1_msg->data)[count] = tempframe[count];
+           }
+           x1_msg->source = RIGHT; // sending message to left
+           x1_msg->datacount = localcount;
+           BgSendPacket(x-1,y,z,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)x1_msg);
+           localcount = 0;
+         }
+       }
+      }
+    }
+  }
+  
+  // initialize message data for x2-planar ghost region
+  if (x != numBgX - 1) {
+    int localcount = 0;
+    int x_size = 1;
+    int y_size = localdata->my_y_size;
+    int z_size = localdata->my_z_size;
+    for (int i=0; i < x_size; i++) {
+      for (int j=0; j < y_size; j++) {
+       for (int k=0; k < z_size; k++) {
+         tempframe[localcount] = 
+           (localdata->maindata)[localdata->my_x_size-1][j][k];
+         localcount++;
+         if ((localcount == NUM_DBLMSG_COUNT) ||
+             ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
+           ghostMsg *x2_msg = new ghostMsg;
+           for (int count=0;count < localcount; count++) {
+             (x2_msg->data)[count] = tempframe[count];
+           }
+           x2_msg->source = LEFT; // sending message to right
+           x2_msg->datacount = localcount;
+           BgSendPacket(x+1,y,z,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)x2_msg);
+           localcount = 0;
+         }
+       }
+      }
+    }
+  }
+
+  // initialize message data for y1-planar ghost region
+  if (y != 0) {
+    int localcount = 0;
+    int x_size = localdata->my_x_size;
+    int y_size = 1;
+    int z_size = localdata->my_z_size;
+    for (int i=0; i < x_size; i++) {
+      for (int j=0; j < y_size; j++) {
+       for (int k=0; k < z_size; k++) {
+         tempframe[localcount] = (localdata->maindata)[i][0][k];
+         localcount++;
+         if ((localcount == NUM_DBLMSG_COUNT) ||
+             ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
+           ghostMsg *y1_msg = new ghostMsg;
+           for (int count=0;count < localcount; count++) {
+             (y1_msg->data)[count] = tempframe[count];
+           }
+           y1_msg->source = ABOVE; // sending message below
+           y1_msg->datacount = localcount;
+           BgSendPacket(x,y-1,z,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)y1_msg);
+           localcount = 0;
+         }
+       }
+      }
+    }
+  }
+
+  // initialize message data for x2-planar ghost region
+  if (y != numBgY - 1) {
+    int localcount = 0;
+    int x_size = localdata->my_x_size;
+    int y_size = 1;
+    int z_size = localdata->my_z_size;
+    for (int i=0; i < x_size; i++) {
+      for (int j=0; j < y_size; j++) {
+       for (int k=0; k < z_size; k++) {
+         tempframe[localcount] = 
+           (localdata->maindata)[i][localdata->my_y_size-1][k];
+         localcount++;
+         if ((localcount == NUM_DBLMSG_COUNT) ||
+             ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
+           ghostMsg *y2_msg = new ghostMsg;
+           for (int count=0;count < localcount; count++) {
+             (y2_msg->data)[count] = tempframe[count];
+           }
+           y2_msg->source = BELOW; // sending message up
+           y2_msg->datacount = localcount;
+           BgSendPacket(x,y+1,z,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)y2_msg);
+           localcount = 0;
+         }
+       }
+      }
+    }
+  }
+  
+  // initialize message data for z1-planar ghost region
+  if (z != 0) {
+    int localcount = 0;
+    int x_size = localdata->my_x_size;
+    int y_size = localdata->my_y_size;
+    int z_size = 1;
+    for (int i=0; i < x_size; i++) {
+      for (int j=0; j < y_size; j++) {
+       for (int k=0; k < z_size; k++) {
+         tempframe[localcount] = (localdata->maindata)[i][j][0];
+         localcount++;
+         if ((localcount == NUM_DBLMSG_COUNT) ||
+             ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
+           ghostMsg *z1_msg = new ghostMsg;
+           for (int count=0; count < localcount; count++) {
+             (z1_msg->data)[count] = tempframe[count];
+           }
+           z1_msg->source = BACK; // sending message to the front
+           z1_msg->datacount = localcount;
+           BgSendPacket(x,y,z-1,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)z1_msg);
+           localcount = 0;
+         }
+       }
+      }
+    }
+  }
+
+  // initialize message data for z2-planar ghost region
+  if (z != numBgZ - 1) {
+    int localcount = 0;
+    int x_size = localdata->my_x_size;
+    int y_size = localdata->my_y_size;
+    int z_size = 1;
+    for (int i=0; i < x_size; i++) {
+      for (int j=0; j < y_size; j++) {
+       for (int k=0; k < z_size; k++) {
+         tempframe[localcount] = 
+           (localdata->maindata)[i][j][localdata->my_z_size-1];
+         localcount++;
+         if ((localcount == NUM_DBLMSG_COUNT) ||
+             ((i == x_size-1) && (j == y_size-1) && (k == z_size-1))) {
+           ghostMsg *z2_msg = new ghostMsg;
+           for (int count=0; count < localcount; count++) {
+             (z2_msg->data)[count] = tempframe[count];
+           }
+           z2_msg->source = FRONT; // sending message to the back
+           z2_msg->datacount = localcount;
+           BgSendPacket(x,y,z+1,ANYTHREAD,ghostID,SMALL_WORK, sizeof(ghostMsg), (char *)z2_msg);
+           localcount = 0;
+         }
+       }
+      }
+    }
+  }
+}
+
+void reduce(char *info) 
+{
+  // only received by processor 0,0,0
+  int x,y,z;
+  BgGetMyXYZ(&x,&y,&z);
+
+  /*
+  CmiPrintf("reduce procedure being activated on Node (%d,%d,%d)\n",
+          x, y, z);
+  */
+  nodeData *localdata = (nodeData *)BgGetNodeData();
+
+  localdata->done &= ((reduceMsg *)info)->done;
+  delete (reduceMsg *)info;
+  localdata->reduction_count--;
+
+  int numBgX, numBgY, numBgZ;
+  BgGetSize(&numBgX, &numBgY, &numBgZ);
+
+  // if that was the last message
+  if (localdata->reduction_count == 0) {
+    if (!localdata->done) {
+      // more work to be done, so ask every node to exchange ghost regions
+      for (int i=0; i<numBgX;i++) {
+       for (int j=0; j<numBgY;j++) {
+         for (int k=0; k<numBgZ;k++) {
+           exchangeMsg *msg = new exchangeMsg;
+           if ((i == 0) && (j == 0) && (k == 0)) {
+             BgSendLocalPacket(ANYTHREAD,exchangeID, LARGE_WORK, sizeof(exchangeMsg),(char *)msg);
+           } else {
+             BgSendPacket(i,j,k,ANYTHREAD,exchangeID,SMALL_WORK, sizeof(exchangeMsg),(char *)msg);
+           }
+         }
+       }
+      }
+      // resetting the value of the reduction count
+      localdata->reduction_count = numBgX * numBgY * numBgZ;
+      // resetting the completion status of the computation to assume true
+      localdata->done = true; 
+    } else {
+      // instruct all computing nodes to output their individual results
+      for (int i=0; i<numBgX;i++) {
+       for (int j=0; j<numBgY;j++) {
+         for (int k=0; k<numBgZ;k++) {
+           outputMsg *msg = new outputMsg;
+           if ((i == 0) && (j == 0) && (k == 0)) {
+             BgSendLocalPacket(ANYTHREAD,outputID, LARGE_WORK, sizeof(outputMsg), (char *)msg);
+           } else {
+             BgSendPacket(i,j,k,ANYTHREAD,outputID,SMALL_WORK, sizeof(outputMsg), (char *)msg);
+           }
+         }
+       }
+      }
+    }
+  } 
+}
+
+void ghostrecv(char *info)
+{
+  int x,y,z;
+  BgGetMyXYZ(&x,&y,&z);
+
+  /*
+  CmiPrintf("Node (%d,%d,%d) receiving exchange message\n",
+          x, y, z);
+  */
+  nodeData *localdata = (nodeData *)BgGetNodeData();
+  
+  int x_size = localdata->my_x_size;
+  int y_size = localdata->my_y_size;
+  int z_size = localdata->my_z_size;
+
+  int x1_total = localdata->ghost_x1_elements_total;
+  int x2_total = localdata->ghost_x2_elements_total;
+  int y1_total = localdata->ghost_y1_elements_total;
+  int y2_total = localdata->ghost_y2_elements_total;
+  int z1_total = localdata->ghost_z1_elements_total;
+  int z2_total = localdata->ghost_z2_elements_total;
+
+  // determine the source of the ghost region and copy the data to the
+  // appropriate local ghost region.
+  if (((ghostMsg *)info)->source == LEFT) {
+    for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
+      int x_offset = 
+       ((localdata->ghost_x1_elements_current)/(y_size*z_size)) % 1;
+      int y_offset =
+       ((localdata->ghost_x1_elements_current)/z_size) % y_size;
+      int z_offset =
+       (localdata->ghost_x1_elements_current) % z_size;
+      (localdata->ghost_x1)[x_offset][y_offset][z_offset] =
+       ((ghostMsg*)info)->data[count];
+      localdata->ghost_x1_elements_current++;
+    }
+  } else if (((ghostMsg *)info)->source == RIGHT) {
+    for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
+      int x_offset = 
+       ((localdata->ghost_x2_elements_current)/(y_size*z_size)) % 1;
+      int y_offset =
+       ((localdata->ghost_x2_elements_current)/z_size) % y_size;
+      int z_offset =
+       (localdata->ghost_x2_elements_current) % z_size;
+      (localdata->ghost_x2)[x_offset][y_offset][z_offset] =
+       ((ghostMsg*)info)->data[count];
+      localdata->ghost_x2_elements_current++;
+    }
+  } else if (((ghostMsg *)info)->source == BELOW) {
+    for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
+      int x_offset = 
+       ((localdata->ghost_y1_elements_current)/(1*z_size)) % x_size;
+      int y_offset =
+       ((localdata->ghost_y1_elements_current)/z_size) % 1;
+      int z_offset =
+       (localdata->ghost_y1_elements_current) % z_size;
+      (localdata->ghost_y1)[x_offset][y_offset][z_offset] =
+       ((ghostMsg*)info)->data[count];
+      localdata->ghost_y1_elements_current++;
+    }
+  } else if (((ghostMsg *)info)->source == ABOVE) {
+    for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
+      int x_offset = 
+       ((localdata->ghost_y2_elements_current)/(1*z_size)) % x_size;
+      int y_offset =
+       ((localdata->ghost_y2_elements_current)/z_size) % 1;
+      int z_offset =
+       (localdata->ghost_y2_elements_current) % z_size;
+      (localdata->ghost_y2)[x_offset][y_offset][z_offset] =
+       ((ghostMsg*)info)->data[count];
+      localdata->ghost_y2_elements_current++;
+    }
+  } else if (((ghostMsg *)info)->source == FRONT) {
+    for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
+      int x_offset = 
+       ((localdata->ghost_z1_elements_current)/(y_size*1)) % x_size;
+      int y_offset =
+       ((localdata->ghost_z1_elements_current)/1) % y_size;
+      int z_offset =
+       (localdata->ghost_z1_elements_current) % 1;
+      (localdata->ghost_z1)[x_offset][y_offset][z_offset] =
+       ((ghostMsg*)info)->data[count];
+      localdata->ghost_z1_elements_current++;
+    }
+  } else if (((ghostMsg *)info)->source == BACK) {
+    for (int count=0;count<(((ghostMsg*)info)->datacount);count++) {
+      int x_offset = 
+       ((localdata->ghost_z2_elements_current)/(y_size*1)) % x_size;
+      int y_offset =
+       ((localdata->ghost_z2_elements_current)/1) % y_size;
+      int z_offset =
+       (localdata->ghost_z2_elements_current) % 1;
+      (localdata->ghost_z2)[x_offset][y_offset][z_offset] =
+       ((ghostMsg*)info)->data[count];
+      localdata->ghost_z2_elements_current++;
+    }
+  }
+
+  // if that was the last message
+  if ((x1_total == localdata->ghost_x1_elements_current) &&
+      (x2_total == localdata->ghost_x2_elements_current) &&
+      (y1_total == localdata->ghost_y1_elements_current) &&
+      (y2_total == localdata->ghost_y2_elements_current) &&
+      (z1_total == localdata->ghost_z1_elements_current) &&
+      (z2_total == localdata->ghost_z2_elements_current)) {
+    // reset exchange counts
+    localdata->ghost_x1_elements_current = 0;
+    localdata->ghost_x2_elements_current = 0;
+    localdata->ghost_y1_elements_current = 0;
+    localdata->ghost_y2_elements_current = 0;
+    localdata->ghost_z1_elements_current = 0;
+    localdata->ghost_z2_elements_current = 0;
+    
+    computeMsg *msg = new computeMsg;
+    BgSendLocalPacket(ANYTHREAD,computeID, LARGE_WORK, sizeof(computeMsg), (char *)msg); // get to work!
+  }
+  delete (ghostMsg*)info;
+}
+
+void outputData(char *info) {
+  delete (outputMsg *)info;
+  int x,y,z;
+  BgGetMyXYZ(&x,&y,&z);
+
+  /*
+  CmiPrintf("Node (%d,%d,%d) printing data.\n",
+          x, y, z);
+  */
+  nodeData *localdata = (nodeData *)BgGetNodeData();
+  
+  int x_size = localdata->my_x_size;
+  int y_size = localdata->my_y_size;
+  int z_size = localdata->my_z_size;
+
+  CmiPrintf("Final output at Node (%d,%d,%d) with iteration count = %d:\n",
+          x, y, z,
+          localdata->iteration_count);
+  printArray(localdata->maindata,x_size,y_size,z_size);
+
+  if ((x == 0) && (y == 0) && (z == 0)) {
+    BgShutdown();
+  }
+}
+
+void initArray(double target[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+              int x, int y, int z) {
+  for (int i=0; i < x; i++) {
+    for (int j=0; j < y; j++) {
+      for (int k=0; k < z; k++) {
+       target[i][j][k] = (i+j+k)*1.0;
+      }
+    }
+  }
+}
+
+void copyArray(double source[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+              double dest[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+              int x, int y, int z,
+              int offset_x, int offset_y, int offset_z) {
+  for (int i=0; i < x; i++) {
+    for (int j=0; j < y; j++) {
+      for (int k=0; k < z; k++) {
+       dest[i][j][k] = source[i+offset_x][j+offset_y][k+offset_z];
+      }
+    }
+  }
+}
+
+void copyXArray(double source[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+               double dest[1][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+               int x, int y, int z,
+               int offset_x, int offset_y, int offset_z) {
+  for (int i=0; i < x; i++) {
+    for (int j=0; j < y; j++) {
+      for (int k=0; k < z; k++) {
+       dest[i][j][k] = source[i+offset_x][j+offset_y][k+offset_z];
+      }
+    }
+  }
+}
+
+void copyYArray(double source[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+               double dest[MAX_ARRAY_SIZE][1][MAX_ARRAY_SIZE],
+               int x, int y, int z,
+               int offset_x, int offset_y, int offset_z) {
+  for (int i=0; i < x; i++) {
+    for (int j=0; j < y; j++) {
+      for (int k=0; k < z; k++) {
+       dest[i][j][k] = source[i+offset_x][j+offset_y][k+offset_z];
+      }
+    }
+  }
+}
+
+void copyZArray(double source[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+               double dest[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][1],
+               int x, int y, int z,
+               int offset_x, int offset_y, int offset_z) {
+  for (int i=0; i < x; i++) {
+    for (int j=0; j < y; j++) {
+      for (int k=0; k < z; k++) {
+       dest[i][j][k] = source[i+offset_x][j+offset_y][k+offset_z];
+      }
+    }
+  }
+}
+
+void printArray(double target[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE][MAX_ARRAY_SIZE],
+               int x, int y, int z) {
+  for (int i=0; i < x; i++) {
+    for (int j=0; j < y; j++) {
+      for (int k=0; k < z; k++) {
+       CmiPrintf("%f ",target[i][j][k]);
+      }
+      CmiPrintf("\n");
+    }
+  }
+}
+
diff --git a/examples/bigsim/emulator/line.C b/examples/bigsim/emulator/line.C
new file mode 100644 (file)
index 0000000..7622335
--- /dev/null
@@ -0,0 +1,219 @@
+#include <stdlib.h>
+#include "blue.h"
+
+int BCastID ;
+int TimeID  ;
+const int X_NODE = 1;
+const int Y_NODE = 2;
+const int Z_NODE = 3;
+const int MAX_ITERATIONS = 50;
+
+extern "C" void BCastHandler(char*);
+extern "C" void TimeHandler(char*);
+
+class Msg {
+  public:
+    char core[CmiBlueGeneMsgHeaderSizeBytes];
+    int    type;
+    double time;
+
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+//for reporting completion of broadcast
+struct userData {
+  int    expected;
+  int    received;
+  double maxTime;
+  int    iter;
+  int    numIterations;
+  double startTime[MAX_ITERATIONS];
+  double endTime  [MAX_ITERATIONS];
+};
+
+void BgEmulatorInit(int argc, char **argv)
+{
+  if (argc < 7) { 
+    CmiPrintf("Usage: line <x> <y> <z> <numCommTh> <numWorkTh> <numIter> \n"); 
+    BgShutdown();
+  }
+
+  BgSetSize(atoi(argv[1]), atoi(argv[2]), atoi(argv[3]));
+  BgSetNumCommThread(atoi(argv[4]));
+  BgSetNumWorkThread(atoi(argv[5]));
+
+}
+
+void BgNodeStart(int argc, char **argv)
+{
+  BCastID = BgRegisterHandler(BCastHandler);
+  TimeID = BgRegisterHandler(TimeHandler);
+
+  //Initialize NodePrivateData
+  userData* ud = new userData;
+    ud->received = 0;
+    ud->expected = 0;
+    ud->maxTime  = 0.;
+    ud->iter     = 0;
+    ud->numIterations = atoi(argv[6]) ;  //6th Argument is number of Iterations
+
+    int x,y,z;
+    BgGetMyXYZ(&x, &y, &z);
+    int numBgX, numBgY, numBgZ;
+    BgGetSize(&numBgX, &numBgY, &numBgZ);
+
+    if(z==numBgZ-1)
+    {
+       if(y==numBgY-1)
+         ud->expected = 1;
+       else if(y!=0 || (y==0 && x==numBgX-1))
+         ud->expected = 2;
+       else
+         ud->expected = 3;
+    }
+
+  ud->startTime[0] = 0.;       //NOTE: This would cause some error in time of first iteration
+  if(x==0 && y==0 && z==0) 
+  {
+    Msg* msg = new Msg;
+       msg->type = X_NODE;
+    BgSendLocalPacket(ANYTHREAD,BCastID, LARGE_WORK, sizeof(Msg),(char*)msg);
+  }
+      //ckout << "BgNodeInit Done for Node" << bgNode->thisIndex.x <<" "
+      //      << bgNode->thisIndex.y <<" "<< bgNode->thisIndex.z << endl ;
+  
+  BgSetNodeData( (char*)ud);
+}
+
+
+void BCastHandler(char* info) {
+  int i, j, k;
+  BgGetMyXYZ(&i, &j, &k);
+
+  //ckout <<"Broadcast Message Received at Node "<<i<<" "<<j<<" "<<k<<endl;
+
+  Msg* msg = (Msg*)info;
+
+  int numBgX, numBgY, numBgZ;
+  BgGetSize(&numBgX, &numBgY, &numBgZ);
+
+  switch(msg->type){
+  case X_NODE: //possible only on X-axis
+       if(i!=numBgX-1)         //send message further down X with type X_NODE
+         BgSendPacket(i+1, j, k, ANYTHREAD,BCastID, LARGE_WORK, sizeof(Msg), (char*)msg);
+
+       //send message down Y with type Y_NODE
+       msg = new Msg;
+       msg->type = Y_NODE;
+       BgSendPacket(i, j+1, k, ANYTHREAD,BCastID, LARGE_WORK, sizeof(Msg), (char*)msg);
+
+       //send message down Z with type Z_NODE
+       msg = new Msg;
+       msg->type = Z_NODE;
+       BgSendPacket(i, j, k+1, ANYTHREAD,BCastID, LARGE_WORK, sizeof(Msg), (char*)msg);
+
+       break;
+  case Y_NODE:
+
+       if(j!=numBgY-1)         //send message furtherdown Y with type Y_NODE
+         BgSendPacket(i, j+1, k, ANYTHREAD,BCastID, LARGE_WORK, sizeof(Msg), (char*)msg);
+
+       //send message down Z with type Z_NODE
+       msg = new Msg;
+       msg->type = Z_NODE;
+       BgSendPacket(i, j, k+1, ANYTHREAD,BCastID, LARGE_WORK, sizeof(Msg), (char*)msg);
+
+       break;
+  case Z_NODE:
+
+       if(k!=numBgZ-1)         //send message furtherdown Z with type Z_NODE
+         BgSendPacket(i, j, k+1, ANYTHREAD,BCastID, LARGE_WORK, sizeof(Msg), (char*)msg);
+       else                            //send timing information
+       {
+         msg->time = BgGetTime();
+         //ckout <<"Time Message sent from Node "<<i<<" "<<j<<" "<<k
+         //      <<" at time "<<msg->time<<endl;
+         BgSendLocalPacket(ANYTHREAD,TimeID, SMALL_WORK, sizeof(Msg), (char*)msg);
+       }
+       break;
+  }
+}
+
+void TimeHandler(char* info) 
+{
+  int i, j, k;
+  BgGetMyXYZ(&i, &j, &k);
+
+
+  userData *ud = (userData*)BgGetNodeData();
+  Msg* msg = (Msg*)info;
+
+  ud->received++;
+  //ckout <<"Time Message Received at Node "<<i<<" "<<j<<" "<<k
+   //     <<" Received "<<ud->received<<" expected "<<ud->expected<< endl;
+  if(ud->maxTime < msg->time)
+       ud->maxTime = msg->time;
+
+  if(ud->expected == ud->received)
+  {
+       msg->time = ud->maxTime;
+       switch(ud->expected) {
+         case 1:
+               BgSendPacket(i, j-1, k, ANYTHREAD,TimeID, SMALL_WORK, sizeof(Msg), (char*)msg);
+               break;
+         case 2:
+               if(j == 0)
+                 BgSendPacket(i-1, j, k, ANYTHREAD,TimeID, SMALL_WORK, sizeof(Msg), (char*)msg);
+               else
+                 BgSendPacket(i, j-1, k, ANYTHREAD,TimeID, SMALL_WORK, sizeof(Msg), (char*)msg);
+               break;
+         case 3:
+               if(i == 0)
+               {
+                 ud->endTime[ud->iter] = ud->maxTime;
+                 //ckout << "Iteration No " << ud->iter << ", end time " << ud->endTime[ud->iter] << endl;
+
+                 if(ud->iter == ud->numIterations-1)
+                 {
+                       //print result
+                       CmiPrintf("-------------------------------------------------------------------------\n");
+                       CmiPrintf("Iter No:     StartTime               EndTime                 TotalTime \n");
+                       CmiPrintf("-------------------------------------------------------------------------\n");
+                       double avg=0;
+                       for(int i=0; i<ud->numIterations; i++)
+                       {
+                       CmiPrintf("%d               %f               %f               %f\n", i, ud->startTime[i], ud->endTime[i], ud->endTime[i]-ud->startTime[i]);
+
+                         avg += ud->endTime[i]-ud->startTime[i];
+                       }
+                       CmiPrintf("-------------------------------------------------------------------------\n");
+                       CmiPrintf("Average BroadCast Time:                      %f\n", avg/ud->numIterations);
+                       CmiPrintf("-------------------------------------------------------------------------\n");
+                       BgShutdown();
+                       return;
+                 }
+                 ud->iter += 1;
+                 ud->startTime[ud->iter] = BgGetTime();
+                 //ckout << "Iteration No " << ud->iter << ", start time " << ud->startTime[ud->iter] << endl;
+                 msg->type = X_NODE;
+                 BgSendPacket(0, 0, 0, ANYTHREAD,BCastID, LARGE_WORK, sizeof(Msg), (char*)msg);
+               }
+               else
+                 BgSendPacket(i-1, j, k, ANYTHREAD,TimeID, SMALL_WORK, sizeof(Msg), (char*)msg);
+               break;
+
+         default:
+               CmiPrintf("\n\n ERROR: More than expected number of message received at node %d %d %d\n\n", i,j,k);
+               break;
+       }
+        //cleanup for next iteration
+        ud->received = 0;
+        ud->maxTime  = 0;      
+
+  }
+  else 
+       delete msg;
+}
+
+
diff --git a/examples/bigsim/emulator/littleMD/Atom.h b/examples/bigsim/emulator/littleMD/Atom.h
new file mode 100644 (file)
index 0000000..f2869b3
--- /dev/null
@@ -0,0 +1,13 @@
+#ifndef __Atom_h__
+#define __Atom_h__
+
+struct Atom {
+  double m;              // Mass in AMU (= 1.66E-27 kg)
+  double q;              // Charge in units of elemental charge ( = 1.602E-19C)
+  double x, y, z;        // Position in Angstroms ( = 1E-10 m)
+  double vx, vy, vz;     // Velocity in angstroms/fs ( = 1E5m/s)
+  double vhx, vhy, vhz;  // Half-step velocity, used for Verlet integration
+  double fx, fy, fz;     // Force in newtons
+};
+
+#endif 
diff --git a/examples/bigsim/emulator/littleMD/Cell.h b/examples/bigsim/emulator/littleMD/Cell.h
new file mode 100644 (file)
index 0000000..6ec955f
--- /dev/null
@@ -0,0 +1,16 @@
+#ifndef __Cell_h__
+#define __Cell_h__
+
+#include "Atom.h"
+
+struct Cell {
+  int    n_atoms;
+  int    x, y, z;                              // index of the cell
+  double min_x, min_y, min_z;
+  double max_x, max_y, max_z;   // probably unnecessary?
+  Atom*  atoms;                 // array to be allocated ;
+
+  Cell() : atoms (NULL) {}
+};
+
+#endif // __Cell_h__
diff --git a/examples/bigsim/emulator/littleMD/Handlers.C b/examples/bigsim/emulator/littleMD/Handlers.C
new file mode 100644 (file)
index 0000000..218b3e4
--- /dev/null
@@ -0,0 +1,620 @@
+/* File                        : Handlers.C 
+ *
+ * Description : Handlers and Utility functions
+ * Note                        : The program is converse-bluegene version of charm-bluegene version written 
+ *                               previously.
+ ******************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "blue.h"
+#include "bgMD.h"
+
+
+/*******************************************************************************
+               UTILITY FUNCTIONS
+ *****************************************************************************/
+
+void PrintStep(int steps, double pe, double ke)
+{
+  CmiPrintf("Step %d PE: %.4e  KE: %.4e  Total: %.4e\n", steps,pe,ke,pe+ke);
+}
+
+void cellPairToNode(const CellCoord* cellCoord1, const CellCoord* cellCoord2, const SimParams* params, int &x, int &y, int &z)
+{
+  const CellCoord* cellCoord = cellCoord2;
+  
+  if(cellCoord1->x<cellCoord2->x)
+    cellCoord =  cellCoord1;
+  else if(cellCoord1->x==cellCoord2->x)
+  {
+    if(cellCoord1->y<cellCoord2->y)
+      cellCoord =  cellCoord1;
+    else if(cellCoord1->y==cellCoord2->y)
+    {
+      if(cellCoord1->z<cellCoord2->z)
+       cellCoord =  cellCoord1;
+      else if(cellCoord1->z==cellCoord2->z)
+       CmiAbort("Error: Two cells have same coordinates\n");
+    }
+  }
+
+  /* find the cellPair node; cellCoord to node conversion */
+  //more: We can store where the pairs are stored.
+  //more: Write why we are not storing
+  cellToNode(cellCoord, params, x, y, z);
+}
+
+void cellToNode(const CellCoord* cellCoord, const SimParams* params, int &x, int &y, int &z)
+{
+  int sX, sY, sZ; BgGetSize(&sX, &sY, &sZ);
+
+  int cellsNodeX = params->cell_dim_x <= sX ? 1 : params->cell_dim_x / sX;
+  int cellsNodeY = params->cell_dim_y <= sY ? 1 : params->cell_dim_y / sY;
+  int cellsNodeZ = params->cell_dim_z <= sZ ? 1 : params->cell_dim_z / sZ;
+
+  x = cellCoord->x / cellsNodeX;
+  if (x>= sX) x = sX-1;
+  y = cellCoord->y / cellsNodeY;
+  if (y>= sY) y = sY-1;
+  z = cellCoord->z / cellsNodeZ;
+  if (z>= sZ) z = sZ-1;
+}
+
+/* Utility Function: init_cell
+ * Initializes atoms in a cell
+ */
+void init_cell(Cell* this_cell, int xi, int yi, int zi, const SimParams* params)
+{
+  // Set some cell values
+  this_cell->x = xi;
+  this_cell->y = yi;
+  this_cell->z = zi;
+    
+  const double pdx = params->max_x - params->min_x;
+  const double pdy = params->max_y - params->min_y;
+  const double pdz = params->max_z - params->min_z;
+
+  this_cell->min_x = params->min_x + (pdx * xi) / params->cell_dim_x;
+  this_cell->max_x = params->min_x + (pdx * (xi+1)) / params->cell_dim_x;
+  this_cell->min_y = params->min_y + (pdy * yi) / params->cell_dim_y;
+  this_cell->max_y = params->min_y + (pdy * (yi+1)) / params->cell_dim_y;
+  this_cell->min_z = params->min_z + (pdz * zi) / params->cell_dim_z;
+  this_cell->max_z = params->min_z + (pdz * (zi+1)) / params->cell_dim_z;
+
+  // Build up the list of atoms, using a uniform random distribution
+  // Since we are doing each cell separately, we don't know exactly
+  // how many atoms to put in it.  Therefore, we will allocate a
+  // random number such that the total number of atoms on average
+  // would come out right.  We'll do this by pretending to assign
+  // all the atoms, but only actually assigning a 1/# of cells fraction
+  // of them.
+  double* x = new double[params->n_atoms];
+  double* y = new double[params->n_atoms];
+  double* z = new double[params->n_atoms];
+  int atom_count = 0;
+
+  const double prob = 1.0 / (params->cell_dim_x * params->cell_dim_y 
+                           * params->cell_dim_z);
+  const double cdx = this_cell->max_x - this_cell->min_x;
+  const double cdy = this_cell->max_y - this_cell->min_y;
+  const double cdz = this_cell->max_z - this_cell->min_z;
+
+  // Give some seed that is unique to this cell
+  unsigned short int seed[3];
+  seed[0] = (unsigned short)xi;
+  seed[1] = (unsigned short)yi;
+  seed[2] = (unsigned short)zi;
+
+  int i;
+  for(i=0; i < params->n_atoms; i++) {
+    if (erand48(seed) < prob) {
+      x[atom_count] = this_cell->min_x + (cdx * erand48(seed));
+      y[atom_count] = this_cell->min_y + (cdy * erand48(seed));
+      z[atom_count] = this_cell->min_z + (cdz * erand48(seed));
+      atom_count++;
+    }
+  }
+  if(atom_count > 40) atom_count = 40 ; 
+
+  // Allocate the atom array for the cell
+  this_cell->atoms = new Atom[atom_count];
+  this_cell->n_atoms = atom_count;
+
+  // Store the positions into the cells.  Also randomly determine
+  // a mass and charge, and zero out the velocity and force
+  for(i=0;i<atom_count;i++) {
+    Atom* this_atom = &(this_cell->atoms[i]);
+    
+    this_atom->m = 10.0 * erand48(seed);
+    this_atom->q = 5.0 * erand48(seed) - 2.5;
+    this_atom->x = x[i];
+    this_atom->y = y[i];
+    this_atom->z = z[i];
+    this_atom->vx = this_atom->vy = this_atom->vz = 0;
+    this_atom->vhx = this_atom->vhy = this_atom->vhz = 0;
+    this_atom->fx = this_atom->fy = this_atom->fz = 0;
+//CmiPrintf("m:%f, q:%f, x:%f y:%f z:%f\n", this_atom->m, this_atom->q, this_atom->x, this_atom->y, this_atom->z);
+  }
+
+  delete [] x;
+  delete [] y;
+  delete [] z;
+}
+
+/* Utility Function: cell_neighbor
+ * Calculate forces between all atoms of cell1 and cell2
+ * and return potential energy
+ */
+double calc_pair_interactions(Cell* cell1, Cell* cell2, const SimParams* params)
+{
+  double potentialEnergy = 0.0;
+  int i,j;
+  for(i=0;i<cell1->n_atoms;i++)
+    for(j=0;j<cell2->n_atoms;j++) {
+      potentialEnergy += calc_force(&(cell1->atoms[i]), &(cell2->atoms[j]), params);
+    }
+  return potentialEnergy;
+}
+
+/* Utility Function: calc_force
+ * calculate force between two atom1 and atom2
+ */
+double calc_force(Atom *atom1 , Atom* atom2, const SimParams* params)
+{
+  double potentialEnergy = 0.0;
+  // Convert values to SI units
+  const double x1 = atom1->x * 1E-10;
+  const double y1 = atom1->y * 1E-10;
+  const double z1 = atom1->z * 1E-10;
+  const double q1 = atom1->q * 1.602E-19;
+
+  const double x2 = atom2->x * 1E-10;
+  const double y2 = atom2->y * 1E-10;
+  const double z2 = atom2->z * 1E-10;
+  const double q2 = atom2->q * 1.602E-19;
+
+  const double dx = x2-x1;
+  const double dy = y2-y1;
+  const double dz = z2-z1;
+  const double r2 = dx*dx + dy*dy + dz * dz;
+
+  const double csi = params->cutoff * 1E-10;
+  const double cutoff2 = csi * csi;
+
+  if (r2 > cutoff2)  // Outside cutoff, ignore
+    return 0.0;
+
+  const double r = sqrt(r2);
+  
+  const double kq1q2 = 9e9*q1*q2;
+  const double f = kq1q2 / (r2*r);
+  const double fx = f*dx;
+  const double fy = f*dy;
+  const double fz = f*dz;
+
+  atom1->fx += fx;
+  atom2->fx -= fx;
+  atom1->fy += fy;
+  atom2->fy -= fy;
+  atom1->fz += fz;
+  atom2->fz -= fz;
+
+  potentialEnergy -= kq1q2 / r;
+  return potentialEnergy;
+}
+
+/*Utility Function: calc_self_interactions
+ */
+double  calc_self_interactions(Cell* this_cell, const SimParams* params)
+{
+  double potentialEnergy = 0.0;
+  int i,j;
+  for(i=0;i<this_cell->n_atoms;i++)
+    for(j=i+1;j<this_cell->n_atoms;j++) {
+      potentialEnergy += calc_force(&(this_cell->atoms[i]),
+                                   &(this_cell->atoms[j]), params);
+    }
+  return potentialEnergy;
+}
+
+/* Utility Function: update_cell
+ */
+double update_cell(Cell* this_cell, const SimParams* params, bool firstStep)
+{
+  double kineticEnergy = 0.0;
+  // This is a Verlet integrator, which uses the half-step velocity
+  // to calculate new positions.  This is a time-reversible integration
+  // scheme that has some nice numerical properties
+  // The equations are:
+  // v(t) = v(t-1/2) + dt/2 * a(t)
+  // v(t+1/2) = v(t) + dt/2 * a(t)
+  // x(t+1) = x(t) + dt * v(t+1/2)
+  // Notice that at the end of this function, you end up with
+  // v(t) and x(t+1), so to get x and v for the same point in time,
+  // you must save x before it gets updated.
+
+  const double dt = params->step_sz * 1E-15;
+  const double dt2 = dt/2;
+
+  int i;
+  for(i=0;i<this_cell->n_atoms; i++) {
+    Atom* this_atom = &(this_cell->atoms[i]);
+    const double mass = this_atom->m * 1.66e-27;
+    const double dt2m = dt2/mass;
+
+    const double tax = dt2m * this_atom->fx;
+    const double tay = dt2m * this_atom->fy;
+    const double taz = dt2m * this_atom->fz;
+
+    //    convert things to SI
+    double vx = this_atom->vx * 1E-10/1E-15;
+    double vy = this_atom->vy * 1E-10/1E-15;
+    double vz = this_atom->vz * 1E-10/1E-15;
+    double vhx = this_atom->vhx * 1E-10/1E-15;
+    double vhy = this_atom->vhy * 1E-10/1E-15;
+    double vhz = this_atom->vhz * 1E-10/1E-15;
+
+    // For the first step, we are given an initial velocity.
+    // After that, we must calculate it.
+    if (!firstStep) {
+      vx = vhx + tax;
+      vy = vhy + tay;
+      vz = vhz + taz;
+    }
+    //    cout << "F = " << this_atom->fx << ", " << this_atom->fy << endl;
+    //    cout << "A = " << tax << ", " << tay << endl;
+    //    cout << "V = " << vx << ", " << vy << endl;
+    
+    kineticEnergy += 0.5 * mass * (vx*vx+vy*vy+vz*vz);
+
+    vhx = vx + tax;
+    vhy = vy + tay;
+    vhz = vz + taz;
+    
+    this_atom->x += dt * vhx / 1E-10; 
+    this_atom->y += dt * vhy / 1E-10;
+    this_atom->z += dt * vhz / 1E-10;
+
+    // Convert other values back from SI and store
+    this_atom->vx = vx * 1E-15 / 1E-10;
+    this_atom->vy = vy * 1E-15 / 1E-10;
+    this_atom->vz = vz * 1E-15 / 1E-10;
+    this_atom->vhx = vhx * 1E-15 / 1E-10;
+    this_atom->vhy = vhy * 1E-15 / 1E-10;
+    this_atom->vhz = vhz * 1E-15 / 1E-10;
+  }
+
+  for(i=0;i<this_cell->n_atoms; i++)
+  {
+       this_cell->atoms[i].fx = 0;
+       this_cell->atoms[i].fy = 0;
+       this_cell->atoms[i].fz = 0;
+  }
+
+  return kineticEnergy;
+}
+
+/*******************************************************************************
+               HANDLER FUNCTIONS
+ ******************************************************************************/
+
+/* Handler: sendCoord
+ * Sends Coordinate data of each cell on the node to the node which has 
+ * corresponding cell pair
+ */
+void sendCoord(char* info)
+{
+  int x, y, z;
+  BgGetMyXYZ(&x, &y, &z);
+  // CmiPrintf("sendCoord in node [%d,%d,%d]\n", x, y, z);
+
+  UserData* ud = (UserData*)BgGetNodeData();
+
+  StepMsg* msg = (StepMsg*)info;
+  if (msg->step < ud->step) {CmiAbort("error!\n");}
+  if (msg->step > ud->step) {  //TODO is this a bug?: shouldn't it be msg->step == ud->step+1
+    ud->step = msg->step;
+    /* clear up remains of previous step, if there are any */
+    /*
+    CellPairData* clear = NULL;
+    while(ud->cellPairData) {
+      clear = ud->cellPairData;
+      ud->cellPairData = clear->next;
+      delete [] clear->cell1->atoms;
+      delete clear->cell1;
+      delete [] clear->cell2->atoms;
+      delete clear->cell2;
+      delete clear;
+    }
+    */
+  }
+
+  /* For each cell and it's each neighbour on this node, create a CoordinateMsg
+   * containing it's coordinates, it's neighbour's coordinates, and coordinate &
+   * charge info of all it's atoms.
+   */
+  // msg only new once and can be re-used !!!
+  for(int i=0; i<ud->cellCount; i++)
+  {
+    Cell *MyCell =  ud->cellData[i].myCell ;
+
+    ud->cellData[i].myPotEnergy = 0.0 ;
+    ud->cellData[i].myKinEnergy = 0.0 ;
+    ud->cellData[i].countForceReceived = 0 ;
+
+    for(int num = 0; num < ud->cellData[i].neighborCount; num++)
+    {
+      //TODO:scope of optimization: store it and copy later.
+      CoordinateMsg* cmsg = new CoordinateMsg ;
+      cmsg->step = ud->step;
+      cmsg->cellCoord1 = ud->cellData[i].myCoord ;
+      cmsg->cellCoord2 = ud->cellData[i].neighborCoord[num] ;
+      //CmiPrintf("  Analysing cellpair %d %d %d - %d %d %d\n",  cmsg->cellCoord1.x , cmsg->cellCoord1.y , cmsg->cellCoord1.z ,  cmsg->cellCoord2.x , cmsg->cellCoord2.y ,  cmsg->cellCoord2.z);
+
+      cmsg->nElem= MyCell->n_atoms ;
+      for(int l = 0,m=0; l < MyCell->n_atoms; l++)
+      {
+           cmsg->dCordBuff[m++] = MyCell->atoms[l].x ;
+           cmsg->dCordBuff[m++] = MyCell->atoms[l].y ;
+           cmsg->dCordBuff[m++] = MyCell->atoms[l].z ;
+           cmsg->dCordBuff[m++] = MyCell->atoms[l].q ; //charge of atom
+      }
+
+      /* The cellPair is mapped to either node of cell1 or cell2 depending 
+       * on which is closer to 0,0,0 */
+      int a=0, b=0, c=0;
+      cellPairToNode(&cmsg->cellCoord1, &cmsg->cellCoord2, ud->spData, a, b, c);
+      //TODO:scope of optimization: if same node then use addMessage
+      BgSendPacket(a, b, c, ANYTHREAD, storeCoordID, SMALL_WORK, sizeof(CoordinateMsg), (char*)cmsg);
+    }
+  }
+  delete msg;
+}
+
+/* Handler: storeCoord
+ * The cellPair, on receiving CoordinateMsg, stores the cell data to node private data.
+ * When both cells of cellPair have arrived, the forces and energy is computed and sent 
+ * back to both the cells.
+ */
+void storeCoord(char* info)
+{
+  int i,j;
+  int x, y, z;
+  BgGetMyXYZ(&x, &y, &z);
+  // CmiPrintf("storeCoord in node [%d,%d,%d]\n", x, y, z);
+
+  UserData*  ud = (UserData*)BgGetNodeData();
+  CoordinateMsg* msg = (CoordinateMsg*)info;
+  CellCoord cellCoord1= msg->cellCoord1;
+  CellCoord cellCoord2= msg->cellCoord2;
+
+  // search if this pair is already arrived one leg
+  CellPairData* temp = ud->cellPairData;
+  CellPairData* prev;
+  bool found = false;
+  while(temp != NULL && !found)
+  {
+       prev = temp;
+       if(temp->cellCoord1 == msg->cellCoord2 && temp->cellCoord2 == msg->cellCoord1)   
+           found = true;
+       else
+           temp = temp->next;
+  }
+
+  CellPairData* cellPairData = NULL;
+  Cell* cell1 = NULL;
+  Cell *cell2 = NULL;
+  if (!found) {
+         /* Extract Cell from the message */
+         cellPairData = new CellPairData ;
+         cellPairData->cellCoord1 = msg->cellCoord1 ;
+         cellPairData->cellCoord2 = msg->cellCoord2 ;
+         cell1    = new Cell ;
+         cell1->n_atoms = msg->nElem ;
+         cell1->atoms   = new Atom[cell1->n_atoms] ;
+         for(i =0, j=0; i < cell1->n_atoms; i++)
+         {
+               cell1->atoms[i].x = msg->dCordBuff[j++] ;
+               cell1->atoms[i].y = msg->dCordBuff[j++] ;
+               cell1->atoms[i].z = msg->dCordBuff[j++] ;
+               cell1->atoms[i].q = msg->dCordBuff[j++] ;
+               cell1->atoms[i].fx = 0 ; 
+               cell1->atoms[i].fy = 0 ; 
+               cell1->atoms[i].fz = 0 ; 
+         }
+         cellPairData->cell1 = cell1 ;
+         cellPairData->cell2= NULL ;
+         cellPairData->next = ud->cellPairData ;  // cellPairData will be non null since not found1
+         ud->cellPairData   = cellPairData ;
+         return;
+  }
+  else 
+    cell2 = prev->cell1;
+
+  if (cell2 == NULL) CmiAbort("god!\n");
+  if (prev->cell2!= NULL) {
+    CmiPrintf("cellCoord1(step:%d): %d %d %d - %d %d %d\n", msg->step, cellCoord1.x, cellCoord1.y, cellCoord1.z, cellCoord2.x, cellCoord2.y, cellCoord2.z);
+    CmiAbort("unexpected error\n");
+  }
+
+  {
+         /* Extract Cell from the message */
+         cell1    = new Cell ;
+         cell1->n_atoms = msg->nElem ;
+         cell1->atoms   = new Atom[cell1->n_atoms] ;
+         for(i =0, j=0; i < cell1->n_atoms; i++)
+         {
+               cell1->atoms[i].x = msg->dCordBuff[j++] ;
+               cell1->atoms[i].y = msg->dCordBuff[j++] ;
+               cell1->atoms[i].z = msg->dCordBuff[j++] ;
+               cell1->atoms[i].q = msg->dCordBuff[j++] ;
+               cell1->atoms[i].fx = 0 ; 
+               cell1->atoms[i].fy = 0 ; 
+               cell1->atoms[i].fz = 0 ; 
+         }
+    //CmiPrintf("set cellCoord1(step:%d): %d %d %d - %d %d %d\n", msg->step, cellCoord1.x, cellCoord1.y, cellCoord1.z, cellCoord2.x, cellCoord2.y, cellCoord2.z);
+         prev->cell2 = cell1;
+  }
+
+  
+  {
+       double potEnergy = calc_pair_interactions(cell1, cell2, ud->spData) ;
+
+       ForceMsg* forcemsg ;
+       int a, b, c ;
+   
+       forcemsg = new ForceMsg ;
+       forcemsg->cellCoord = cellCoord1 ;
+       forcemsg->nElem=cell1->n_atoms ;
+       forcemsg->potEnergy = potEnergy ;       //<- only count potential once
+       for(i =0, j=0; i < cell1->n_atoms ; i++)
+       {
+         forcemsg->dForceBuff[j++] =  cell1->atoms[i].fx ;
+         forcemsg->dForceBuff[j++] =  cell1->atoms[i].fy ;
+         forcemsg->dForceBuff[j++] =  cell1->atoms[i].fz ;
+       }
+       cellToNode(&forcemsg->cellCoord, ud->spData, a, b, c);
+       BgSendPacket(a, b, c, ANYTHREAD, retrieveForcesID, SMALL_WORK, sizeof(ForceMsg), (char*)forcemsg);
+    
+       forcemsg = new ForceMsg ;
+       forcemsg->cellCoord = cellCoord2 ;
+       forcemsg->nElem=cell2->n_atoms ;
+       forcemsg->potEnergy = 0.0 ;             //<- only count potential once
+       for(i =0, j=0; i < cell2->n_atoms ; i++)
+       {
+         forcemsg->dForceBuff[j++] =  cell2->atoms[i].fx ;
+         forcemsg->dForceBuff[j++] =  cell2->atoms[i].fy ;
+         forcemsg->dForceBuff[j++] =  cell2->atoms[i].fz ;
+       }
+       cellToNode(&forcemsg->cellCoord, ud->spData, a, b, c);
+       BgSendPacket(a, b, c, ANYTHREAD, retrieveForcesID, SMALL_WORK, sizeof(ForceMsg), (char*)forcemsg);
+
+  }
+  delete msg;
+}
+
+/* Handler: retrieveForces
+ * Extracts forces and potential energy from message sent from cellPair and
+ * update atom coordinates in each cell
+ */
+void retrieveForces(char *info)
+{
+  int x, y, z;
+  BgGetMyXYZ(&x, &y, &z);
+  // CmiPrintf("retrieveForces in node [%d,%d,%d]\n", x, y, z);
+
+  UserData*  ud = (UserData*)BgGetNodeData();
+  ForceMsg*  forcemsg = (ForceMsg*)info;
+
+  /* search for which cell on this node has the forces arrived */
+  CellData* cellDataArray = ud->cellData ;
+  CellCoord* cellCoord = &forcemsg->cellCoord; 
+
+  int cellIndex = 0;
+  bool Found = false;
+  while(cellIndex < ud->cellCount && !(Found))
+  {
+       if(cellDataArray[cellIndex].myCoord == *cellCoord)
+          Found = true;
+       else
+          cellIndex++;
+  }
+  
+  if(!Found)
+       CmiAbort("Error> Cell not found in NodePrivateData\n");
+
+  CellData* cellData = &cellDataArray[cellIndex];
+
+  cellData->countForceReceived++ ;
+
+  /* update data of the cell */
+  cellData->myPotEnergy += forcemsg->potEnergy ;
+  Cell *cell = cellData->myCell ;
+
+  for(int j=0, i=0; i<cell->n_atoms; i++)
+  {
+       cell->atoms[i].fx += forcemsg->dForceBuff[j++] ; 
+       cell->atoms[i].fy += forcemsg->dForceBuff[j++] ; 
+       cell->atoms[i].fz += forcemsg->dForceBuff[j++] ;        
+  }
+
+  /* check if all pair computations are done, update coordinates of atoms */
+  if(cellData->countForceReceived==cellData->neighborCount)
+  {
+       /* calculate self interactions */
+       cellData->myPotEnergy += calc_self_interactions(cellData->myCell, ud->spData);
+       cellData->myKinEnergy = update_cell(cell, ud->spData, cellData->firstStep);
+       cellData->firstStep = false;
+
+       /* start reduction */
+       EnergyMsg* energyMsg = new EnergyMsg;
+       energyMsg->potEnergy = cellData->myPotEnergy;
+       energyMsg->kinEnergy = cellData->myKinEnergy;
+
+       ud->lmdData->nreported++;
+       if (ud->lmdData->nreported == ud->cellCount) {
+          CellPairData* clear = NULL;
+          while(ud->cellPairData) {
+            clear = ud->cellPairData;
+            ud->cellPairData = clear->next;
+            delete [] clear->cell1->atoms;
+            delete clear->cell1;
+            delete [] clear->cell2->atoms;
+            delete clear->cell2;
+            delete clear;
+          }
+          ud->lmdData->nreported=0;
+         /*
+          int x,y,z; BgGetMyXYZ(&x, &y, &z);
+          CmiPrintf("report %d %d %d\n", x, y, z);
+          */
+       }
+
+//    CmiPrintf("report %d %d %d\n", cellCoord->x, cellCoord->y, cellCoord->z);
+       BgSendPacket(0, 0, 0, ANYTHREAD, reduceID, LARGE_WORK, sizeof(EnergyMsg), (char*)energyMsg);
+
+  }
+  delete forcemsg;
+}
+
+/* Handler: reduce
+ */
+void reduce(char *info)
+{
+  UserData*  ud = (UserData*)BgGetNodeData();
+  EnergyMsg*  energyMsg = (EnergyMsg*)info;
+
+  ud->lmdData->systemKinEnergy += energyMsg->kinEnergy ;
+  ud->lmdData->systemPotEnergy += energyMsg->potEnergy ;
+  ud->lmdData->numberOfCellsDone++ ;
+
+// CmiPrintf("done: %d %d\n", ud->lmdData->numberOfCellsDone, ud->spData->total_cells);
+  if(ud->lmdData->numberOfCellsDone == ud->spData->total_cells)
+  {
+    ud->lmdData->numberOfCellsDone = 0;
+    ud->lmdData->systemStepsCompleted++ ;
+    PrintStep(ud->lmdData->systemStepsCompleted, 
+              ud->lmdData->systemPotEnergy,
+              ud->lmdData->systemKinEnergy);
+    ud->lmdData->systemPotEnergy = ud->lmdData->systemKinEnergy = 0.0;
+
+    if (ud->lmdData->systemStepsCompleted >= ud->spData->steps)  {
+      CmiPrintf("TIMING: %fs/step\n", (BgGetTime()-startTime)/ud->spData->steps);
+      BgShutdown();
+    }
+
+    /* start next step */
+    int sX, sY, sZ; BgGetSize(&sX, &sY, &sZ);
+    // new StepMsg once and reuse it
+    for (int i=0; i<sX; i++)
+      for (int j=0; j<sY; j++)
+        for (int k=0; k<sZ; k++) {
+          StepMsg *msg = new StepMsg(ud->lmdData->systemStepsCompleted);
+          BgSendPacket(i, j, k, ANYTHREAD, sendCoordID, LARGE_WORK, sizeof(StepMsg), (char*)msg);
+        }
+  }
+  delete energyMsg;
+}
+
diff --git a/examples/bigsim/emulator/littleMD/Makefile b/examples/bigsim/emulator/littleMD/Makefile
new file mode 100644 (file)
index 0000000..8aef87e
--- /dev/null
@@ -0,0 +1,22 @@
+CHARMC=../../../../bin/charmc $(OPTS)
+
+all: bgMD
+
+bgMD: bgMD.o Handlers.o 
+       $(CHARMC) -language bluegene -o bgMD bgMD.o Handlers.o
+
+bgMD.o: bgMD.C bgMD.h Atom.h Cell.h SimParams.h
+       $(CHARMC) -c bgMD.C
+
+Handlers.o: bgMD.h Atom.h Cell.h SimParams.h Handlers.C
+       $(CHARMC) -c Handlers.C
+
+test: jacobi3D maxReduce
+       ./charmrun ./jacobi3D +p2  3 3 3 2 10
+       ./charmrun ./maxReduce +p2 +cth3 +wth10
+
+clean:
+       rm -f core *.cpm.h
+       rm -f TAGS *.o
+       rm -f bgMD
+       rm -f conv-host charmrun
diff --git a/examples/bigsim/emulator/littleMD/SimParams.h b/examples/bigsim/emulator/littleMD/SimParams.h
new file mode 100644 (file)
index 0000000..95d3f45
--- /dev/null
@@ -0,0 +1,61 @@
+#ifndef _SIM_PARAMS_H
+#define _SIM_PARAMS_H
+
+typedef struct {
+  int    steps ;
+  int    n_atoms ;
+  double min_x, min_y, min_z ;
+  double max_x, max_y, max_z ;  // Angstroms
+  int    cell_dim_x, cell_dim_y, cell_dim_z, total_cells ;
+  double cutoff ;               // Angstroms
+  double margin ;               // Angstroms
+  double step_sz ;              // femtoseconds (=1E-15s)
+
+  void setParamsDebug() ;
+  void setParams270K() ;
+  void setParams40K() ;
+} SimParams ;
+
+inline void SimParams::setParamsDebug() {
+  steps = 10 ;
+  n_atoms = 300 ;
+  min_x = 0 ;
+  min_y = 0 ;
+  min_z = 0 ;
+  max_x = 40 ;
+  max_y = 40 ;
+  max_z = 20 ;
+  cutoff = 15 ;
+  margin = 2 ;
+  step_sz = 1.0 ;
+}
+
+inline void SimParams::setParams270K() {
+  steps = 10 ;
+  n_atoms = 270000 ;  // 5000 old
+  min_x = 0 ;
+  min_y = 0 ;
+  min_z = 0 ;
+  max_x = 450 ;  // 300 old
+  max_y = 450 ;  // 300 old
+  max_z = 450 ;  // 300 old
+  cutoff = 15 ;
+  margin = 2 ;
+  step_sz = 1.0 ;
+}
+
+inline void SimParams::setParams40K() {
+  steps = 10 ;
+  n_atoms = 40000 ;  // 5000 old
+  min_x = 0 ;
+  min_y = 0 ;
+  min_z = 0 ;
+  max_x = 105 ;  // 300 old
+  max_y = 105 ;  // 300 old
+  max_z = 105 ;  // 300 old
+  cutoff = 15 ;
+  margin = 2 ;
+  step_sz = 1.0 ;
+}
+
+#endif
diff --git a/examples/bigsim/emulator/littleMD/bgMD.C b/examples/bigsim/emulator/littleMD/bgMD.C
new file mode 100644 (file)
index 0000000..f14a250
--- /dev/null
@@ -0,0 +1,192 @@
+/* File                        : bgMD.C 
+ * Author               : Arun Singla, Neelam Saboo, Joshua Unger,
+ *                        Gengbin Zheng
+ * Version             : 1.0: 3/24/2001
+ *                        2.0: 3/27/2001 - make this an actual working version 
+ *                             on converse bluegene by Gengbin Zheng
+ *
+ * Description : Sample application for Blue Gene emulator (converse version)
+ *                               Prototype Molecular Dynamics: LittleMD
+ * Note                        : The program is converse-bluegene version of 
+ *                        charm-bluegene version written previously.
+*******************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream.h>
+#include "blue.h"
+#include "bgMD.h"
+
+int sendCoordID;
+int storeCoordID;
+int retrieveForcesID;
+int reduceID;
+
+double startTime;
+
+void BgEmulatorInit(int argc, char **argv) 
+{
+  CmiPrintf("Initializing littleMD\n");
+
+  if (argc < 6) {
+    CmiPrintf("Usage: bgMD <x> <y> <z> <numCommTh> <numWorkTh>\n");
+    BgShutdown();
+  }
+  
+  /* set machine configuration */
+  BgSetSize(atoi(argv[1]), atoi(argv[2]), atoi(argv[3]));
+  BgSetNumCommThread(atoi(argv[4]));
+  BgSetNumWorkThread(atoi(argv[5]));
+
+}
+
+void BgNodeStart(int argc, char **argv)
+{
+  /* register handlers */
+  sendCoordID      = BgRegisterHandler(sendCoord);
+  storeCoordID     = BgRegisterHandler(storeCoord);
+  retrieveForcesID = BgRegisterHandler(retrieveForces);
+  reduceID         = BgRegisterHandler(reduce);
+
+  int x, y, z;
+  BgGetMyXYZ(&x, &y, &z);
+  //CmiPrintf("BgNodeStart [%d,%d,%d]\n", x, y, z);
+
+  int sX, sY, sZ;
+  BgGetSize(&sX, &sY, &sZ);
+
+  /* Initialize simulation space  */
+  SimParams *spData = new SimParams;
+  // spData->setParams270K();
+  // spData->setParams40K();
+   spData->setParamsDebug();
+
+  double pdx = spData->max_x - spData->min_x;
+  double pdy = spData->max_y - spData->min_y;
+  double pdz = spData->max_z - spData->min_z;
+
+  spData->cell_dim_x = 1 + (int)(pdx / (spData->cutoff + spData->margin));
+  spData->cell_dim_y = 1 + (int)(pdy / (spData->cutoff + spData->margin));
+  spData->cell_dim_z = 1 + (int)(pdz / (spData->cutoff + spData->margin));
+  spData->total_cells = spData->cell_dim_x * spData->cell_dim_y * spData->cell_dim_z;
+
+  // CmiPrintf("BgNodeStart [%d,%d,%d] SimParams intialized %d %d %d\n", x, y, z, spData->cell_dim_x, spData->cell_dim_y, spData->cell_dim_z);
+
+  // Initialize cells per Blue Gene Node
+  // Find initial(x1, y1, z1) and final(x2, y2, z2) coordinates of cells mapped
+  //   to each node
+  // If the number of nodes exceeds the number of cells in any dimension, 
+  //   make a one-to-one mapping (Some Blue Gene nodes might be ununsed if 
+  //    simulation space is too small)
+  // Find number of cells mapped (x2-x1+1)*(y2-y1+1)*(z2-z1+1) to each node
+  
+  int cellsNodeX = spData->cell_dim_x <= sX ? 1 : spData->cell_dim_x / sX;
+  int cellsNodeY = spData->cell_dim_y <= sY ? 1 : spData->cell_dim_y / sY;
+  int cellsNodeZ = spData->cell_dim_z <= sZ ? 1 : spData->cell_dim_z / sZ;
+
+//CmiPrintf("%d %d %d cellsNodeX:%d cellsNodeY:%d cellsNodeZ:%d\n", spData->cell_dim_x, spData->cell_dim_y, spData->cell_dim_z,cellsNodeX,cellsNodeY,cellsNodeZ);
+
+  int numCellMapped = 0;
+  int x1, y1, z1, x2, y2, z2;
+  if(x<spData->cell_dim_x && y<spData->cell_dim_y && z<spData->cell_dim_z)
+  {
+    x1 = cellsNodeX * x;
+    y1 = cellsNodeY * y;
+    z1 = cellsNodeZ * z;
+    
+    x2 = x==sX-1 ? (spData->cell_dim_x - 1) : (cellsNodeX*(x+1)-1);
+    y2 = y==sY-1 ? (spData->cell_dim_y - 1) : (cellsNodeY*(y+1)-1);
+    z2 = z==sZ-1 ? (spData->cell_dim_z - 1) : (cellsNodeZ*(z+1)-1);
+
+    numCellMapped = (x2-x1+1)*(y2-y1+1)*(z2-z1+1);
+    if(numCellMapped<0) {
+       CmiPrintf("Error> Number of cells mapped to Node [%d,%d,%d] is less than 0\n", x, y, z);
+       CmiAbort("\n");
+    }
+  }
+
+  // CmiPrintf("Mapping to Node [%d,%d,%d], [%d,%d,%d] to [%d,%d,%d] cells %d\n", x, y, z, x1, y1, z1, x2, y2, z2, numCellMapped);
+
+  CellData *cellData = new CellData[numCellMapped];
+
+  /* Store Neighbour Information for each Blue Gene Node, 
+   * If no cells are mapped to this node, don't use it.
+   */
+  if(numCellMapped>0)
+  {
+    int cellIndex = 0;
+    int neighborCount;
+    for(int i=x1; i<=x2; i++)
+     for(int j=y1; j<=y2; j++)
+      for(int k=z1; k<=z2; k++)
+      {
+       neighborCount = 0;
+
+       cellData[cellIndex].myCoord = CellCoord(i,j,k);
+       cellData[cellIndex].firstStep = true;
+       
+       //Initialize the cell
+       cellData[cellIndex].myCell = new Cell;
+       init_cell(cellData[cellIndex].myCell, i, j, k, spData);
+
+       //find the neighbors
+       //more: scope of optimization
+       cellData[cellIndex].neighborCoord = new CellCoord[(2*DISTANCE+1)*(2*DISTANCE+1)*(2*DISTANCE+1)];
+       int nx, ny, nz;
+       for(int l = -DISTANCE; l <= DISTANCE; l++)
+         for(int m = -DISTANCE; m <= DISTANCE; m++)
+          for(int n = -DISTANCE; n <= DISTANCE; n++)
+         {
+          nx = i + l ;
+          ny = j + m ;
+          nz = k + n ;
+
+          if ((nx < spData->cell_dim_x && nx >= 0) && 
+              (ny < spData->cell_dim_y && ny >= 0) &&
+              (nz < spData->cell_dim_z && nz >= 0) &&
+              !(l==0 && m==0 && n==0)) //<- you aren't your own neighbor
+          {
+               CellCoord* cd = &cellData[cellIndex].neighborCoord[neighborCount];
+               cd->x = nx;
+               cd->y = ny;
+               cd->z = nz;
+               neighborCount++; 
+          }
+         }
+
+        cellData[cellIndex].neighborCount = neighborCount;
+        cellData[cellIndex].myPotEnergy = 0.0;
+        cellData[cellIndex].myKinEnergy = 0.0;
+        cellData[cellIndex].countForceReceived = 0;
+
+        cellIndex++;
+    }
+  }
+//  CmiPrintf("BgNodeStart [%d,%d,%d] cells mapped\n", x, y, z);
+
+  /* declare and set node private data */
+  LittleMDData *lmdData = new LittleMDData;
+  lmdData->systemPotEnergy      = 0.0;
+  lmdData->systemKinEnergy      = 0.0;
+  lmdData->numberOfCellsDone    = 0;
+  lmdData->systemStepsCompleted = 0;
+  lmdData->nreported = 0;
+
+  UserData *ud  = new UserData;
+  ud->lmdData   = lmdData;
+  ud->cellData  = cellData;
+  ud->cellCount = numCellMapped;
+  ud->cellPairData = NULL;
+  ud->spData    = spData;
+  ud->step      = 0;
+
+  BgSetNodeData((char*)ud);
+
+  //CmiPrintf("BgNodeStart [%d,%d,%d] creating micro-task\n", x, y, z);
+  /* trigger computation at each node */
+  StepMsg *msg = new StepMsg(0);
+  BgSendLocalPacket(ANYTHREAD, sendCoordID, LARGE_WORK, sizeof(StepMsg), (char*)msg);
+
+//  CmiPrintf("BgNodeStart Done [%d,%d,%d]\n", x, y, z);
+  if (x==0 && y==0 && z==0) startTime = BgGetTime();
+}
diff --git a/examples/bigsim/emulator/littleMD/bgMD.h b/examples/bigsim/emulator/littleMD/bgMD.h
new file mode 100644 (file)
index 0000000..e07d465
--- /dev/null
@@ -0,0 +1,173 @@
+/* File                        : bgMD.h 
+ * Version             : 1.0; March 24, 2001
+ *
+ * Description : Sample application for Blue Gene emulator (converse version)
+ *                               Prototype Molecular Dynamics: LittleMD
+ * Note                        : The program is converse-bluegene version of charm-bluegene version written 
+ *                               previously.
+ ******************************************************************************/
+
+#ifndef __BG_MD_h__
+#define __BG_MD_h__
+
+#include "SimParams.h"
+#include "Cell.h"
+#include "Atom.h"
+
+/* global constants */
+#define DISTANCE                       3   /* 3-away neighbors */
+#define MAX_ATOMS_PER_CELL     50 
+
+/* Handler function Identifiers */
+extern int sendCoordID;
+extern int storeCoordID;
+extern int retrieveForcesID;
+extern int reduceID;
+
+extern double startTime;
+
+//TODO describe each handler function
+/* Handler function prototypes */
+extern "C" void sendCoord(char *);
+extern "C" void storeCoord(char *);
+extern "C" void retrieveForces(char *);
+extern "C" void reduce(char *);
+
+/* forward declarations */
+class  CellCoord;
+
+//TODO describe each utility function
+/* Utilities functions */
+extern "C" void   init_cell(Cell* cell, int xi, int yi, int zi, const SimParams* params);
+extern "C" double calc_self_interactions(Cell* this_cell, const SimParams* params);
+extern "C" double calc_pair_interactions(Cell* cell1, Cell* cell2, const SimParams* params);
+extern "C" double calc_force(Atom *atom1 , Atom* atom2, const SimParams* params);
+extern "C" double update_cell(Cell* this_cell, const SimParams* params, bool firstStep);
+extern "C" void   cellPairToNode(const CellCoord* cellCoord1, const CellCoord* cellCoord2, const SimParams* params, int &x, int &y, int &z);
+extern "C" void   cellToNode(const CellCoord* cellCoord, const SimParams* params, int &x, int &y, int &z);
+
+
+/* forward declarations */
+struct CellPairDataStruct;
+
+typedef unsigned int Byte;
+
+/* Node Private Data and internal structures */
+class CellCoord
+{
+ public:
+  Byte x,y,z;
+  CellCoord(void) {}
+  CellCoord(Byte Nx,Byte Ny,Byte Nz)
+  {
+       x=Nx; y=Ny; z=Nz;
+  }
+  CellCoord(const CellCoord& cell)
+  {
+       this->x = cell.x; this->y = cell.y; this->z = cell.z;
+  }
+  /*
+  CellCoord& operator=(const CellCoord& cell)
+  {
+       this->x = cell.x; this->y = cell.y; this->z = cell.z;
+  }
+  */
+  friend CellCoord operator+(const CellCoord &a,const CellCoord &b) 
+  {
+    return CellCoord(a.x+b.x,a.y+b.y,a.z+b.z);
+  }
+  bool operator==(const CellCoord &a) 
+  {
+    return (x==a.x && y==a.y && z==a.z);
+  }
+};
+
+typedef struct
+{
+  Cell              *myCell;
+  CellCoord  myCoord;
+  int       neighborCount;     //basically a subset of neighbors for optimization
+  CellCoord* neighborCoord;
+  double     myPotEnergy;
+  double     myKinEnergy;
+  int        countForceReceived;
+  bool      firstStep;
+} CellData;
+
+
+typedef struct CellPairDataStruct
+{
+  Cell                              *cell1, *cell2;
+  CellCoord                 cellCoord1;
+  CellCoord                 cellCoord2;
+  struct CellPairDataStruct* next;
+} CellPairData;
+
+typedef struct 
+{
+  double systemPotEnergy;
+  double systemKinEnergy;
+  int    numberOfCellsDone;
+  int    systemStepsCompleted; /* Number of steps completed */
+  int           nreported;
+} LittleMDData;
+
+/* Node Private Data */
+typedef struct
+{
+  LittleMDData* lmdData;
+  int              cellCount;
+  CellData*     cellData;
+  CellPairData* cellPairData;
+  SimParams*    spData;
+  int           step;
+} UserData;
+
+
+
+/* Message Definitions */
+class StepMsg
+{
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  int step;
+  StepMsg(int s) : step(s) {}
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+class CoordinateMsg    // Message to send the data from Cells to Cellpairs
+{
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  CellCoord cellCoord1, cellCoord2;            // Coordinates of Cells
+  int       nElem;                             // Number of atoms in Cell1
+  double    dCordBuff[MAX_ATOMS_PER_CELL*4];   // Coordinates and charge of all the cells in Cell1
+  int       step;
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+class ForceMsg
+{
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  CellCoord cellCoord;
+  int       nElem;
+  double    potEnergy;
+  double    dForceBuff[MAX_ATOMS_PER_CELL*3];
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+class EnergyMsg
+{
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  double potEnergy;
+  double kinEnergy;
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+#endif
diff --git a/examples/bigsim/emulator/maxReduce.C b/examples/bigsim/emulator/maxReduce.C
new file mode 100644 (file)
index 0000000..5a62a3a
--- /dev/null
@@ -0,0 +1,236 @@
+#include <stdlib.h>
+#include "blue.h"
+
+#define A_SIZE 4
+
+#define X_DIM 3
+#define Y_DIM 3
+#define Z_DIM 3
+
+int REDUCE_HANDLER_ID;
+int COMPUTATION_ID;
+
+extern "C" void reduceHandler(char *);
+extern "C" void computeMax(char *);
+
+class ReductionMsg {
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  int max;
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+class ComputeMsg {
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  int dummy;
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+void BgEmulatorInit(int argc, char **argv)
+{
+  BgSetSize(X_DIM, Y_DIM, Z_DIM);
+}
+
+void BgNodeStart(int argc, char **argv) {
+  int x, y, z;
+  BgGetMyXYZ(&x, &y, &z);
+
+  REDUCE_HANDLER_ID = BgRegisterHandler(reduceHandler);
+  COMPUTATION_ID = BgRegisterHandler(computeMax);
+
+  CmiPrintf("Initializing node %d, %d, %d\n", x,y,z);
+
+  ComputeMsg *msg = new ComputeMsg;
+  BgSendLocalPacket(ANYTHREAD, COMPUTATION_ID, LARGE_WORK, sizeof(ComputeMsg), (char *)msg);
+}
+
+void reduceHandler(char *info) {
+  // assumption: THey are initialized to zero?
+  static int max[X_DIM][Y_DIM][Z_DIM];
+  static int num_msg[X_DIM][Y_DIM][Z_DIM];
+
+  int i,j,k;
+  int external_max;
+
+  BgGetMyXYZ(&i,&j,&k);
+  external_max = ((ReductionMsg *)info)->max;
+  delete (ReductionMsg *)info;
+  num_msg[i][j][k]++;
+
+  if ((i == 0) && (j == 0) && (k == 0)) {
+    // master node expects 4 messages:
+    // 1 from itself;
+    // 1 from the i dimension;
+    // 1 from the j dimension; and
+    // 1 from the k dimension
+    if (num_msg[i][j][k] < 4) {
+      // not ready yet, so just find the max
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+    } else {
+      // done. Can report max data after making last comparison
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+      CmiPrintf("The maximal value is %d \n", max[i][j][k]);
+      BgShutdown();
+      return;
+    }
+  } else if ((i == 0) && (j == 0) && (k != Z_DIM - 1)) {
+    // nodes along the k-axis other than the last one expects 4 messages:
+    // 1 from itself;
+    // 1 from the i dimension;
+    // 1 from the j dimension; and
+    // 1 from the k dimension
+    if (num_msg[i][j][k] < 4) {
+      // not ready yet, so just find the max
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+    } else {
+      // done. Forwards max data to node i,j,k-1 after making last comparison
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+      ReductionMsg *msg = new ReductionMsg;
+      msg->max = max[i][j][k];
+      BgSendPacket(i,j,k-1,ANYTHREAD,REDUCE_HANDLER_ID,LARGE_WORK, sizeof(ReductionMsg), (char *)msg);
+    }
+  } else if ((i == 0) && (j == 0) && (k == Z_DIM - 1)) {
+    // the last node along the k-axis expects 3 messages:
+    // 1 from itself;
+    // 1 from the i dimension; and
+    // 1 from the j dimension
+    if (num_msg[i][j][k] < 3) {
+      // not ready yet, so just find the max
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+    } else {
+      // done. Forwards max data to node i,j,k-1 after making last comparison
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+      ReductionMsg *msg = new ReductionMsg;
+      msg->max = max[i][j][k];
+      BgSendPacket(i,j,k-1,ANYTHREAD,REDUCE_HANDLER_ID,LARGE_WORK, sizeof(ReductionMsg), (char *)msg);
+    }
+  } else if ((i == 0) && (j != Y_DIM - 1)) {
+    // for nodes along the j-k plane except for the last and first row of j,
+    // we expect 3 messages:
+    // 1 from itself;
+    // 1 from the i dimension; and
+    // 1 from the j dimension
+    if (num_msg[i][j][k] < 3) {
+      // not ready yet, so just find the max
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+    } else {
+      // done. Forwards max data to node i,j-1,k after making last comparison
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+      ReductionMsg *msg = new ReductionMsg;
+      msg->max = max[i][j][k];
+      BgSendPacket(i,j-1,k,ANYTHREAD,REDUCE_HANDLER_ID,LARGE_WORK, sizeof(ReductionMsg), (char *)msg);
+    }
+  } else if ((i == 0) && (j == Y_DIM - 1)) {
+    // for nodes along the last row of j on the j-k plane,
+    // we expect 2 messages:
+    // 1 from itself;
+    // 1 from the i dimension;
+    if (num_msg[i][j][k] < 2) {
+      // not ready yet, so just find the max
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+    } else {
+      // done. Forwards max data to node i,j-1,k after making last comparison
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+      ReductionMsg *msg = new ReductionMsg;
+      msg->max = max[i][j][k];
+      BgSendPacket(i,j-1,k,ANYTHREAD,REDUCE_HANDLER_ID,LARGE_WORK, sizeof(ReductionMsg), (char *)msg);
+    }
+  } else if (i != X_DIM - 1) {
+    // for nodes anywhere the last row of i,
+    // we expect 2 messages:
+    // 1 from itself;
+    // 1 from the i dimension;
+    if (num_msg[i][j][k] < 2) {
+      // not ready yet, so just find the max
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+    } else {
+      // done. Forwards max data to node i-1,j,k after making last comparison
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+      ReductionMsg *msg = new ReductionMsg;
+      msg->max = max[i][j][k];
+      BgSendPacket(i-1,j,k,ANYTHREAD,REDUCE_HANDLER_ID,LARGE_WORK, sizeof(ReductionMsg), (char *)msg);
+    }
+  } else if (i == X_DIM - 1) {
+    // last row of i, we expect 1 message:
+    // 1 from itself;
+    if (num_msg[i][j][k] < 1) {
+      // not ready yet, so just find the max
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+    } else {
+      // done. Forwards max data to node i-1,j,k after making last comparison
+      if (max[i][j][k] < external_max) {
+       max[i][j][k] = external_max;
+      }
+      ReductionMsg *msg = new ReductionMsg;
+      msg->max = max[i][j][k];
+      BgSendPacket(i-1,j,k,-1,REDUCE_HANDLER_ID,LARGE_WORK, sizeof(ReductionMsg), (char *)msg);
+    }
+  }
+}
+
+void computeMax(char *info) {
+  int A[A_SIZE][A_SIZE];
+  int i, j;
+  int max = 0;
+
+  int x,y,z; // test variables
+  BgGetMyXYZ(&x,&y,&z);
+
+  // Initialize
+  for (i=0;i<A_SIZE;i++) {
+    for (j=0;j<A_SIZE;j++) {
+      A[i][j] = i*j;
+    }
+  }
+
+  CmiPrintf("Finished Initializing %d %d %d!\n",  x , y , z);
+
+  // Find Max
+  for (i=0;i<A_SIZE;i++) {
+    for (j=0;j<A_SIZE;j++) {
+      if (max < A[i][j]) {
+       max = A[i][j];
+      }
+    }
+  }
+
+  CmiPrintf ("Finished Finding Max\n");
+
+  // prepare to reduce
+  ReductionMsg *msg = new ReductionMsg;
+  msg->max = max;
+  BgSendLocalPacket(ANYTHREAD, REDUCE_HANDLER_ID, LARGE_WORK, sizeof(ReductionMsg), (char *)msg);
+
+  CmiPrintf("Sent reduce message to myself with max value %d\n", max);
+}
+
+
diff --git a/examples/bigsim/emulator/octo.C b/examples/bigsim/emulator/octo.C
new file mode 100644 (file)
index 0000000..37bf28b
--- /dev/null
@@ -0,0 +1,850 @@
+#include <stdlib.h>
+#include "blue.h"
+
+int BCastID;
+int ReduceID;
+const int MAX_MESSAGES = 10;
+const int NumBroadcasts = 10;  
+
+extern "C" void BCastOcto(char*);
+extern "C" void ReduceOcto(char*);
+
+class OctoMsg;
+
+void Octo(int, int, int, OctoMsg*, char*);
+void SplitXYZ(int, int, int, OctoMsg*, char*);
+bool HandleSpecial(int, int, int, OctoMsg*, char*);
+// special cases
+void Handle112(int, int, int, OctoMsg*, char*);
+void Handle11X(int, int, int, OctoMsg*, char*);
+void Handle122(int, int, int, OctoMsg*, char*);
+void Handle12X(int, int, int, OctoMsg*, char*);
+void Handle1XX(int, int, int, OctoMsg*, char*);
+void Handle222(int, int, int, OctoMsg*, char*);
+void Handle22X(int, int, int, OctoMsg*, char*);
+void Handle2XX(int, int, int, OctoMsg*, char*);
+
+class OctoMsg { 
+  public:
+    char core[CmiBlueGeneMsgHeaderSizeBytes];
+    int x_min, y_min, z_min;
+    int x_max, y_max, z_max;
+    bool v;  // true if min_coord already visited
+    int sender_x, sender_y, sender_z;
+
+    OctoMsg() { }
+    OctoMsg(int x1, int x2, int y1, int y2, int z1, int z2, bool visit) :
+      x_min   (x1), 
+      x_max   (x2),
+      y_min   (y1),
+      y_max   (y2),
+      z_min   (z1),
+      z_max   (z2),
+      v       (visit)
+    { 
+    }
+
+    /*
+    friend CkOutStream& operator<<(CkOutStream& os, const OctoMsg& msg)
+    {
+      os <<"OctoMsg from " <<msg.sender_x <<" " <<msg.sender_y <<" " 
+         <<msg.sender_z <<"->x(" <<msg.x_min <<"," <<msg.x_max <<")->y(" 
+         <<msg.y_min <<"," <<msg.y_max <<")->z(" <<msg.z_min <<"," 
+         <<msg.z_max <<")";
+      return os;
+    }
+    */
+    void *operator new(size_t s) { return CmiAlloc(s); }
+    void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+class TimeMsg {
+  public :
+    char core[CmiBlueGeneMsgHeaderSizeBytes];
+    double time;
+    TimeMsg(double t) : time(t) { }
+    void *operator new(size_t s) { return CmiAlloc(s); }
+    void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+// for storing results of different tests
+struct TimeRecord 
+{
+  int    test_num;
+  double* max_time;
+  double* start_time;
+
+  TimeRecord() : test_num(0) {
+    max_time = new double[NumBroadcasts];
+    start_time = new double[NumBroadcasts];
+    for (int i=0; i<NumBroadcasts; i++) {
+      max_time[i] = start_time[i] = 0.0;
+    }
+  }
+
+  ~TimeRecord() {
+    delete [] max_time;
+    delete [] start_time;
+  }
+     
+  void Print(char* info) {
+    //print result
+    double average = 0.0;
+    int sizeX, sizeY, sizeZ;
+    BgGetSize(&sizeX, &sizeY, &sizeZ);
+    int numComm = BgGetNumCommThread();
+    int numWork = BgGetNumWorkThread();
+
+    CmiPrintf("\nResults for %d by %d by %d with %d comm %d work\n\n",
+              sizeX, sizeY, sizeZ, numComm, numWork);
+    CmiPrintf("-------------------------------------------------------------\n"
+              "Iter No:    StartTime           EndTime          TotalTime   \n"
+              "-------------------------------------------------------------\n");
+    for (int i=0; i<NumBroadcasts; i++) {
+      CmiPrintf("    %d         %f               %f           %f\n",
+                i, start_time[i], max_time[i], max_time[i] - start_time[i]);
+      average += max_time[i] - start_time[i];
+    }
+    CmiPrintf("-------------------------------------------------------------\n"
+              "Average BroadCast Time:                     %f\n"
+              "-------------------------------------------------------------\n",
+              average/NumBroadcasts);
+    BgShutdown();
+  }
+};
+
+struct OctoData 
+{
+  OctoMsg     messages[MAX_MESSAGES];
+  int         dest_x[MAX_MESSAGES];
+  int         dest_y[MAX_MESSAGES];
+  int         dest_z[MAX_MESSAGES];
+  int         root_x, root_y, root_z;
+  int         parent_x, parent_y, parent_z;
+  TimeRecord* record;
+  
+  OctoData(int x, int y, int z): 
+    root_x       (x),
+    root_y       (y),
+    root_z       (z),
+    record       (NULL)
+ { }
+};
+
+BnvStaticDeclare(int, num_messages);
+BnvStaticDeclare(double, max_time);
+BnvStaticDeclare(int, reduce_count);
+
+void BgEmulatorInit(int argc, char** argv)
+{
+  if (argc < 6) { 
+    CmiPrintf("Usage: octo <x> <y> <z> <numComm> <numWork> \n"); 
+    BgShutdown();
+  }
+    
+  BgSetSize(atoi(argv[1]), atoi(argv[2]), atoi(argv[3]));
+  BgSetNumCommThread(atoi(argv[4]));
+  BgSetNumWorkThread(atoi(argv[5]));
+
+}
+
+void BgWorkThreadInit(int argc, char** argv)
+{
+}
+
+void BgNodeStart(int argc, char** argv)
+{
+  BCastID = BgRegisterHandler(BCastOcto);
+  ReduceID  = BgRegisterHandler(ReduceOcto);
+
+  int x, y, z;
+  BgGetMyXYZ(&x, &y, &z);
+  int numBgX, numBgY, numBgZ;
+  BgGetSize(&numBgX, &numBgY, &numBgZ);
+
+  BnvInitialize(int, num_messages);
+  BnvAccess(num_messages) = -1;
+  BnvInitialize(double, max_time);
+  BnvInitialize(int, reduce_count);
+  int center_x = (numBgX == 2) ? 0 : numBgX/2;
+  int center_y = (numBgY == 2) ? 0 : numBgY/2;
+  int center_z = (numBgZ == 2) ? 0 : numBgZ/2;
+  OctoData* data = new OctoData(center_x, center_y, center_z);
+  BgSetNodeData((char*)data);
+
+  // check to see if center node
+  if (x == center_x && y == center_y && z == center_z) 
+  {
+    data->record = new TimeRecord();
+    OctoMsg* msg = new OctoMsg(0, numBgX-1, 0, numBgY-1, 0, numBgZ-1, false);
+    BgSendLocalPacket(ANYTHREAD, BCastID, LARGE_WORK, sizeof(OctoMsg), 
+                      (char*)msg);
+  }
+}
+
+void BCastOcto(char* info) {
+  int x, y, z;
+  BgGetMyXYZ(&x, &y, &z);
+
+  OctoMsg* m = (OctoMsg*)info;
+  OctoData* d = (OctoData*)BgGetNodeData();
+  
+  // CmiPrintf("BCastOcto at %d %d %d range x=(%d,%d) y=(%d,%d) z=(%d,%d)\n", 
+  //           x, y, z, m->x_min, m->x_max, m->y_min, m->y_max, m->z_min, 
+  //           m->z_max);
+
+  if (x == d->root_x && y == d->root_y && z == d->root_z) {
+    if (d->record == NULL) {
+      CmiPrintf("Error, root node should have timing info\n");
+      BgShutdown();
+      return;
+    }
+    TimeRecord* r = d->record;
+    r->start_time[r->test_num] = BgGetTime();
+    // CmiPrintf("Starting new broadcast at %f\n", r->start_time[r->test_num]);
+  }
+
+  // if first time through, calculate correct nodes to send
+  if (BnvAccess(num_messages) == -1) { 
+    d->parent_x = m->sender_x;
+    d->parent_y = m->sender_y;
+    d->parent_z = m->sender_z;
+  
+    SplitXYZ(x, y, z, m, info); 
+  }
+
+  // once calculated, send to nodes
+  if (BnvAccess(num_messages) <= 0) {
+    // this is a terminal node, so start reduction
+    TimeMsg* parent = new TimeMsg(BgGetTime());
+    // CmiPrintf("  Terminal node, time %f to %d %d %d\n", 
+    //           parent->time, d->parent_x, d->parent_y, d->parent_z);
+    BgSendPacket(d->parent_x, d->parent_y, d->parent_z, 
+                 ANYTHREAD, ReduceID, SMALL_WORK, sizeof(TimeMsg), 
+                 (char*)parent);
+  }
+  else {
+    // CmiPrintf("  There are %d messages\n", d->num_messages);
+    for (int i=0; i<BnvAccess(num_messages); i++) {
+      OctoMsg* msg = new OctoMsg(d->messages[i]);
+      // CmiPrintf("    Sending message %d/%d to %d %d %d\n", i+1, 
+      //           d->num_messages, d->dest_x[i], d->dest_y[i], d->dest_z[i]);
+      BgSendPacket(d->dest_x[i], d->dest_y[i], d->dest_z[i], 
+                   ANYTHREAD, BCastID, SMALL_WORK, sizeof(OctoMsg), (char*)msg);
+    }
+  }
+}
+
+void ReduceOcto(char* info) 
+{
+  int x, y, z;
+  BgGetMyXYZ(&x, &y, &z);
+
+  OctoData* d = (OctoData*)BgGetNodeData();
+  TimeMsg* msg = (TimeMsg*)info;
+  BnvAccess(reduce_count)++;
+  // CmiPrintf("ReduceOcto at node %d %d %d message %d/%d with time %f "
+  //           "max time is %f\n", x, y, z, d->reduce_count, d->num_messages,
+  //           msg->time, d->max_time);
+  if (msg->time > BnvAccess(max_time)) { BnvAccess(max_time) = msg->time; }
+
+  // check to see if all children have been heard from
+  if (BnvAccess(reduce_count) >= BnvAccess(num_messages)) {
+    BnvAccess(reduce_count) = 0;
+
+    if (x == d->root_x && y == d->root_y && z == d->root_z) {
+      // this is the root node
+      TimeRecord* r = d->record;
+      r->max_time[r->test_num] = BnvAccess(max_time);
+      r->test_num++;
+      if (r->test_num < NumBroadcasts) {
+        // start bcast all over again
+        int numBgX, numBgY, numBgZ;
+        BgGetSize(&numBgX, &numBgY, &numBgZ);
+        OctoMsg* start = 
+          new OctoMsg(0, numBgX-1, 0, numBgY-1, 0, numBgZ-1, false);
+        BgSendLocalPacket(ANYTHREAD, BCastID, SMALL_WORK, sizeof(OctoMsg), 
+                          (char*)start);
+      }
+      else {
+        // print results and quit
+        r->Print(info);
+        BgShutdown();
+        return;
+      }
+    }
+    else {
+      // this is not the root node
+      TimeMsg* parent = new TimeMsg(BnvAccess(max_time));
+      BgSendPacket(d->parent_x, d->parent_y, d->parent_z, ANYTHREAD, ReduceID,
+                   SMALL_WORK, sizeof(TimeMsg), (char*)parent);
+      BnvAccess(max_time) = 0;
+    }
+  }
+}
+
+void SplitXYZ
+(
+  int         x,     // x pos of node
+  int         y,     // y pos of node
+  int         z,     // z pos of node
+  OctoMsg*    m,     // message with limit data
+  char* info
+)
+{
+  // check terminal case
+  if (x == m->x_min && x == m->x_max &&
+      y == m->y_min && y == m->y_max &&
+      z == m->z_min && z == m->z_max)
+  {
+    return;
+  }
+  
+  if (HandleSpecial(x, y, z, m, info)) { return; }
+
+  int del_x = m->x_max - m->x_min;
+  int del_y = m->y_max - m->y_min;
+  int del_z = m->z_max - m->z_min;
+
+  // calculate new midpoints
+  int x_lo = (m->x_min + x - 1)/2;  int x_hi = (m->x_max + x)/2;
+  int y_lo = (m->y_min + y - 1)/2;  int y_hi = (m->y_max + y)/2;
+  int z_lo = (m->z_min + z - 1)/2;  int z_hi = (m->z_max + z)/2;
+  
+  OctoMsg* new_msg = NULL;
+
+  // x_lo, y_lo, z_lo
+  new_msg = new OctoMsg(m->x_min, x-1, m->y_min, y-1, m->z_min, z-1, m->v);
+  if (x_lo == m->x_min && y_lo == m->y_min && z_lo == m->z_min && m->v) {
+    SplitXYZ(x_lo, y_lo, z_lo, new_msg, info);
+    delete new_msg;
+  }
+  else { Octo(x_lo, y_lo, z_lo, new_msg, info); }
+  // x_lo, y_lo, z_hi
+  new_msg = new OctoMsg(m->x_min, x-1, m->y_min, y-1, z, m->z_max, false);
+  Octo(x_lo, y_lo, z_hi, new_msg, info);
+  // x_lo, y_hi, z_lo
+  new_msg = new OctoMsg(m->x_min, x-1, y, m->y_max, m->z_min, z-1, false);
+  Octo(x_lo, y_hi, z_lo, new_msg, info);
+  // x_lo, y_hi, z_hi
+  new_msg = new OctoMsg(m->x_min, x-1, y, m->y_max, z, m->z_max, false);
+  Octo(x_lo, y_hi, z_hi, new_msg, info);
+  // x_hi, y_lo, z_lo
+  new_msg = new OctoMsg(x, m->x_max, m->y_min, y-1, m->z_min, z-1, false);
+  Octo(x_hi, y_lo, z_lo, new_msg, info);
+  // x_hi, y_lo, z_hi
+  new_msg = new OctoMsg(x, m->x_max, m->y_min, y-1, z, m->z_max, false);
+  Octo(x_hi, y_lo, z_hi, new_msg, info);
+  // x_hi, y_hi, z_lo
+  new_msg = new OctoMsg(x, m->x_max, y, m->y_max, m->z_min, z-1, false);
+  Octo(x_hi, y_hi, z_lo, new_msg, info);
+  // x_hi, y_hi, z_hi (x, y, z) 
+  new_msg = new OctoMsg(x, m->x_max, y, m->y_max, z, m->z_max, true);
+  if (x_hi == x && y_hi == y && z_hi == z) {
+    // just call local function instead of passing message
+    SplitXYZ(x_hi, y_hi, z_hi, new_msg, info);
+    delete new_msg;
+  }
+  else { Octo(x_hi, y_hi, z_hi, new_msg, info); }
+}
+
+void Octo
+(
+  int         x,     // x pos of node
+  int         y,     // y pos of node
+  int         z,     // z pos of node
+  OctoMsg*    m,     // message with limit data
+  char* info
+)
+{
+  OctoData* d = (OctoData*)BgGetNodeData();
+  BgGetMyXYZ(&m->sender_x, &m->sender_y, &m->sender_z);
+
+  // save the message if desired
+  if (BnvAccess(num_messages) == -1) { BnvAccess(num_messages)++; }
+  if (BnvAccess(num_messages) < MAX_MESSAGES) {
+    int i = BnvAccess(num_messages);
+    d->dest_x[i] = x;
+    d->dest_y[i] = y;
+    d->dest_z[i] = z;
+    d->messages[i] = *m;
+    BnvAccess(num_messages)++;
+  }
+  else { 
+    CmiPrintf("ERR: No room to write message on node %d %d %d\n",
+              m->sender_x, m->sender_y, m->sender_z);
+    BgShutdown();
+    return;
+  }
+}
+
+bool HandleSpecial
+(
+  int         x,     // x pos of node
+  int         y,     // y pos of node
+  int         z,     // z pos of node
+  OctoMsg*    m,     // message with limit data
+  char* info
+)
+{
+  int del_x = m->x_max - m->x_min;
+  int del_y = m->y_max - m->y_min;
+  int del_z = m->z_max - m->z_min;
+
+  if (del_x >= 2 && del_y >= 2 && del_z >= 2) {
+    return false;
+  }
+
+  if ((del_x == 1 && del_y == 0 && del_z == 0) ||
+      (del_x == 0 && del_y == 1 && del_z == 0) ||
+      (del_x == 0 && del_y == 0 && del_z == 1))
+  {
+    Handle112(x, y, z, m, info);
+  }
+  else if ((del_x >  1 && del_y == 0 && del_z == 0) ||
+           (del_x == 0 && del_y  > 1 && del_z == 0) ||
+           (del_x == 0 && del_y == 0 && del_z  > 1))
+  {
+    Handle11X(x, y, z, m, info);
+  }
+  else if ((del_x == 1 && del_y == 1 && del_z == 0) ||
+           (del_x == 1 && del_y == 0 && del_z == 1) ||
+           (del_x == 0 && del_y == 1 && del_z == 1))
+  {
+    Handle122(x, y, z, m, info);
+  }
+  else if ((del_x == 0 && del_y == 1 && del_z  > 1) ||
+           (del_x == 0 && del_y  > 1 && del_z == 1) ||
+           (del_x == 1 && del_y == 0 && del_z  > 1) ||
+           (del_x == 1 && del_y  > 1 && del_z == 0) ||
+           (del_x  > 1 && del_y == 0 && del_z == 1) ||
+           (del_x  > 1 && del_y == 1 && del_z == 0)) 
+  {
+    Handle12X(x, y, z, m, info);
+  }
+  else if ((del_x == 0 && del_y  > 1 && del_z  > 1) ||
+           (del_x  > 1 && del_y == 0 && del_z  > 1) ||
+           (del_x  > 1 && del_y  > 1 && del_z == 0))
+  {
+    Handle1XX(x, y, z, m, info);
+  }
+  else if (del_x == 1 && del_y == 1 && del_z == 1) {
+    Handle222(x, y, z, m, info);
+  }
+  else if ((del_x == 1 && del_y == 1 && del_z  > 1) ||
+           (del_x == 1 && del_y  > 1 && del_z == 1) ||
+           (del_x  > 1 && del_y == 1 && del_z == 1)) 
+  {
+    Handle22X(x, y, z, m, info);
+  }
+  else if ((del_x == 1 && del_y  > 1 && del_z  > 1) ||
+           (del_x  > 1 && del_y == 1 && del_z  > 1) ||
+           (del_x  > 1 && del_y  > 1 && del_z == 1)) 
+  {
+    Handle2XX(x, y, z, m, info);
+  }
+  else {
+    CmiPrintf("ERR: COULDN'T HANDLE %d x %d x %d\n", del_x+1, del_y+1, del_z+1);
+    BgShutdown();
+    return false;
+  }
+  
+  return true;
+}
+
+void Handle112(int x, int y, int z, OctoMsg* m, char* info)
+{
+  // 2x1x1 -> only 2nd x hasn't been touched
+  // 1x2x1 -> only 2nd y hasn't been touched
+  // 1x1x2 -> only 2nd z hasn't been touched
+
+  int del_x = m->x_max - m->x_min;
+  int del_y = m->y_max - m->y_min;
+  int del_z = m->z_max - m->z_min;
+
+  OctoMsg* new_msg = 
+    new OctoMsg(x+del_x, x+del_x, y+del_y, y+del_y, z+del_z, z+del_z, false);
+  Octo(x+del_x, y+del_y, z+del_z, new_msg, info);
+}
+
+void Handle11X(int x, int y, int z, OctoMsg* m, char* info)
+{
+  int x_mod = ((m->x_max - m->x_min) > 1) ? 1 : 0;
+  int y_mod = ((m->y_max - m->y_min) > 1) ? 1 : 0;
+  int z_mod = ((m->z_max - m->z_min) > 1) ? 1 : 0; 
+    
+  // calculate new midpoints
+  int x_lo = (m->x_min + x - x_mod)/2;  int x_hi = (m->x_max + x)/2;
+  int y_lo = (m->y_min + y - y_mod)/2;  int y_hi = (m->y_max + y)/2;
+  int z_lo = (m->z_min + z - z_mod)/2;  int z_hi = (m->z_max + z)/2;
+
+  // hit high node
+  OctoMsg* new_msg = new OctoMsg(x, m->x_max, y, m->y_max, z, m->z_max, true);
+  if (x == x_hi && y == y_hi && z == z_hi) {
+    SplitXYZ(x_hi, y_hi, z_hi, new_msg, info);
+    delete new_msg;
+  }
+  else { Octo(x_hi, y_hi, z_hi, new_msg, info); }
+
+  // hit low node if it hasn't been hit yet
+  int new_x_width = (x - x_mod - m->x_min + 1);
+  int new_y_width = (y - y_mod - m->y_min + 1);
+  int new_z_width = (z - z_mod - m->z_min + 1);
+  
+  new_msg = new OctoMsg(m->x_min, x-x_mod, m->y_min, y-y_mod, 
+                        m->z_min, z-z_mod, m->v);
+
+  // error handle all different cases
+  if (new_x_width * new_y_width * new_z_width == 1) {
+    if (!m->v) { Octo(x_lo, y_lo, z_lo, new_msg, info); }
+    else { delete new_msg; }
+  }
+  else if (x == x_lo && y == y_lo && z == z_lo) {
+    SplitXYZ(x_lo, y_lo, z_lo, new_msg, info);
+    delete new_msg;
+  }
+  else if (new_x_width * new_y_width * new_z_width != 0) { 
+    if (x_lo == m->x_min && y_lo == m->y_min && z_lo == m->z_min && m->v) {
+      SplitXYZ(x_lo, y_lo, z_lo, new_msg, info);
+      delete new_msg;
+    }
+    else { Octo(x_lo, y_lo, z_lo, new_msg, info); }
+  }
+}
+
+void Handle122(int x, int y, int z, OctoMsg* m, char* info)
+{
+  int del_x = m->x_max - m->x_min;
+  int del_y = m->y_max - m->y_min;
+  int del_z = m->z_max - m->z_min;
+
+  OctoMsg* new_msg = NULL;
+  
+  // 2x2x1 -> split along x axiz and handle it
+  // 2x1x2 -> split along x axis and handle it
+  if (del_z == 0 || del_y == 0)
+  {
+    new_msg = new OctoMsg(x, x, m->y_min, m->y_max, m->z_min, m->z_max, true);
+    SplitXYZ(x, y, z, new_msg, info); 
+    delete new_msg;
+
+    new_msg = 
+      new OctoMsg(x+1, x+1, m->y_min, m->y_max, m->z_min, m->z_max, false);
+    Octo(x+1, y, z, new_msg, info);
+  }
+  // 1x2x2 -> split along y axis and handle it
+  else if (del_x == 0) {
+    new_msg = new OctoMsg(x, x, y, y, m->z_min, m->z_max, true);
+    SplitXYZ(x, y, z, new_msg, info);
+    delete new_msg;
+
+    new_msg = new OctoMsg(x, x, y+1, y+1, m->z_min, m->z_max, false);
+    Octo(x, y+1, z, new_msg, info);
+  }
+}
+
+void Handle12X(int x, int y, int z, OctoMsg* m, char* info)
+{
+  int x_mod = ((m->x_max - m->x_min) > 1) ? 1 : 0;
+  int y_mod = ((m->y_max - m->y_min) > 1) ? 1 : 0;
+  int z_mod = ((m->z_max - m->z_min) > 1) ? 1 : 0; 
+    
+  // calculate new midpoints
+  int x_lo = (m->x_min + x - x_mod)/2;  int x_hi = (m->x_max + x)/2;
+  int y_lo = (m->y_min + y - y_mod)/2;  int y_hi = (m->y_max + y)/2;
+  int z_lo = (m->z_min + z - z_mod)/2;  int z_hi = (m->z_max + z)/2;
+
+  // first split high side
+  OctoMsg* new_msg = new OctoMsg(x, m->x_max, y, m->y_max, z, m->z_max, true);
+  if (x == x_hi && y == y_hi && z == z_hi) {
+    SplitXYZ(x_hi, y_hi, z_hi, new_msg, info);
+    delete new_msg;
+  }
+  else { Octo(x_hi, y_hi, z_hi, new_msg, info); }
+
+  int del_x = m->x_max - m->x_min;
+  int del_y = m->y_max - m->y_min;
+  int del_z = m->z_max - m->z_min;
+
+  // hit low node if it hasn't been hit yet
+  int new_x_width = del_x > 1 ? (x - x_mod - m->x_min + 1) : del_x + 1;
+  int new_y_width = del_y > 1 ? (y - y_mod - m->y_min + 1) : del_y + 1;
+  int new_z_width = del_z > 1 ? (z - z_mod - m->z_min + 1) : del_z + 1;
+  
+  // adjust mods
+  x_mod = ((m->x_max - m->x_min) == 1) ? -1 : x_mod;
+  y_mod = ((m->y_max - m->y_min) == 1) ? -1 : y_mod;
+  z_mod = ((m->z_max - m->z_min) == 1) ? -1 : z_mod; 
+  new_msg = new OctoMsg(m->x_min, x-x_mod, m->y_min, y-y_mod, 
+                        m->z_min, z-z_mod, m->v);
+
+  // error handle all different cases
+  if (new_x_width * new_y_width * new_z_width == 2) {
+    if (!m->v) { Octo(x_lo, y_lo, z_lo, new_msg, info); }
+    else { 
+      SplitXYZ(x_lo, y_lo, z_lo, new_msg, info);
+      delete new_msg; 
+    }
+  }
+  else if (x == x_lo && y == y_lo && z == z_lo) {
+    SplitXYZ(x_lo, y_lo, z_lo, new_msg, info);
+    delete new_msg;
+  }
+  else if (new_x_width * new_y_width * new_z_width != 0) { 
+    if (x_lo == m->x_min && y_lo == m->y_min && z_lo == m->z_min && m->v) {
+      SplitXYZ(x_lo, y_lo, z_lo, new_msg, info);
+      delete new_msg;
+    }
+    else { Octo(x_lo, y_lo, z_lo, new_msg, info); }
+  }
+}
+
+void Handle1XX(int x, int y, int z, OctoMsg* m, char* info)
+{
+  int del_x = m->x_max - m->x_min;
+  int del_y = m->y_max - m->y_min;
+  int del_z = m->z_max - m->z_min;
+
+  int x_mod = ((m->x_max - m->x_min) > 1) ? 1 : 0;
+  int y_mod = ((m->y_max - m->y_min) > 1) ? 1 : 0;
+  int z_mod = ((m->z_max - m->z_min) > 1) ? 1 : 0; 
+
+  // calculate new midpoints
+  int x_lo = (m->x_min + x - x_mod)/2;  int x_hi = (m->x_max + x)/2;
+  int y_lo = (m->y_min + y - y_mod)/2;  int y_hi = (m->y_max + y)/2;
+  int z_lo = (m->z_min + z - z_mod)/2;  int z_hi = (m->z_max + z)/2;
+
+  OctoMsg* new_msg = 
+    new OctoMsg(x, m->x_max, y, m->y_max, z, m->z_max, true);
+  if (x == x_hi && y == y_hi && z == z_hi) {
+    SplitXYZ(x_hi, y_hi, z_hi, new_msg, info);
+    delete new_msg;
+  }
+  else { Octo(x_hi, y_hi, z_hi, new_msg, info); }
+    
+  // hit low node if it hasn't been hit yet
+  int new_x_width = (x-x_mod - m->x_min + 1);
+  int new_y_width = (y-y_mod - m->y_min + 1);
+  int new_z_width = (z-z_mod - m->z_min + 1);
+  
+  new_msg = new OctoMsg(m->x_min, x-x_mod, m->y_min, y-y_mod, 
+                        m->z_min, z-z_mod, m->v);
+
+  // error handle all different cases
+  if (new_x_width * new_y_width * new_z_width == 1) {
+    if (!m->v) { Octo(x_lo, y_lo, z_lo, new_msg, info); }
+    else { delete new_msg; }
+  }
+  else if (x == x_lo && y == y_lo && z == z_lo) {
+    SplitXYZ(x_lo, y_lo, z_lo, new_msg, info);
+    delete new_msg;
+  }
+  else { 
+    if (x_lo == m->x_min && y_lo == m->y_min && z_lo == m->z_min && m->v) {
+      SplitXYZ(x_lo, y_lo, z_lo, new_msg, info);
+      delete new_msg;
+    }
+    else { Octo(x_lo, y_lo, z_lo, new_msg, info); }
+  }
+
+  // calculate new mods
+  x_mod = (x - m->x_min - 1);  x_mod = x_mod > 0 ? x_mod : 0;
+  y_mod = (y - m->y_min - 1);  y_mod = y_mod > 0 ? y_mod : 0;
+  z_mod = (z - m->z_min - 1);  z_mod = z_mod > 0 ? z_mod : 0;
+
+  // 1xXxX
+  if (del_x == 0) {
+    new_msg = 
+      new OctoMsg(x, x, y, m->y_max, m->z_min, m->z_min + z_mod, false);
+    Octo(x, y_hi, z_lo, new_msg, info);
+    
+    new_msg = 
+      new OctoMsg(x, x, m->y_min, m->y_min + y_mod, z, m->z_max, false);
+    Octo(x, y_lo, z_hi, new_msg, info);
+  }
+  // Xx1xX
+  else if (del_y == 0) {
+    new_msg = 
+      new OctoMsg(x, m->x_max, y, y, m->z_min, m->z_min + z_mod, false);
+    Octo(x_hi, y, z_lo, new_msg, info);
+    
+    new_msg = 
+      new OctoMsg(m->x_min, m->x_min + x_mod, y, y, z, m->z_max, false);
+    Octo(x_lo, y, z_hi, new_msg, info);
+  }
+  // XxXx1
+  else { // (del_z == 0)
+    new_msg = 
+      new OctoMsg(x, m->x_max, m->y_min, m->y_min + y_mod, z, z, false);
+    Octo(x_hi, y_lo, z, new_msg, info);
+    
+    new_msg = 
+      new OctoMsg(m->x_min, m->x_min + x_mod, y, m->y_max, z, z, false);
+    Octo(x_lo, y_hi, z, new_msg, info);
+  }
+}
+
+void Handle222(int x, int y, int z, OctoMsg* m, char* info)
+{
+  int del_x = m->x_max - m->x_min;
+  int del_y = m->y_max - m->y_min;
+  int del_z = m->z_max - m->z_min;
+
+  // 2x2x2 -> split along z axis and handle it
+  OctoMsg* new_msg = 
+    new OctoMsg(m->x_min, m->x_max, m->y_min, m->y_max, z, z, true);
+  SplitXYZ(x, y, z, new_msg, info);
+  delete new_msg;
+
+  new_msg = 
+    new OctoMsg(m->x_min, m->x_max, m->y_min, m->y_max, z+1, z+1, m->v);
+  Octo(x, y, z+1, new_msg, info);
+}
+
+void Handle22X(int x, int y, int z, OctoMsg* m, char* info)
+{
+  int del_x = m->x_max - m->x_min;
+  int del_y = m->y_max - m->y_min;
+  int del_z = m->z_max - m->z_min;
+
+  // calculate new midpoints
+  int x_lo = (m->x_min + x - 1)/2;  int x_hi = (m->x_max + x)/2;
+  int y_lo = (m->y_min + y - 1)/2;  int y_hi = (m->y_max + y)/2;
+  int z_lo = (m->z_min + z - 1)/2;  int z_hi = (m->z_max + z)/2;
+
+  OctoMsg* new_msg = NULL;
+
+  // 2x2xX
+  if (del_z != 1) {
+    new_msg = 
+      new OctoMsg(m->x_min, m->x_max, m->y_min, m->y_max, z, m->z_max, true);
+    if (z == z_hi) {
+      SplitXYZ(x, y, z_hi, new_msg, info);
+      delete new_msg;
+    }
+    else { Octo(x, y, z_hi, new_msg, info); }
+
+    new_msg = 
+      new OctoMsg(m->x_min, m->x_max, m->y_min, m->y_max, m->z_min, z-1, m->v);
+    if (z_lo == m->z_min && m->v) {
+      SplitXYZ(x, y, z_lo, new_msg, info);
+      delete new_msg;
+    }
+    else { Octo(x, y, z_lo, new_msg, info); }
+  }
+  // 2xXx2
+  else if (del_y != 1) {
+    new_msg = 
+      new OctoMsg(m->x_min, m->x_max, y, m->y_max, m->z_min, m->z_max, true);
+    if (y == y_hi) {
+      SplitXYZ(x, y_hi, z, new_msg, info); 
+      delete new_msg;
+    }
+    else { Octo(x, y_hi, z, new_msg, info); }
+
+    new_msg = 
+      new OctoMsg(m->x_min, m->x_max, m->y_min, y-1, m->z_min, m->z_max, m->v);
+    if (y_lo == m->y_min && m->v) {
+      SplitXYZ(x, y_lo, z, new_msg, info);
+      delete new_msg;
+    }
+    else { Octo(x, y_lo, z, new_msg, info); }
+  }
+  // Xx2x2
+  else { // (del_x != 1)
+    new_msg = 
+      new OctoMsg(x, m->x_max, m->y_min, m->y_max, m->z_min, m->z_max, true);
+    if (x == x_hi) {
+      SplitXYZ(x_hi, y, z, new_msg, info); 
+      delete new_msg;
+    }
+    else { Octo(x_hi, y, z, new_msg, info); }
+
+    new_msg = 
+      new OctoMsg(m->x_min, x-1, m->y_min, m->y_max, m->z_min, m->z_max, m->v);
+    if (x_lo == m->x_min && m->v) {
+      SplitXYZ(x_lo, y, z, new_msg, info);
+      delete new_msg;
+    }
+    else { Octo(x_lo, y, z, new_msg, info); }
+  }
+}
+
+void Handle2XX(int x, int y, int z, OctoMsg* m, char* info)
+{
+  int del_x = m->x_max - m->x_min;
+  int del_y = m->y_max - m->y_min;
+  int del_z = m->z_max - m->z_min;
+
+  // calculate new midpoints
+  int x_lo = (m->x_min + x - 1)/2;  int x_hi = (m->x_max + x)/2;
+  int y_lo = (m->y_min + y - 1)/2;  int y_hi = (m->y_max + y)/2;
+  int z_lo = (m->z_min + z - 1)/2;  int z_hi = (m->z_max + z)/2;
+
+  OctoMsg* new_msg = NULL;
+
+  // get upper corner
+  new_msg = new OctoMsg(x, m->x_max, y, m->y_max, z, m->z_max, true);
+  if (x == x_hi && y == y_hi && z == z_hi) {
+    SplitXYZ(x_hi, y_hi, z_hi, new_msg, info);
+    delete new_msg;
+  }
+  else { Octo(x_hi, y_hi, z_hi, new_msg, info); }
+
+  // 2xXxX 
+  if (del_x == 1) {
+    // lower corner
+    new_msg = new OctoMsg(x, x+1, m->y_min, y-1, m->z_min, z-1, m->v);
+    if (x == m->x_min && y_lo == m->y_min && z_lo == m->z_min && m->v) {
+      SplitXYZ(x, y_lo, z_lo, new_msg, info); 
+      delete new_msg;
+    }
+    else { Octo(x, y_lo, z_lo, new_msg, info); }
+
+    // y_hi, z_lo
+    new_msg = new OctoMsg(x, x+1, y, m->y_max, m->z_min, z-1, false);
+    Octo(x, y_hi, z_lo, new_msg, info);
+
+    // y_lo, z_hi
+    new_msg = new OctoMsg(x, x+1, m->y_min, y-1, z, m->z_max, false);
+    Octo(x, y_lo, z_hi, new_msg, info);
+  }
+  // Xx2xX
+  else if (del_y == 1) {
+    // lower corner
+    new_msg = new OctoMsg(m->x_min, x-1, y, y+1, m->z_min, z-1, m->v);
+    if (x_lo == m->x_min && y == m->y_min && z_lo == m->z_min && m->v) {
+      SplitXYZ(x_lo, y, z_lo, new_msg, info);
+      delete new_msg;
+    }
+    else { Octo(x_lo, y, z_lo, new_msg, info); }
+
+    // x_hi, z_lo
+    new_msg = new OctoMsg(x, m->x_max, y, y+1, m->z_min, z-1, false);
+    Octo(x_hi, y, z_lo, new_msg, info);
+
+    // x_lo, z_hi
+    new_msg = new OctoMsg(m->x_min, x-1, y, y+1, z, m->z_max, false);
+    Octo(x_lo, y, z_hi, new_msg, info);
+  }
+  // XxXx2
+  else if (del_z == 1) {
+    // lower corner
+    new_msg = new OctoMsg(m->x_min, x-1, m->y_min, y-1, z, z+1, m->v);
+    if (x_lo == m->x_min && y_lo == m->y_min && z == m->z_min && m->v) {
+      SplitXYZ(x_lo, y_lo, z, new_msg, info);
+      delete new_msg;
+    }
+    else { Octo(x_lo, y_lo, z, new_msg, info); }
+
+    // x_hi, y_lo
+    new_msg = new OctoMsg(x, m->x_max, m->y_min, y-1, z, z+1, false);
+    Octo(x_hi, y_lo, z, new_msg, info); 
+
+    // x_lo, y_hi
+    new_msg = new OctoMsg(m->x_min, x-1, y, m->y_max, z, z+1, false);
+    Octo(x_lo, y_hi, z, new_msg, info);
+  }
+}
+
diff --git a/examples/bigsim/emulator/ping.C b/examples/bigsim/emulator/ping.C
new file mode 100644 (file)
index 0000000..6e91de8
--- /dev/null
@@ -0,0 +1,145 @@
+#include <stdlib.h>
+#include "blue.h"
+
+int sendID;
+#define MAX_NUM_CASES          4       
+
+extern "C" void sendHandler(char *) ;
+
+const int NUM_ITERATIONS=1<<8;
+
+class Msg
+{
+public:
+char core[CmiBlueGeneMsgHeaderSizeBytes];
+char c[100];
+void *operator new(size_t s) { return CmiAlloc(s); }
+void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+struct userData
+{
+  double st;
+  int    iter;
+
+  int    caseCount;
+  double pingTime[MAX_NUM_CASES] ;
+};
+
+void BgEmulatorInit(int argc, char **argv) 
+{
+  if (argc < 6) { 
+    CmiPrintf("Usage: <program> <x> <y> <z> <numCommTh> <numWorkTh> \n"); 
+    BgShutdown();
+  }
+
+  BgSetSize(atoi(argv[1]), atoi(argv[2]), atoi(argv[3]));
+  BgSetNumCommThread(atoi(argv[4]));
+  BgSetNumWorkThread(atoi(argv[5]));
+
+}
+
+void BgNodeStart(int argc, char **argv) 
+{
+  //ckout << "Initializing node " << bgNode->thisIndex.x << "," 
+  //      << bgNode->thisIndex.y << "," << bgNode->thisIndex.z << endl; 
+  sendID = BgRegisterHandler(sendHandler);
+
+  int x,y,z;
+  BgGetMyXYZ(&x, &y, &z);
+
+  //declare node private data
+  userData *ud = new userData ;
+    ud->st   = 0.;
+    ud->iter = -1;
+
+    ud->caseCount = 1;
+    for(int i=0; i<MAX_NUM_CASES; i++)
+       ud->pingTime[i] = 0.;
+
+  BgSetNodeData((char*)ud) ;
+
+  if(x == 0 && y == 0 && z == 0)
+  {
+       Msg *msg = new Msg;
+       BgSendLocalPacket(ANYTHREAD,sendID, SMALL_WORK, sizeof(Msg), (char *)msg);
+  }
+
+}
+
+void sendHandler(char *info) 
+{
+  int x,y,z;
+  BgGetMyXYZ(&x,&y,&z);
+
+  int numBgX, numBgY, numBgZ;
+  BgGetSize(&numBgX, &numBgY, &numBgZ);
+
+  userData* ud = (userData*)BgGetNodeData();
+  Msg* msg = (Msg*)info;
+
+  if(x==0 && y==0 && z==0)
+  {
+       //ckout <<"Iteration no "<<ud->iter<<endl;
+       ud->iter++;
+       if(ud->iter==0) 
+         ud->st = BgGetTime();
+
+       if(ud->iter==NUM_ITERATIONS)
+       {
+         ud->pingTime[ud->caseCount-1] = (BgGetTime() - ud->st)/NUM_ITERATIONS;
+
+
+         if(ud->caseCount==MAX_NUM_CASES)
+         {
+               CmiPrintf("Pingpong time averaged over %d iterations\n",NUM_ITERATIONS);
+               CmiPrintf("---------------------------------------------------------\n");
+               CmiPrintf("case                 Time(RRT)\n");
+               CmiPrintf("---------------------------------------------------------\n");
+               for(int i=0; i<MAX_NUM_CASES; i++)
+               switch(i+1){
+               case 1:
+                       CmiPrintf("0,0,0 <--> 0,0,%d          %f\n", numBgZ-1,ud->pingTime[0]);
+                       break;
+               case 2:
+                       CmiPrintf("0,0,0 <--> %d,0,0          %f\n",numBgX-1, ud->pingTime[1]);
+                       break;
+               case 3:
+                       CmiPrintf("0,0,0 <--> 0,%d,0          %f\n",numBgY-1, ud->pingTime[2]);
+                       break;
+               case 4:
+                       CmiPrintf("0,0,0 <--> %d,%d,%d          %f\n",numBgX-1,numBgY-1,numBgZ-1,ud->pingTime[3]);
+                       break;
+               }
+               CmiPrintf("---------------------------------------------------------\n");
+
+               BgShutdown();
+               return;
+         }
+         ud->caseCount++;
+         ud->iter = 0;
+         ud->st = BgGetTime();
+
+       }
+
+
+       switch(ud->caseCount) {
+       case 1:
+               BgSendPacket(0,0,numBgZ-1, ANYTHREAD,sendID, SMALL_WORK, sizeof(Msg), (char *)msg);
+               break;
+       case 2:
+               BgSendPacket(numBgX-1,0,0, ANYTHREAD,sendID, SMALL_WORK, sizeof(Msg), (char *)msg);
+               break;
+       case 3:
+               BgSendPacket(0,numBgY-1,0, ANYTHREAD,sendID, SMALL_WORK, sizeof(Msg), (char *)msg);
+               break;
+       case 4:
+               BgSendPacket(numBgX-1,numBgY-1,numBgZ-1, ANYTHREAD,sendID, SMALL_WORK, sizeof(Msg), (char *)msg);
+               break;
+       }
+  }
+  else
+       BgSendPacket(0,0,0, ANYTHREAD,sendID, SMALL_WORK, sizeof(Msg), (char *)msg);
+}
+
+
diff --git a/examples/bigsim/emulator/prime.C b/examples/bigsim/emulator/prime.C
new file mode 100644 (file)
index 0000000..4a8ae48
--- /dev/null
@@ -0,0 +1,248 @@
+#include <stdlib.h>
+#include "blue.h"
+
+int computeNumPrimeID;
+int contributeID;
+int reduceID;
+
+extern "C" void reduce(char *) ;
+extern "C" void computeNumPrime(char *) ;
+extern "C" void contribute(char *) ;
+
+class Msg
+{
+public:
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  int pc;      //number of primes
+  int min, max;        //range of numbers
+  void *operator new(size_t s) { return CmiAlloc(s); }
+  void operator delete(void* ptr) { CmiFree(ptr); }
+};
+
+BnvStaticDeclare(int, pc);
+BnvStaticDeclare(int, count);
+
+void BgEmulatorInit(int argc, char **argv)
+{
+  if (argc < 7) { 
+    CmiPrintf("Usage: <program> <x> <y> <z> <numCommTh> <numWorkTh> <range>\n");
+    BgShutdown();
+  }
+
+  BgSetSize(atoi(argv[1]), atoi(argv[2]), atoi(argv[3]));
+  BgSetNumCommThread(atoi(argv[4]));
+  BgSetNumWorkThread(atoi(argv[5]));
+
+}
+
+void BgNodeStart(int argc, char **argv) 
+{
+  reduceID = BgRegisterHandler(reduce);
+  computeNumPrimeID = BgRegisterHandler(computeNumPrime);
+  contributeID = BgRegisterHandler(contribute);
+
+  int x,y,z;
+  BgGetMyXYZ(&x, &y, &z);
+
+  //CmiPrintf("BgNodeStart in %d %d %d\n", x,y,z);
+
+  //declare node private data
+  BnvInitialize(int, pc);
+  BnvInitialize(int, count);
+  BnvAccess(count) = 0;
+  BnvAccess(pc)    = 0;
+  const int RANGE = atoi(argv[6]); 
+
+  //Decide range of number for this node
+  int min, max;
+  int numBgX, numBgY, numBgZ;
+  BgGetSize(&numBgX, &numBgY, &numBgZ);
+  int slot = RANGE/(numBgX * numBgY * numBgZ);
+  if(0==slot) { slot = 1; }
+
+  min = slot*(x*numBgY*numBgZ + y*numBgZ + z);
+  max = min + slot;
+
+  if(x==numBgX-1 &&
+     y==numBgZ-1 &&
+     z==numBgY-1) 
+       max = RANGE;    
+
+  if(max>RANGE)         { max = RANGE; }
+  if(min>=RANGE) { min=max=0; }
+
+  //ckout << "Range for node " <<x<< ", " <<y<< ", " <<z<< " is " << min <<", "<< max << endl;
+
+  //Divide range to each worker thread and Fire....
+  int numWTh = BgGetNumWorkThread();
+  if(max-min>0)
+  {
+       int numWTh = BgGetNumWorkThread();
+
+       slot = slot/numWTh;
+       if(0==slot) { slot = 1; }
+
+       for(int i=0; i<numWTh; i++)
+       {       
+         Msg *msg = new Msg;
+         msg->min = min + i*slot;
+         msg->max = msg->min + slot;
+
+         if(i==numWTh-1)
+               msg->max = max;
+
+         if(msg->max>max) { msg->max=max; }                    
+         if(msg->min>=max) { msg->min=msg->max=0; }                    
+
+         BgSendLocalPacket(ANYTHREAD,computeNumPrimeID, LARGE_WORK, sizeof(Msg), (char *)msg);
+
+       }
+  }
+  else
+  {
+       for(int i=0; i<numWTh; i++)
+       {
+         Msg *msg = new Msg;
+         msg->min=msg->max=0;
+         BgSendLocalPacket(ANYTHREAD,computeNumPrimeID, LARGE_WORK, sizeof(Msg), (char *)msg);
+       }
+  }
+
+}
+
+
+/* Utility: isPrime
+ * Checks whether a given number is prime or not
+ * Assumption: Any number<=1 is not prime
+ */
+int isPrime(const int number)
+{
+  if(number<=1)        return 0;
+
+  for(int i=2; i<number; i++)
+  {
+       if(0 == number%i)
+         return 0;
+  }
+
+  return 1;
+}
+
+/* Utility: computeNumberOfPrimesIn
+ * Computes number of prime number in a given range
+ */
+int computeNumberOfPrimesIn(const int min, const int max)      
+{
+  int count=0;
+  for(int i=min; i<max; i++)
+  {
+     if(isPrime(i))
+     {
+       count++;
+       //ckout << i << " is prime"<< endl;
+     }
+     else
+       //ckout << i << " is not prime"<< endl;
+       ;
+  }
+
+  return count;
+}
+
+
+/* Handler: computeNumPrime
+ * Compute the number of primes in the assigned range and contribute
+ */
+void computeNumPrime(char *info) 
+{
+  int x,y,z;
+  BgGetMyXYZ(&x,&y,&z);
+  //CmiPrintf("computeNumPrime in %d %d %d\n", x,y,z);
+
+  Msg *msg = (Msg*)info;
+
+  //compute number of primes in the assigned range of numbers
+  int min = msg->min;
+  int max = msg->max;
+  int pc = computeNumberOfPrimesIn(min, max);
+
+  // contribute the results for reduction : reusing the same message for optimization
+  msg->pc = pc;
+  BgSendLocalPacket(ANYTHREAD,contributeID,LARGE_WORK, sizeof(Msg), (char *)msg);
+}
+
+/* Handler: contribute
+ * If the number of contributions received is equal to expected number
+ *    call the reduction handler
+ * Else contribute the contribution to node private data and update counter
+ */
+void contribute(char *info)
+{
+  int x,y,z;
+  BgGetMyXYZ(&x,&y,&z);
+  //ckout << "contribute in " << x << ", " << y << ", " << z << endl ;
+
+  BnvAccess(pc) += ((Msg*)info)->pc;
+  BnvAccess(count)++;
+
+  //compute expected contributions
+  //more: This computation need not be repeated everytime. Change it after this pgm works
+  int reqCount ;
+  int numBgX, numBgY, numBgZ;
+  BgGetSize(&numBgX, &numBgY, &numBgZ);
+  if(z==numBgZ-1)
+    reqCount = 0 ;  
+  else if(z>0 || (z==0 && x==numBgX-1))
+    reqCount = 1 ;
+  else if(x>0 || (x==0 && y==numBgY-1))
+    reqCount = 2 ;  
+  else
+    reqCount = 3 ;
+
+  reqCount += BgGetNumWorkThread();
+
+  if(BnvAccess(count)==reqCount)  //if data for reduction is ready
+  {
+    Msg *msg = (Msg*)info ;
+    msg->pc = BnvAccess(pc);
+    BgSendLocalPacket(ANYTHREAD,reduceID, LARGE_WORK, sizeof(Msg), (char *)msg);
+    return ;
+  }
+  else
+    delete (Msg*)info;
+}
+
+/* Handler: reduce
+ * If reduction has finished, print the number of primes
+ * Else send the number of primes to next node in reduction sequence
+ */
+void reduce(char *info) 
+{
+  int pc = BnvAccess(pc);
+
+  int x,y,z;
+  BgGetMyXYZ(&x,&y,&z);
+  //ckout << "reduce in " << x << ", " << y << ", " << z 
+  //      << " number of primes " << pc << endl;
+
+  if(x==0 && y==0 && z==0)
+  {
+  delete (Msg*)info;
+  CmiPrintf("Finished : The number of primes is %d\n",pc);
+  CmiPrintf("Emulation Time : %f seconds\n", BgGetTime());
+  BgShutdown();
+  return;
+  }
+  //send pc to destination(decide) 
+  if(z>0)
+    z--;   
+  else if(x>0)
+    x--;
+  else
+    y--;
+
+  Msg *msg = (Msg*)info;
+  msg->pc = pc;
+  BgSendPacket(x,y,z, ANYTHREAD,contributeID, LARGE_WORK, sizeof(Msg), (char *)msg);
+}
+
diff --git a/examples/bigsim/emulator/ring.C b/examples/bigsim/emulator/ring.C
new file mode 100644 (file)
index 0000000..26400e1
--- /dev/null
@@ -0,0 +1,62 @@
+#include <stdio.h>
+#include "blue.h"
+
+typedef struct {
+  char core[CmiBlueGeneMsgHeaderSizeBytes];
+  int data;
+} RingMsg;
+
+const int MAXITER = 1;
+int passRingID;
+int iter = 0;
+
+void passRing(char *msg);
+
+void nextxyz(int x, int y, int z, int *nx, int *ny, int *nz)
+{
+  int sx, sy, sz;
+  BgGetSize(&sx, &sy, &sz);
+  *nz = z+1;  *ny = y; *nx = x;
+  if (*nz == sz) {
+    *nz = 0;
+    (*ny) ++;
+    if (*ny == sy) {
+      *ny = 0;
+      (*nx) ++;
+    }
+    if (*nx == sx) {
+      (*nx) = 0;
+    }
+  }
+}
+
+void BgEmulatorInit(int argc, char **argv)
+{
+   if (argc < 6)     CmiAbort("Usage: <ring> <x> <y> <z> <numCommTh> <numWorkTh>\n");
+   BgSetSize(atoi(argv[1]), atoi(argv[2]), atoi(argv[3]));
+   BgSetNumCommThread(atoi(argv[4]));        
+   BgSetNumWorkThread(atoi(argv[5]));
+}
+
+void BgNodeStart(int argc, char **argv)  {
+    int x,y,z, nx, ny, nz;  
+    RingMsg *msg = new RingMsg;                             
+    msg->data =  888;
+    passRingID = BgRegisterHandler(passRing);
+    BgGetMyXYZ(&x, &y, &z);           
+    nextxyz(x, y, z, &nx, &ny, &nz);
+    if (x == 0 && y==0 && z==0) {
+      printf("%d %d %d => %d %d %d\n", x,y,z,nx,ny,nz);
+      BgSendPacket(nx, ny, nz, ANYTHREAD, passRingID, LARGE_WORK, sizeof(RingMsg), (char *)msg);
+    }
+}
+
+void passRing(char *msg)  {
+     int x, y, z,  nx, ny, nz;
+     BgGetMyXYZ(&x, &y, &z);            
+     nextxyz(x, y, z, &nx, &ny, &nz);
+     if (x==0 && y==0 && z==0)     if (++iter == MAXITER) BgShutdown();
+     printf("%d %d %d => %d %d %d\n", x,y,z,nx,ny,nz);
+     BgSendPacket(nx, ny, nz, ANYTHREAD, passRingID, LARGE_WORK, sizeof(RingMsg), msg);
+}
+