//$ mpicc mpitest.c -o mpitest -std=c99 -lhdf5 -DMPI_DEBUG
//$ NDOMAINS=2 mpirun -np 4 ./mpitest
// http://www.hdfgroup.org/ftp/HDF5/examples/examples-by-api/hdf5-examples/1_8/C/H5D/

#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <assert.h>
#include <unistd.h>
#include <stdio.h>
#include <mpi.h>
#include <hdf5.h>

typedef struct
{
  int rank;
  int size;
  int *members;
  MPI_Comm comm;
  MPI_Group group;
  void *parent;
} MPI_Object;

typedef struct
{
  int start;
  int stop;
  int step;
} slice;

int take(slice s, int idx)
{
  int count = s.stop-s.start;
  assert(0 <= idx && idx < count);
  return s.start + idx*s.step;
}

// ------------------------------------------------------------------

const char datfile[] = "parallelc.hdf5";
const char logfile[] = "parallelc.log";
static int devnull;

#define H5RUN(x,y) if((x=y) < 0) { fprintf(stderr, "H5 %s:%d returned %d.\n", __FILE__, __LINE__, x); return -1; }
#define MPIRUN(x,y) if((x=y) != MPI_SUCCESS) { fprintf(stderr, "MPI %s:%d returned %d.\n", __FILE__, __LINE__, x); return -1; }

#define MPI_WAIT_ORDERED 10000
#define MPI_PRE_ORDERED(comm_obj) \
  for(int i=0; i<(comm_obj).rank; i++) { \
      MPI_Barrier((comm_obj).comm); \
      usleep(MPI_WAIT_ORDERED); \
      }
#define MPI_POST_ORDERED(comm_obj) \
  for(int i=(comm_obj).rank; i<(comm_obj).size; i++) { \
      usleep(MPI_WAIT_ORDERED); \
      MPI_Barrier((comm_obj).comm); \
      }

// ------------------------------------------------------------------

typedef struct
{
  int nbands;
  int mynbands;
  MPI_Object comm;
} BandDescriptor;

BandDescriptor new_band_descriptor(MPI_Object band_comm)
{
  BandDescriptor bd;
  bd.nbands = 24;

  // Determine band descriptors (from band_descriptor.py)
  bd.mynbands = bd.nbands / band_comm.size;
  assert(bd.mynbands * band_comm.size == bd.nbands);
  bd.comm = band_comm;

  return bd;
}

// ------------------------------------------------------------------

typedef struct
{
  int N_c[3]; // global number of grid points per axis incl. padding
  bool pbc_c[3];
  int parsize_c[3];
  int parpos_c[3];
  double h_c[3];
  int beg_c[3];
  int end_c[3];
  int n_c[3]; // local number of grid points per axis
  int gpts_c[3]; // global number of grid points per axis
  MPI_Object comm;
} GridDescriptor;

GridDescriptor new_grid_descriptor(MPI_Object domain_comm)
{
  GridDescriptor gd;
  gd.N_c[0] = gd.N_c[1] = gd.N_c[2] = 12;
  gd.pbc_c[0] = gd.pbc_c[1] = gd.pbc_c[2] = false;

  // Determine domain decomposition (from domain.py)
  gd.parsize_c[0] = domain_comm.size;
  gd.parsize_c[1] = gd.parsize_c[2] = 1;
  assert(gd.parsize_c[0] * gd.parsize_c[1] * gd.parsize_c[2] == domain_comm.size);

  gd.parpos_c[0] = domain_comm.rank / (gd.parsize_c[1] * gd.parsize_c[2]);
  gd.parpos_c[1] = (domain_comm.rank % (gd.parsize_c[1] * gd.parsize_c[2])) / gd.parsize_c[2];
  gd.parpos_c[2] = domain_comm.rank % gd.parsize_c[2];
  assert(gd.parpos_c[0] * gd.parsize_c[1] * gd.parsize_c[2] \
       + gd.parpos_c[1] * gd.parsize_c[2] \
       + gd.parpos_c[2] == domain_comm.rank);

  // Determine grid descriptors (from grid_descriptor.py)
  gd.h_c[0] = gd.h_c[1] = gd.h_c[2] = 0.2;
  for(int c=0; c<3; c++) {
      double tmp;
      int n = gd.parsize_c[c]+1, *p = malloc(n*sizeof(int));
      for(int i=0; i<n; i++) {
          tmp = i * gd.N_c[c];
          tmp /= gd.parsize_c[c];
          tmp += 0.4999;
          if(i==0 && !gd.pbc_c[c])
              p[0] = 1;
          else
              p[i] = tmp;
          //printf("c=%d, i=%d, p=%d\n", c, i, p[i]);
          }
      gd.beg_c[c] = p[gd.parpos_c[c]];
      gd.end_c[c] = p[gd.parpos_c[c] + 1];
      gd.n_c[c] = gd.end_c[c] - gd.beg_c[c];
      gd.gpts_c[c] = gd.N_c[c] - 1 + gd.pbc_c[c];
      free(p);
      }

#ifdef VERBOSE
  // Print ordered grid information
  MPI_Object world = *((MPI_Object*) domain_comm.parent);
  MPI_PRE_ORDERED(world);
  int ndigits = 1 + floor(log10(world.size));
  printf("GRID[%0*d]: gd: %d/%d, pbc=%d,%d,%d; parsize=%d,%d,%d; parpos=%d,%d,%d; beg=%2d,%2d,%2d; end=%2d,%2d,%2d; n=%2d,%2d,%2d\n",
         ndigits, world.rank, domain_comm.rank, domain_comm.size,
         gd.pbc_c[0], gd.pbc_c[1], gd.pbc_c[2],
         gd.parsize_c[0], gd.parsize_c[1], gd.parsize_c[2], 
         gd.parpos_c[0], gd.parpos_c[1], gd.parpos_c[2],
         gd.beg_c[0], gd.beg_c[1], gd.beg_c[2],
         gd.end_c[0], gd.end_c[1], gd.end_c[2],
         gd.n_c[0], gd.n_c[1], gd.n_c[2]);
  fflush(stdout);
  MPI_POST_ORDERED(world);
#endif

  gd.comm = domain_comm;

  return gd;
}

// ------------------------------------------------------------------

int check_fspace(MPI_Object world, hid_t fspace_id, BandDescriptor bd, GridDescriptor gd,
                 int ndims, const hsize_t* refsize)
{
  htri_t valid, simple;
  H5RUN(valid, H5Sselect_valid(fspace_id));
  H5RUN(simple, H5Sis_simple(fspace_id));

  assert((int) H5Sget_simple_extent_ndims(fspace_id) == ndims);
  hsize_t *count = (hsize_t*) malloc(ndims*sizeof(hsize_t));
  H5RUN(devnull, H5Sget_simple_extent_dims(fspace_id, (hsize_t*) count, (hsize_t*) NULL));
  int npts = H5Sget_select_npoints(fspace_id);

#ifdef VERBOSE
  char tbuf[32];
  H5S_sel_type ftype;
  H5RUN(ftype, H5Sget_select_type(fspace_id));

  switch(ftype) {
      case H5S_SEL_ERROR:
          strncpy(tbuf, "ERROR", 32);
          break;
      case H5S_SEL_NONE:
          strncpy(tbuf, "NONE", 32);
          break;
      case H5S_SEL_POINTS:
          strncpy(tbuf, "POINTS", 32);
          break;
      case H5S_SEL_HYPERSLABS:
          strncpy(tbuf, "HYPERSLABS", 32);
          break;
      case H5S_SEL_ALL:
          strncpy(tbuf, "ALL", 32);
          break;
      case H5S_SEL_N:
          strncpy(tbuf, "N", 32);
          break;
      default:
          strncpy(tbuf, "UNKNOWN", 32);
      }

  char cbuf[32];
  switch(ndims) {
      case 3:
          snprintf(cbuf, 32, "{%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2]);
          break;
      case 4:
          snprintf(cbuf, 32, "{%d,%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2], (int) count[3]);
          break;
      default:
          strncpy(cbuf, "UNKNOWN", 32);
      }

  MPI_PRE_ORDERED(world);
  int ndigits = 1 + floor(log10(world.size));
  printf("FILE[%0*d]: valid: %d, simple: %d, npts: %d, H5S: %s, count: %s\n", ndigits, world.rank, valid, simple, npts, tbuf, cbuf);
  MPI_POST_ORDERED(world);
#endif

  assert((int) valid>0);
  assert((int) simple>0);

  int refnpts = 1;
  for(int i=0; i<ndims; i++) {
      assert(count[i] == refsize[i]);
      refnpts *= refsize[i];
      }
  //assert(npts == refnpts); //XXX BAD XXX get selected extent instead!!!

  free(count);
  return 0;
}


// ------------------------------------------------------------------

int check_mspace(MPI_Object world, hid_t mspace_id, BandDescriptor bd, GridDescriptor gd,
                 int ndims, const hsize_t* refsize)
{
  htri_t valid, simple;
  H5RUN(valid, H5Sselect_valid(mspace_id));
  H5RUN(simple, H5Sis_simple(mspace_id));

  assert((int) H5Sget_simple_extent_ndims(mspace_id) == ndims);
  hsize_t *count = (hsize_t*) malloc(ndims*sizeof(hsize_t));
  H5RUN(devnull, H5Sget_simple_extent_dims(mspace_id, (hsize_t*) count, (hsize_t*) NULL));
  int npts = H5Sget_select_npoints(mspace_id);

#ifdef VERBOSE
  char tbuf[32];
  H5S_sel_type mtype;
  H5RUN(mtype, H5Sget_select_type(mspace_id));

  switch(mtype) {
      case H5S_SEL_ERROR:
          strncpy(tbuf, "ERROR", 32);
          break;
      case H5S_SEL_NONE:
          strncpy(tbuf, "NONE", 32);
          break;
      case H5S_SEL_POINTS:
          strncpy(tbuf, "POINTS", 32);
          break;
      case H5S_SEL_HYPERSLABS:
          strncpy(tbuf, "HYPERSLABS", 32);
          break;
      case H5S_SEL_ALL:
          strncpy(tbuf, "ALL", 32);
          break;
      case H5S_SEL_N:
          strncpy(tbuf, "N", 32);
          break;
      default:
          strncpy(tbuf, "UNKNOWN", 32);
      }

  char cbuf[32];
  switch(ndims) {
      case 3:
          snprintf(cbuf, 32, "{%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2]);
          break;
      case 4:
          snprintf(cbuf, 32, "{%d,%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2], (int) count[3]);
          break;
      default:
          strncpy(cbuf, "UNKNOWN", 32);
      }

  MPI_PRE_ORDERED(world);
  int ndigits = 1 + floor(log10(world.size));
  printf("MEMO[%0*d]: valid: %d, simple: %d, npts: %d, H5S: %s, count: %s\n", ndigits, world.rank, valid, simple, npts, tbuf, cbuf);
  fflush(stdout);
  MPI_POST_ORDERED(world);
#endif

  assert((int) valid>0);
  assert((int) simple>0);

  int refnpts = 1;
  for(int i=0; i<ndims; i++) {
      assert(count[i] == refsize[i]);
      refnpts *= refsize[i];
      }
  //assert(npts == refnpts); //XXX BAD XXX can be zero if none selected!!!

  free(count);
  return 0;
}

// ------------------------------------------------------------------

int check_dset(MPI_Object world, hid_t dset_id, BandDescriptor bd, GridDescriptor gd,
               int ndims, const hsize_t* count, int typesize)
{
  int nbytes = H5Dget_storage_size(dset_id), refnbytes = typesize;
  for(int i=0; i<ndims; i++) {
      refnbytes *= count[i];
      }

#ifdef VERBOSE
  hid_t dcpl_id;
  H5RUN(dcpl_id, H5Dget_create_plist(dset_id));

  char abuf[32];
  H5D_alloc_time_t alloctime;
  H5RUN(devnull, H5Pget_alloc_time(dcpl_id, &alloctime));

  switch(alloctime) {
      case H5D_ALLOC_TIME_ERROR:
          strncpy(abuf, "ERROR", 32);
          break;
      case H5D_ALLOC_TIME_DEFAULT:
          strncpy(abuf, "DEFAULT", 32);
          break;
      case H5D_ALLOC_TIME_EARLY:
          strncpy(abuf, "EARLY", 32);
          break;
      case H5D_ALLOC_TIME_LATE:
          strncpy(abuf, "LATE", 32);
          break;
      case H5D_ALLOC_TIME_INCR:
          strncpy(abuf, "INCR", 32);
          break;
      default:
          strncpy(abuf, "UNKNOWN", 32);
      }

  char fbuf[32];
  H5D_fill_time_t filltime;
  H5RUN(devnull, H5Pget_fill_time(dcpl_id, &filltime));

  switch(filltime) {
      case H5D_FILL_TIME_ERROR:
          strncpy(fbuf, "ERROR", 32);
          break;
      case H5D_FILL_TIME_ALLOC:
          strncpy(fbuf, "ALLOC", 32);
          break;
      case H5D_FILL_TIME_NEVER:
          strncpy(fbuf, "NEVER", 32);
          break;
      case H5D_FILL_TIME_IFSET:
          strncpy(fbuf, "IFSET", 32);
          break;
      default:
          strncpy(fbuf, "UNKNOWN", 32);
      }

  char lbuf[32];
  H5D_layout_t layout;
  H5RUN(layout, H5Pget_layout(dcpl_id));

  switch(layout) {
      case H5D_COMPACT:
          strncpy(lbuf, "COMPACT", 32);
          break;
      case H5D_CONTIGUOUS:
          strncpy(lbuf, "CONTIGUOUS", 32);
          break;
      case H5D_CHUNKED:
          strncpy(lbuf, "CHUNKED", 32);
          break;
      default:
          strncpy(lbuf, "UNKNOWN", 32);
      }

  H5Pclose(dcpl_id);

  char sbuf[32];
  H5D_space_status_t stat;
  H5RUN(devnull, H5Dget_space_status(dset_id, &stat));

  switch(stat) {
      case H5D_SPACE_STATUS_ERROR:
          strncpy(sbuf, "ERROR", 32);
          break;
      case H5D_SPACE_STATUS_NOT_ALLOCATED:
          strncpy(sbuf, "NOT_ALLOCATED", 32);
          break;
      case H5D_SPACE_STATUS_PART_ALLOCATED:
          strncpy(sbuf, "PART_ALLOCATED", 32);
          break;
      case H5D_SPACE_STATUS_ALLOCATED:
          strncpy(sbuf, "ALLOCATED", 32);
          break;
      default:
          strncpy(sbuf, "UNKNOWN", 32);
      }

  char cbuf[32];
  switch(ndims) {
      case 3:
          snprintf(cbuf, 32, "{%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2]);
          break;
      case 4:
          snprintf(cbuf, 32, "{%d,%d,%d,%d}", (int) count[0], (int) count[1], (int) count[2], (int) count[3]);
          break;
      default:
          strncpy(cbuf, "UNKNOWN", 32);
      }

  MPI_PRE_ORDERED(world);
  int ndigits = 1 + floor(log10(world.size));
  printf("DSET[%0*d]: H5D: A=%s, F=%s, L=%s, S=%s, size: %d (%s*%d=%d)\n", ndigits, world.rank, abuf, fbuf, lbuf, sbuf, nbytes, cbuf, typesize, refnbytes);
  MPI_POST_ORDERED(world);
#endif

  assert(nbytes == refnbytes);

  return 0;
}

// ------------------------------------------------------------------

int run_write(MPI_Object world, MPI_Object domain_comm, MPI_Object band_comm)
{
  // Print ordered processor distribution
  MPI_PRE_ORDERED(world);
  printf("WRITE: %d/%d, gd: %d/%d, bd: %d/%d\n", world.rank, world.size, domain_comm.rank, domain_comm.size, band_comm.rank, band_comm.size);
  fflush(stdout);
  MPI_POST_ORDERED(world);

  GridDescriptor gd = new_grid_descriptor(domain_comm);
  BandDescriptor bd = new_band_descriptor(band_comm);

  // Set up slices
  slice *indices = malloc(4*sizeof(slice));
  indices[0].start = band_comm.rank * bd.mynbands;
  indices[0].stop = indices[0].start + bd.mynbands;
  indices[0].step = 1;
  for(int c=0; c<3; c++) {
      indices[c+1].start = gd.beg_c[c] - 1 + gd.pbc_c[c];
      indices[c+1].stop = gd.end_c[c] - 1 + gd.pbc_c[c];
      indices[c+1].step = 1;
      }

  // Fill in bogus but reproducible data in local domain
  int myngpts = gd.n_c[0] * gd.n_c[1] * gd.n_c[2];
  double *psit_nG = malloc(bd.mynbands * myngpts * sizeof(double));
  if(psit_nG) {
      double tmp, tmp2, r0_v[3];
      for(int myn=0; myn<bd.mynbands; myn++) {
          int n = take(indices[0], myn), myg = 0;
          for(int c=0; c<3; c++)
              r0_v[c] = gd.pbc_c[c]?n*gd.h_c[c]*gd.N_c[c]:0;
          for(int gx=gd.beg_c[0]; gx<gd.end_c[0]; gx++) {
              tmp = pow(gd.h_c[0] * gx - r0_v[0], 2);
              for(int gy=gd.beg_c[1]; gy<gd.end_c[1]; gy++) {
                  tmp2 = tmp + pow(gd.h_c[1] * gy - r0_v[1], 2);
                  for(int gz=gd.beg_c[2]; gz<gd.end_c[2]; gz++) {
                      psit_nG[myn*myngpts+myg] = sqrt(tmp2 + pow(gd.h_c[2] * gz - r0_v[2], 2));
                      myg++;
                      }
                  }
              }
          assert(myg==myngpts);
          }
      }

  // Fill in bogus but reproducible data on domain masters
  int ngpts = gd.gpts_c[0] * gd.gpts_c[1] * gd.gpts_c[2];
  double *big_rhot_G = malloc(ngpts * sizeof(double));
  if(big_rhot_G && domain_comm.rank == 0) {
      double tmp, tmp2, r0_v[3];
      int g = 0;
      //XXX lazy global indexing should be offset by 1-pbc_c
      for(int c=0; c<3; c++)
          r0_v[c] = 0.5*gd.h_c[c]*gd.N_c[c];
      for(int gx=0; gx<gd.gpts_c[0]; gx++) {
          tmp = exp(-pow(gd.h_c[0] * gx - r0_v[0], 2));
          for(int gy=0; gy<gd.gpts_c[1]; gy++) {
              tmp2 = tmp * exp(-pow(gd.h_c[1] * gy - r0_v[1], 2));
              for(int gz=0; gz<gd.gpts_c[2]; gz++) {
                  big_rhot_G[g] = tmp2 * exp(-pow(gd.h_c[2] * gz - r0_v[2], 2));
                  g++;
                  }
              }
          }
      assert(g==ngpts);
      }      

  hid_t fapl_id = 0, f_id = 0, dxpl_id = 0, dcpl_id = 0;

  // File access property list (from io/hdf5.py)
  MPI_Info info;
  H5RUN(fapl_id, H5Pcreate(H5P_FILE_ACCESS));
#ifndef MPI_DEBUG
  H5RUN(devnull, H5Pset_fclose_degree(fapl_id, H5F_CLOSE_STRONG));
#endif
  if(world.size > 1) {
      // Copy MPI information to the file access property list (driver)
      char buf[20];
      snprintf(buf, 20, "%d", world.size);
      MPIRUN(devnull, MPI_Info_create(&info));
      MPIRUN(devnull, MPI_Info_set(info, "cb_nodes", buf));
      //H5RUN(devnull, H5FD_mpio_init());
      H5RUN(devnull, H5Pset_fapl_mpio(fapl_id, world.comm, info));
      }
  else {
      H5RUN(devnull, H5Pset_fapl_stdio(fapl_id));
      //H5RUN(devnull, H5Pset_fapl_log(fapl_id, logfile, H5FD_LOG_ALL, 256*1024));
      }

  // Create file in write-mode using the file access property list from above
  H5RUN(f_id, H5Fcreate(datfile, H5F_ACC_TRUNC /*'w'*/, H5P_DEFAULT /*fcpl*/, fapl_id /*fapl*/));

  // Speficy collective I/O for data transfer property list
  H5RUN(dxpl_id, H5Pcreate(H5P_DATASET_XFER));
  if(world.size > 1) {
      H5RUN(devnull, H5Pset_dxpl_mpio(dxpl_id, H5FD_MPIO_COLLECTIVE));
      //H5RUN(devnull, H5Pset_dxpl_mpio_chunk_opt(dxpl_id, H5FD_MPIO_CHUNK_MULTI_IO)); //new in 1.8x?
      //H5RUN(devnull, H5Pset_buffer(dxpl_id, 32*1024*1024)); //set transfer buffer to 32mb
      }

  // Specify early allocation for dataset creation property list if serial
  H5RUN(dcpl_id, H5Pcreate(H5P_DATASET_CREATE));
  H5RUN(devnull, H5Pset_fill_time(dcpl_id, H5D_FILL_TIME_NEVER));
  if(world.size == 1) {
      H5RUN(devnull, H5Pset_alloc_time(dcpl_id, H5D_ALLOC_TIME_EARLY));
      }

  hid_t dset_id = 0, dtype_id = 0, fspace_id = 0, mspace_id = 0;

  if(psit_nG) {
      if(world.rank == 0) {
          printf("\npsit_nG\n");
          fflush(stdout);
          }

      hsize_t count[] = {bd.mynbands, gd.n_c[0], gd.n_c[1], gd.n_c[2]},
              start[] = {indices[0].start, indices[1].start, indices[2].start, indices[3].start},
              fshape[] = {bd.nbands, gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]};

      // Create reduced-size memory dataspace
      H5RUN(mspace_id, H5Screate_simple(4, (const hsize_t*) count, NULL));
      check_mspace(world, mspace_id, bd, gd, 4, (const hsize_t*) count);

      // Create fixed-size file dataspace
      H5RUN(fspace_id, H5Screate_simple(4, (const hsize_t*) fshape, (const hsize_t*) NULL)); // max_shape???

      // Select the hyperslab that memory dataspace represents in file dataspace
      H5RUN(devnull, H5Sselect_hyperslab(fspace_id, H5S_SELECT_SET /*op*/, (const hsize_t*) start, (const hsize_t*) NULL /*stride*/, (const hsize_t*) count, (const hsize_t*) NULL /*block*/));
      //H5RUN(devnull, H5Sselect_none(fspace_id));
      //H5RUN(devnull, H5Soffset_simple(fspace_id, (const hsize_t*) zeros));
      check_fspace(world, fspace_id, bd, gd, 4, (const hsize_t*) fshape);

      // Create datatype and verify size and order
      H5RUN(dtype_id, H5Tcopy(H5T_NATIVE_DOUBLE));
      H5RUN(devnull, H5Tset_order(dtype_id, H5T_ORDER_LE));
      assert(H5Tget_size(dtype_id) == sizeof(double));

      // Create psit_nG dataset
      H5RUN(dset_id, H5Dcreate(f_id, "PseudoWaveFunctions", dtype_id /*ftype*/, fspace_id /*fspace*/, dcpl_id /*dcpl*/));

      check_dset(world, dset_id, bd, gd, 4, (const hsize_t*) fshape, sizeof(double));

      H5RUN(devnull, H5Dwrite(dset_id, dtype_id /*mtype*/, mspace_id /*mspace*/, fspace_id /*fspace*/, dxpl_id /*dxpl*/, psit_nG));

      // Close all open ids
      H5RUN(devnull, H5Dclose(dset_id));
      H5RUN(devnull, H5Tclose(dtype_id));
      H5RUN(devnull, H5Sclose(fspace_id));
      H5RUN(devnull, H5Sclose(mspace_id));
      }

  if(big_rhot_G) {
      if(world.rank == 0) {
          printf("\nbig_rhot_G\n");
          fflush(stdout);
          }

      hsize_t count[] = {gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]},
              start[] = {0, 0, 0},
              fshape[] = {gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]};

      // Create reduced-size memory dataspace
      H5RUN(mspace_id, H5Screate_simple(3, (const hsize_t*) count, NULL));
      if(domain_comm.rank != 0) {
          H5RUN(devnull, H5Sselect_none(mspace_id));
          }
      check_mspace(world, mspace_id, bd, gd, 3, (const hsize_t*) count);

      // Create fixed-size file dataspace
      H5RUN(fspace_id, H5Screate_simple(3, (const hsize_t*) fshape, (const hsize_t*) NULL)); // max_shape???

      // Select the hyperslab that memory dataspace represents in file dataspace
      if(domain_comm.rank == 0) {
          H5RUN(devnull, H5Sselect_hyperslab(fspace_id, H5S_SELECT_SET /*op*/, (const hsize_t*) start, (const hsize_t*) NULL /*stride*/, (const hsize_t*) count, (const hsize_t*) NULL /*block*/));
          }
      else {
          H5RUN(devnull, H5Sselect_none(fspace_id));
          }
      //H5RUN(devnull, H5Soffset_simple(fspace_id, (const hsize_t*) zeros));
      check_fspace(world, fspace_id, bd, gd, 3, (const hsize_t*) fshape);

      // Create datatype and verify size and order
      H5RUN(dtype_id, H5Tcopy(H5T_NATIVE_DOUBLE));
      H5RUN(devnull, H5Tset_order(dtype_id, H5T_ORDER_LE));
      assert(H5Tget_size(dtype_id) == sizeof(double));

      // Create psit_nG dataset
      H5RUN(dset_id, H5Dcreate(f_id, "PseudoChargeDensity", dtype_id /*ftype*/, fspace_id /*fspace*/, dcpl_id /*dcpl*/));
      check_dset(world, dset_id, bd, gd, 3, (const hsize_t*) fshape, sizeof(double));

      H5RUN(devnull, H5Dwrite(dset_id, dtype_id /*mtype*/, mspace_id /*mspace*/, fspace_id /*fspace*/, dxpl_id /*dxpl*/, big_rhot_G));

      // Close all open ids
      H5RUN(devnull, H5Dclose(dset_id));
      H5RUN(devnull, H5Tclose(dtype_id));
      H5RUN(devnull, H5Sclose(fspace_id));
      H5RUN(devnull, H5Sclose(mspace_id));
      }

  // Close remaining ids
  H5RUN(devnull, H5Pclose(dcpl_id));
  H5RUN(devnull, H5Pclose(dxpl_id));
  H5RUN(devnull, H5Pclose(fapl_id));
  if(world.size > 1)
      MPIRUN(devnull, MPI_Info_free(&info));
  H5RUN(devnull, H5Fclose(f_id));

  free(psit_nG);
  free(big_rhot_G);
  fflush(stdout);
  MPI_Barrier(world.comm);

  return 0;
}

// ------------------------------------------------------------------

int run_read(MPI_Object world, MPI_Object domain_comm, MPI_Object band_comm)
{
  // Print ordered processor distribution
  MPI_PRE_ORDERED(world);
  printf("READ: %d/%d, gd: %d/%d, bd: %d/%d\n", world.rank, world.size, domain_comm.rank, domain_comm.size, band_comm.rank, band_comm.size);
  fflush(stdout);
  MPI_POST_ORDERED(world);

  GridDescriptor gd = new_grid_descriptor(domain_comm);
  BandDescriptor bd = new_band_descriptor(band_comm);

  // Set up slices
  slice *indices = malloc(4*sizeof(slice));
  indices[0].start = band_comm.rank * bd.mynbands;
  indices[0].stop = indices[0].start + bd.mynbands;
  indices[0].step = 1;
  for(int c=0; c<3; c++) {
      indices[c+1].start = gd.beg_c[c] - 1 + gd.pbc_c[c];
      indices[c+1].stop = gd.end_c[c] - 1 + gd.pbc_c[c];
      indices[c+1].step = 1;
      }

  // Allocate buffer for bogus but reproducible data in local domain
  int myngpts = gd.n_c[0] * gd.n_c[1] * gd.n_c[2];
  double *psit_nG = malloc(bd.mynbands * myngpts * sizeof(double));
  if(!psit_nG)
      return -1;

  // Allocate buffer for bogus but reproducible data on domain masters
  int ngpts = gd.gpts_c[0] * gd.gpts_c[1] * gd.gpts_c[2];
  double *big_rhot_G = malloc(ngpts * sizeof(double));
  if(!big_rhot_G)
      return -1;

  hid_t fapl_id = 0, f_id = 0, dxpl_id = 0;

  // File access property list (from io/hdf5.py)
  MPI_Info info;
  H5RUN(fapl_id, H5Pcreate(H5P_FILE_ACCESS));
#ifndef MPI_DEBUG
  H5RUN(devnull, H5Pset_fclose_degree(fapl_id, H5F_CLOSE_STRONG));
#endif
  if(world.size > 1) {
      // Copy MPI information to the file access property list (driver)
      char buf[20];
      snprintf(buf, 20, "%d", world.size);
      MPIRUN(devnull, MPI_Info_create(&info));
      MPIRUN(devnull, MPI_Info_set(info, "cb_nodes", buf));
      //H5RUN(devnull, H5FD_mpio_init());
      H5RUN(devnull, H5Pset_fapl_mpio(fapl_id, world.comm, info));
      }
  else
      H5RUN(devnull, H5Pset_fapl_stdio(fapl_id));
      //H5RUN(devnull, H5Pset_fapl_log(fapl_id, logfile, H5FD_LOG_ALL, 256*1024));

  // Create file in write-mode using the file access property list from above
  H5RUN(f_id, H5Fopen(datfile, H5F_ACC_RDONLY /*'r'*/, fapl_id /*fapl*/));

  // Speficy collective I/O for data transfer property list
  H5RUN(dxpl_id, H5Pcreate(H5P_DATASET_XFER));
  if(world.size > 1) {
      H5RUN(devnull, H5Pset_dxpl_mpio(dxpl_id, H5FD_MPIO_COLLECTIVE));
      //H5RUN(devnull, H5Pset_dxpl_mpio_chunk_opt(dxpl_id, H5FD_MPIO_CHUNK_MULTI_IO)); //new in 1.8x?
      //H5RUN(devnull, H5Pset_buffer(dxpl_id, 32*1024*1024)); //set transfer buffer to 32mb
      }

  hid_t dset_id = 0, dtype_id = 0, fspace_id = 0, mspace_id = 0;

  if(psit_nG) {
      if(world.rank == 0) {
          printf("\npsit_nG\n");
          fflush(stdout);
          }
      hsize_t refsize[] = {bd.nbands, gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]},
              count[] = {bd.mynbands, gd.n_c[0], gd.n_c[1], gd.n_c[2]},
              start[] = {indices[0].start, indices[1].start, indices[2].start, indices[3].start};

      // Open psit_nG dataset
      H5RUN(dset_id, H5Dopen(f_id, "PseudoWaveFunctions"));
      check_dset(world, dset_id, bd, gd, 4, (const hsize_t*) refsize, sizeof(double));

      // Open datatype and verify size and order
      H5RUN(dtype_id, H5Dget_type(dset_id));
      assert((int) H5Tequal(dtype_id, H5T_NATIVE_DOUBLE) > 0);
      assert(H5Tget_size(dtype_id) == sizeof(double));
      assert(H5Tget_order(dtype_id) == H5T_ORDER_LE);

      // Open file dataspace and verify fixed size
      H5RUN(fspace_id, H5Dget_space(dset_id));
      assert(H5Sget_simple_extent_ndims(fspace_id) == 4);
      hsize_t fshape[4];

      // Select the hyperslab that memory dataspace represents in file dataspace
      H5RUN(devnull, H5Sselect_hyperslab(fspace_id, H5S_SELECT_SET /*op*/, (const hsize_t*) start, (const hsize_t*) NULL /*stride*/, (const hsize_t*) count, (const hsize_t*) NULL /*block*/));
      //H5RUN(devnull, H5Sselect_none(fspace_id));
      //H5RUN(devnull, H5Soffset_simple(fspace_id, (const hsize_t*) zeros));
      check_fspace(world, fspace_id, bd, gd, 4, (const hsize_t*) refsize);

      // Create reduced-size memory dataspace
      H5RUN(mspace_id, H5Screate_simple(4, (const hsize_t*) count, NULL));
      check_mspace(world, mspace_id, bd, gd, 4, (const hsize_t*) count);

      H5RUN(devnull, H5Dread(dset_id, dtype_id /*mtype*/, mspace_id /*mspace*/, fspace_id /*fspace*/, dxpl_id /*dxpl*/, psit_nG));

      // Close all open ids
      H5RUN(devnull, H5Dclose(dset_id));
      H5RUN(devnull, H5Tclose(dtype_id));
      H5RUN(devnull, H5Sclose(fspace_id));
      H5RUN(devnull, H5Sclose(mspace_id));
      }

  if(big_rhot_G) {
      if(world.rank == 0) {
          printf("\nbig_rhot_G\n");
          fflush(stdout);
          }
      hsize_t refsize[] = {gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]},
              count[] = {gd.gpts_c[0], gd.gpts_c[1], gd.gpts_c[2]},
              start[] = {0, 0, 0};

      // Open psit_nG dataset
      H5RUN(dset_id, H5Dopen(f_id, "PseudoChargeDensity"));
      check_dset(world, dset_id, bd, gd, 3, (const hsize_t*) refsize, sizeof(double));

      // Open datatype and verify size and order
      H5RUN(dtype_id, H5Dget_type(dset_id));
      assert((int) H5Tequal(dtype_id, H5T_NATIVE_DOUBLE) > 0);
      assert(H5Tget_size(dtype_id) == sizeof(double));
      assert(H5Tget_order(dtype_id) == H5T_ORDER_LE);

      // Open file dataspace and verify fixed size
      H5RUN(fspace_id, H5Dget_space(dset_id));
      assert(H5Sget_simple_extent_ndims(fspace_id) == 3);
      hsize_t fshape[3];

      // Select the hyperslab that memory dataspace represents in file dataspace
      if(domain_comm.rank == 0) {
          H5RUN(devnull, H5Sselect_hyperslab(fspace_id, H5S_SELECT_SET /*op*/, (const hsize_t*) start, (const hsize_t*) NULL /*stride*/, (const hsize_t*) count, (const hsize_t*) NULL /*block*/));
          }
      else {
          H5RUN(devnull, H5Sselect_none(fspace_id));
          }
      //H5RUN(devnull, H5Soffset_simple(fspace_id, (const hsize_t*) zeros));
      check_fspace(world, fspace_id, bd, gd, 3, (const hsize_t*) refsize);

      // Create reduced-size memory dataspace
      H5RUN(mspace_id, H5Screate_simple(3, (const hsize_t*) count, NULL));
      if(domain_comm.rank != 0) {
          H5RUN(devnull, H5Sselect_none(mspace_id));
          }
      check_mspace(world, mspace_id, bd, gd, 3, (const hsize_t*) count);

      H5RUN(devnull, H5Dread(dset_id, dtype_id /*mtype*/, mspace_id /*mspace*/, fspace_id /*fspace*/, dxpl_id /*dxpl*/, big_rhot_G));

      // Close all open ids
      H5RUN(devnull, H5Dclose(dset_id));
      H5RUN(devnull, H5Tclose(dtype_id));
      H5RUN(devnull, H5Sclose(fspace_id));
      H5RUN(devnull, H5Sclose(mspace_id));
      }

  H5RUN(devnull, H5Pclose(dxpl_id));
  H5RUN(devnull, H5Pclose(fapl_id));
  if(world.size > 1)
      MPIRUN(devnull, MPI_Info_free(&info));
  H5RUN(devnull, H5Fclose(f_id));

  if(psit_nG) {
      // Verify read data against bogus but reproducible equivalent
      double tmp, tmp2, r0_v[3];
      for(int myn=0; myn<bd.mynbands; myn++) {
          int n = take(indices[0], myn), myg = 0;
          for(int c=0; c<3; c++)
              r0_v[c] = gd.pbc_c[c]?n*gd.h_c[c]*gd.N_c[c]:0;
          for(int gx=gd.beg_c[0]; gx<gd.end_c[0]; gx++) {
              tmp = pow(gd.h_c[0] * gx - r0_v[0], 2);
              for(int gy=gd.beg_c[1]; gy<gd.end_c[1]; gy++) {
                  tmp2 = tmp + pow(gd.h_c[1] * gy - r0_v[1], 2);
                  for(int gz=gd.beg_c[2]; gz<gd.end_c[2]; gz++) {
                      assert(abs(psit_nG[myn*myngpts+myg] - sqrt(tmp2 + pow(gd.h_c[2] * gz - r0_v[2], 2))) < 1e-9);
                      myg++;
                      }
                  }
              }
          assert(myg==myngpts);
          }
      free(psit_nG);
      }

  if(big_rhot_G) {
      MPIRUN(devnull, MPI_Bcast(big_rhot_G, ngpts, MPI_DOUBLE, 0, domain_comm.comm));

      // Verify read data against bogus but reproducible equivalent
      double tmp, tmp2, r0_v[3];
      int g = 0;
      //XXX lazy global indexing should be offset by 1-pbc_c
      for(int c=0; c<3; c++)
          r0_v[c] = 0.5*gd.h_c[c]*gd.N_c[c];
      for(int gx=0; gx<gd.gpts_c[0]; gx++) {
          tmp = exp(-pow(gd.h_c[0] * gx - r0_v[0], 2));
          for(int gy=0; gy<gd.gpts_c[1]; gy++) {
              tmp2 = tmp * exp(-pow(gd.h_c[1] * gy - r0_v[1], 2));
              for(int gz=0; gz<gd.gpts_c[2]; gz++) {
                  assert(abs(big_rhot_G[g] - tmp2 * exp(-pow(gd.h_c[2] * gz - r0_v[2], 2))) < 1e-9);
                  g++;
                  }
              }
          }
      assert(g==ngpts);
      free(big_rhot_G);
      }

  fflush(stdout);
  MPI_Barrier(world.comm);

  return 0;
}

// ------------------------------------------------------------------

int main(int argc, char **argv)
{
  int status;
  MPI_Init(&argc, &argv);

#ifdef MPI_DEBUG
  // Default Errhandler is MPI_ERRORS_ARE_FATAL
  MPI_Errhandler_set(MPI_COMM_WORLD, MPI_ERRORS_RETURN);
#endif

  MPI_Barrier(MPI_COMM_WORLD);

  // Basic info on world communicator
  int rank, nprocs, root;
  MPI_Object world;
  world.comm = MPI_COMM_WORLD;
  MPI_Comm_group(world.comm, &(world.group));
  MPI_Comm_rank(world.comm, &(world.rank));
  MPI_Comm_size(world.comm, &(world.size));
  world.members = NULL; //XXX not needed, right?
  world.parent = NULL; //XXX not needed, right?

  // Determine communicator sizes
  rank = world.rank;
  nprocs = world.size;
  int ndomains = world.size, ngroups = 1;
  char *p = getenv("NDOMAINS");
  if(p) {
      ndomains = atoi(p);
      if(ndomains <= 0 || nprocs < ndomains || nprocs % ndomains != 0) {
          if(rank == 0)
              fprintf(stderr, "Invalid number of domains (%d).\n", ndomains);
          MPI_Finalize();
          return -1;
          }
      ngroups = nprocs / ndomains;
      }
  assert(ndomains * ngroups == nprocs);

  // Create domain communicator
  MPI_Object domain_comm;
  domain_comm.members = (int*) malloc(ndomains*sizeof(int));
  root = rank - rank % ndomains;
  for(int i=0; i<ndomains; i++)
      domain_comm.members[i] = root + i;
  MPI_Group_incl(world.group, ndomains, domain_comm.members, &(domain_comm.group));
  MPI_Comm_create(world.comm, domain_comm.group, &(domain_comm.comm)); // has a memory leak!
  MPI_Comm_rank(domain_comm.comm, &(domain_comm.rank));
  MPI_Comm_size(domain_comm.comm, &(domain_comm.size));
  domain_comm.parent = (void*) &world;

  // Create band communicator
  MPI_Object band_comm;
  band_comm.members = (int*) malloc(ngroups*sizeof(int));
  root = rank % ndomains;
  for(int i=0; i<ngroups; i++)
      band_comm.members[i] = root + i * ndomains;
  MPI_Group_incl(world.group, ngroups, band_comm.members, &(band_comm.group));
  MPI_Comm_create(world.comm, band_comm.group, &(band_comm.comm)); // has a memory leak!
  MPI_Comm_rank(band_comm.comm, &(band_comm.rank));
  MPI_Comm_size(band_comm.comm, &(band_comm.size));
  band_comm.parent = (void*) &world;

  bool write = true, read = true;
  p = getenv("ACTION");
  if(p) {
      assert(strlen(p) == 1);
      switch(p[0]) {
          case 'w':
          case 'W':
              write = true;
              read = false;
              break;
          case 'r':
          case 'R':
              write = false;
              read = true;
              break;
          case 'b':
          case 'B':
              write = true;
              read = true;
              break;
          default:
              if(world.rank == 0) {
                  fprintf(stderr, "\nWARNING: Action '%c' not understood.\n\n", p[0]);
                  fflush(stderr);
                  }
          }
      }

  if(write)
      status = run_write(world, domain_comm, band_comm);

  if(write && read && world.rank == 0) {
      printf("-----------------------------------------------\n");
      fflush(stdout);
    }

  if(read)
      status = run_read(world, domain_comm, band_comm);

  free(domain_comm.members);
  free(band_comm.members);

  MPI_Finalize();
  return status;
}