/**
# Read X-Y-Z formatted binary data.

This function reades thee-dimensional (cubic $N^3$) data in file `fname` at level `dlev` ($N = 2^{\mathrm{dlev}}$) into the scalar field `s`. 

It works well and fast for large data sets and is compatible with MPI or OpenMP on Multigrid or full octrees  Unfortuantely, it is limited to run on a single node. (use `dump()`->`restore()`)

The first function reads in data stored in so-called single precision
*/
int read_xyz_float(char * fname, scalar s, int dlev){
  unsigned long long int size = (1 << (dimension*dlev));
  /**
The MPI parallel strategy requires special attention. We use the marvalous *shared-memory* functionality that is facilitated by modern MPI implementations.   
  */
#if _MPI 
  MPI_Win win;
  float * a;
  /**
  The root allocated an array and a MPI window is created for the other ranks.
  */
  if (pid() == 0){
    MPI_Win_allocate_shared (size*sizeof(float), sizeof(float), MPI_INFO_NULL,
			    MPI_COMM_WORLD, &a, &win);
  }
  else{ // Slaves obtain the location of the pid()=0 allocated array    
    int disp_unit;
    MPI_Aint  ssize; 
    MPI_Win_allocate_shared (0, sizeof(float), MPI_INFO_NULL,
			    MPI_COMM_WORLD, &a, &win);
    MPI_Win_shared_query (win, 0, &ssize, &disp_unit, &a);
  }
  MPI_Barrier (MPI_COMM_WORLD);
  /**
  The root is also tasked with reading *all* the data. Notice that this is quite fast because it reads contiguous data and true parallel IO is, well... 
  */
  if (pid() == 0){ 
    MPI_Win_lock (MPI_LOCK_EXCLUSIVE, 0, MPI_MODE_NOCHECK, win);
    FILE * fp = fopen (fname, "rb");
    fread (a, sizeof(float), size,fp);
    MPI_Win_unlock (0, win);
  }
  MPI_Barrier (MPI_COMM_WORLD);
  /**
  In serial, life is a bit easier.
  */
#else 
  float * a = (float*) malloc (sizeof(float)*size);
  FILE * fp = fopen (fname, "rb");
  fread (a, sizeof(float), size, fp);
#endif
/** We may need to take case of some specifics of MG parralelism and work with an offset. */
#if MULTIGRID_MPI 
  int nz = pid() % (mpi_dims[2]); 
  int ny = (pid()/mpi_dims[2]) % (mpi_dims[1]);
  int nx = pid()/(mpi_dims[2]*mpi_dims[1]);
  int nxoffset = ((1 << depth())*nx);
  int nyoffset = ((1 << depth())*ny);
  int nzoffset = ((1 << depth())*nz);
#else // Non MG-MPI, no offset 
  int nxoffset = 0;
  int nyoffset = 0;
  int nzoffset = 0;
#endif
  unsigned long long int CR = (1 << dlev);
  int o = -BGHOSTS - 1;
  unsigned long long int index;
/** Loading the data itself is now straightforward*/
  foreach(){
    index = ((nxoffset + point.i + o) + (CR*(nyoffset + point.j + o)) + (sq(CR)*(nzoffset + point.k + o)));
    s[] = (double)a[index];
  }
  return 0;
}
/**
This simular function reads double-precision data:
*/
int read_xyz_double (char * fname,scalar s,int dlev){
  unsigned long long int size= (1<<(dimension*dlev));
#if _MPI
  MPI_Win win;
  double * a;
  if (pid() == 0){
    MPI_Win_allocate_shared(size*sizeof(double), sizeof(double), MPI_INFO_NULL,
			    MPI_COMM_WORLD, &a, &win);
  }
  else{ // Slave    
    int disp_unit;
    MPI_Aint  ssize; 
    MPI_Win_allocate_shared(0, sizeof(double), MPI_INFO_NULL,
			    MPI_COMM_WORLD, &a, &win);
    MPI_Win_shared_query(win, 0, &ssize, &disp_unit, &a);
  }
  MPI_Barrier(MPI_COMM_WORLD);
  if (pid() == 0){
    MPI_Win_lock(MPI_LOCK_EXCLUSIVE, 0, MPI_MODE_NOCHECK, win);
    FILE * fp = fopen(fname,"rb");
    fread(a,sizeof(double),size,fp);
    MPI_Win_unlock(0,win);
  }
  MPI_Barrier(MPI_COMM_WORLD);
#else //not _MPI, life is easier
  double * a= (double*) malloc(sizeof(double)*size);
  FILE * fp = fopen(fname,"rb");
  fread(a,sizeof(double),size,fp);
#endif
  
#if MULTIGRID_MPI
  int nz = pid() % (mpi_dims[2]);
  int ny = (pid()/mpi_dims[2]) % (mpi_dims[1]);
  int nx = pid()/(mpi_dims[2]*mpi_dims[1]);
  unsigned long long int nxoffset = ((1<<depth())*nx);
  unsigned long long int nyoffset = ((1<<depth())*ny);
  unsigned long long int nzoffset = ((1<<depth())*nz);
#else // Non MG-MPI, no offsets 
  int nxoffset = 0;
  int nyoffset = 0;
  int nzoffset = 0;
#endif
  unsigned long long int CR = (1<<dlev);
  int o = -BGHOSTS-1;
  unsigned long long int ind;
  
  foreach(){
    ind =((nxoffset+point.i+o)+(CR*(nyoffset+point.j+o)) + (sq(CR)*(nzoffset+point.k+o)));
    s[]=a[ind];
  }
  return 0;
}