/*@@
   @file      DumpVar.c
   @date      01 Oct 1999
   @author    Jonghyun Lee
   @desc      Do the actual writing of a 3D grid function,
              for output or for checkpointing
   @enddesc 
   @history
   @hendhistory
 @@*/

#include <stdio.h>
#include <stdlib.h>
#ifdef SGI
#include <time.h>
#endif

#include "cctk.h"
#include "cctk_DefineThorn.h"
#include "cctk_Parameters.h"
#ifdef CACTUSPUGH_PUGH
#include "CactusPUGH/PUGH/src/include/pugh.h"
#endif
#include "CactusBase/IOUtil/src/ioGH.h"
#include "ioPandaGH.h"


#define IOTAGBASE 20000 /* This may break on more than 2000 processors */


static char *char_time_date = NULL;


void IOPanda_getDumpData (cGH *GH, int index, int timelevel, void **outme,
                         int *free_outme, CCTK_INT4 bnd [9], int element_size) 
{
  DECLARE_CCTK_PARAMETERS
  int i;
  int myproc;
  ioGH *ioUtilGH;
  pGH *pughGH;
  CCTK_REAL4 *single_ptr;
  CCTK_REAL *real_ptr;
  CCTK_CHAR *char_ptr;
  CCTK_INT *int_ptr;
  void *data = CCTK_VarDataPtrI (GH, timelevel, index);

  /* to make the compiler happy */
  single_ptr = NULL;
  real_ptr = NULL;
  char_ptr = NULL;
  int_ptr = NULL;

  ioUtilGH = (ioGH *) GH->extensions [CCTK_GHExtensionHandle ("IO")];
  pughGH = (pGH *) GH->extensions [CCTK_GHExtensionHandle ("PUGH")];

  myproc = CCTK_MyProc (GH);

  if (ioUtilGH->downsample_x == 1 &&
      ioUtilGH->downsample_y == 1 &&
      ioUtilGH->downsample_z == 1) {

    if (ioUtilGH->out_single) {
      single_ptr = (CCTK_REAL4 *) malloc (pughGH->npoints*sizeof (CCTK_REAL4));

      for (i = 0; i < pughGH->npoints; i++)
        single_ptr [i] = (CCTK_REAL4) ((CCTK_REAL *) data) [i];

      *outme = single_ptr;
      *free_outme = 1;
    } else {
      *outme = data;
      *free_outme = 0;
    }

    for (i = 0; i < 3; i++) {
      bnd [i] = GH->cctk_lbnd[i];     /* the bounds */
      bnd [i+3] = GH->cctk_lsh[i];   /* the sizes */
      bnd [i+6] = GH->cctk_gsh[i];  /* the global space */
    }

  } else {

    int start [3], end [3];
    int i, j, k, l;

    /* Downsampling code ... */
    bnd [6] = GH->cctk_gsh[0] / ioUtilGH->downsample_x;
    if (GH->cctk_gsh[0] % ioUtilGH->downsample_x)
      bnd [6]++;
    bnd [7] = GH->cctk_gsh[1] / ioUtilGH->downsample_y;
    if (GH->cctk_gsh[1] % ioUtilGH->downsample_y)
      bnd [7]++;
    bnd [8] = GH->cctk_gsh[2] / ioUtilGH->downsample_z;
    if (GH->cctk_gsh[2] % ioUtilGH->downsample_z)
      bnd [8]++;

    if (verbose) 
      printf ("Downsampled sizes (%d, %d, %d) -> (%d, %d, %d)\n",
              GH->cctk_gsh[0], GH->cctk_gsh[1], GH->cctk_gsh[2],
              (int) bnd [6], (int) bnd [7], (int) bnd [8]);
    
    /* Now figure out the local downsampling */
    /* The local starts are the lb modded into the downsample */
    for (i = 0; i < 3; i++) {
      int downsample;

      if (i == 0)
        downsample = ioUtilGH->downsample_x;
      else if (i == 1)
        downsample = ioUtilGH->downsample_y;
      else
        downsample = ioUtilGH->downsample_z;

      bnd [i]   = GH->cctk_lbnd[i] / downsample;
      start [i] = bnd [i] * downsample;
      if (start [i] <
          GH->cctk_lbnd[i] + pughGH->ownership [PUGH_VERTEXCTR][i][0]) {
        start [i] += downsample;
        bnd [i] ++;
      }
      end [i]   = ((GH->cctk_lbnd [i] +
                   pughGH->ownership [PUGH_VERTEXCTR][i][1] - 1) / downsample)
                   * downsample;
      bnd [i+3] = (end [i] - start [i]) / downsample + 1;
    }

    if (verbose) {
      printf ("Downsample ranges (%d, %d, %d) -> (%d, %d, %d)\n",
        start [0], start [1], start [2],
        end [0], end [1], end [2]);
      printf ("Local size/bound  (%d, %d, %d)  (%d, %d, %d)\n",
        (int) bnd [3], (int) bnd [4], (int) bnd [5],
        (int) bnd [0], (int) bnd [1], (int) bnd [2]);
    }

    /* compute local ranges */
    for (i = 0; i < 3; i++) {
      start [i] -= GH->cctk_lbnd [i];
      end [i] -= GH->cctk_lbnd [i];
    }

    *outme = malloc (bnd [3] * bnd [4] * bnd [5] * element_size);
    *free_outme = 1;

    /* I hate it to repeat the loops for each case label
       but that way produces much more efficient code */
    l = 0;
    switch (CCTK_VarTypeI (index)) {
      case CCTK_VARIABLE_CHAR:
        char_ptr = (CCTK_CHAR *) *outme;
        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample_z)
          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample_y)
            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample_x)
              char_ptr [l++] = ((CCTK_CHAR *) data) [DI (pughGH, i, j, k)];
        break;

      case CCTK_VARIABLE_INT:
        int_ptr = (CCTK_INT *) *outme;
        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample_z)
          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample_y)
            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample_x)
              int_ptr [l++] = ((CCTK_INT *) data) [DI (pughGH, i, j, k)];
        break;

      case CCTK_VARIABLE_REAL:
        if (ioUtilGH->out_single)
          single_ptr = (CCTK_REAL4 *) *outme;
        else
          real_ptr = (CCTK_REAL *) *outme;
        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample_z)
          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample_y)
            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample_x)
              if (ioUtilGH->out_single)
                single_ptr [l++] = (CCTK_REAL4)
                                  (((CCTK_REAL *) data) [DI (pughGH, i, j, k)]);
              else
                real_ptr [l++] = ((CCTK_REAL *) data) [DI (pughGH, i, j, k)];
        break;

      default:
        CCTK_WARN (1, "Unsupported variable type in IOPanda_getDumpData");
        return;
    }
  }

  if (verbose) {
    printf ("Global size: %d %d %d\n",
            (int) bnd [6], (int) bnd [7], (int) bnd [8]);
    printf ("Lower bound: %d %d %d\n",
            (int) bnd [0], (int) bnd [1], (int) bnd [2]);
    printf ("Chunk size : %d %d %d\n",
            (int) bnd [3], (int) bnd [4], (int) bnd [5]);
  }
}