testxmdf.c

testxmdf.c provides a Parallel MPI C code example and performance measurements for scalar dataset writes.

/******************************************************************************
 * testxmdf.c : Defines the entry point for the parallel XMDF tests.
 ******************************************************************************/

/*
 * INCLUDE HEADERS
 */

#include "xmdf_putil.h"
#include <mpi.h>
#include <stdlib.h>
#include <stdio.h>
#include "xmdf_par.h"

#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>             /* ,lseek */
#include <string.h>

#include "closeh5.h"            /* get close library */

/*
 * DEFINITIONS
 */

#define ERRORD(A,B,C) ERROR(A,B,C)
#define PRINTM(A) PRINTS("msg",(A))

/* FUNCTION  fileinit
 * PURPOSE   initializes xmdf file
 * NOTES
 */
xid fileinit( char *a_fname, char *a_gname )
{
  xid fileID, dataSetID;
  closeObjPtr closeobj = closeCreate();
  xid ret;

  xfiSetComm( MPI_COMM_WORLD );
  
  ret = xfCreateFilePar( a_fname, &fileID, XTRUE, MPI_INFO_NULL );
  if (0 > closePush( closeobj, xfCloseFile, ret < 0 ? ret : fileID ))
    {
      PRINTM( "xfCreateFilePar returned error" );
      return ret;
    }

  /* xfOpenFilePar( a_fname, &fileId, XFALSE, MPI_INFO_NULL ); */

  ret = xfCreateScalarDataset( fileID, a_gname, "", TS_SECONDS, -1, &dataSetID );
  if (0 > closePush( closeobj, xfCloseGroup, ret < 0 ? ret : dataSetID ))
    {
      PRINTM( "xfCreateScalarDataset returned error" );
      return ret;
    }

  closePush( closeobj, NULL, -1 ); /* close all */
  return 0;                     /* return success */
}

/* FUNCTION  setup
 * PURPOSE   reopens file and returns dataset id and totalcount
 * NOTES
 */
xid setup (closeObjPtr a_closeobj, /* stack close calls */
           const char *a_filename,
           const char *a_groupname,
           const int a_mynum,     /* number of values /process/timestep */
           int *a_totalCount)     /* number of values /timestep */
{
  xid ret;
  xid fileID, dataSetID;

  /* open file for write */
  ret = xfOpenFilePar(a_filename, &fileID, XFALSE, MPI_INFO_NULL);
  if (0 > closePush( a_closeobj, xfCloseFile, ret < 0 ? ret : fileID ))
    {
      ERRORD( "setup", "xfOpenFilePar returned error", ret );
      return ret;
    }

  PRINTS( "filename", a_filename );
  PRINTS( "groupname", a_groupname );

  /* open group */
  ret = xfOpenGroup(fileID, a_groupname, &dataSetID);
  if (0 > closePush( a_closeobj, xfCloseGroup, ret < 0 ? ret : dataSetID ))
    {
      ERRORD( "setup", "xfOpenGroup returned error", ret );
      return ret;
    }

  /* get totalCount sum of all mynum values */
  {
    int num = a_mynum;
    MPI_Allreduce( &num,
                   a_totalCount,
                   1,
                   MPI_INT,
                   MPI_SUM,
                   MPI_COMM_WORLD );
  }
  return dataSetID;
}

/* FUNCTION  setup2
 * PURPOSE   reopens file and returns dataset id and totalcount
 * NOTES
 */
xid setup2 (closeObjPtr a_closeobj, /* stack close calls */
           const char *a_filename,
           const char *a_groupname,
           const int a_mynum,     /* number of values /process/timestep */
           int *a_totalCount)     /* number of values /timestep */
{
  xid ret;
  xid fileID, dataSetID;

  /* open file for write */
  ret = xfOpenFilePar(a_filename, &fileID, XFALSE, MPI_INFO_NULL);
  if (0 > closePush( a_closeobj, xfCloseFile, ret < 0 ? ret : fileID ))
    {
      ERRORD( "setup", "xfOpenFilePar returned error", ret );
      return ret;
    }

  PRINTS( "filename", a_filename );
  PRINTS( "groupname", a_groupname );

  /* open group */
  ret = xfOpenGroup(fileID, a_groupname, &dataSetID);
  if (0 > closePush( a_closeobj, xfCloseGroup, ret < 0 ? ret : dataSetID ))
    {
      ERRORD( "setup", "xfOpenGroup returned error", ret );
      return ret;
    }

  /* get totalCount sum of all mynum values */
  {
    int num = a_mynum;
    MPI_Allreduce( &num,
                   a_totalCount,
                   1,
                   MPI_INT,
                   MPI_SUM,
                   MPI_COMM_WORLD );
  }
  return dataSetID;
}

/* FUNCTION  abs
 * PURPOSE   return absolute value of val
 * NOTES     
 */
double absd( double a_val )
{
  return a_val < 0 ? -a_val : a_val;
}

/* FUNCTION  testminmax
 * PURPOSE   compares array of values against xfGetDatasetMins/Maxs array
 * NOTES     returns <0 on miss compare, and prints error
 */
int testminmax( int a_timesteps,  /* number of timesteps in file */
                xid a_dataSetID,  /* dataset in file */
                int a_minnotmax,  /* XTRUE is minimum test, not max */
                double a_tolerance, /* tolerance of compare */
                double *a_expectminmaxs ) /* array of expected mins or maxs */
{
  xid ret;
  int ii;
  float *minmaxs = malloc(sizeof(float)*a_timesteps);

  if (!minmaxs)
    {
      ERRORD( "testScalar", "malloc", 0 );
      return -1;
    }

  TIMER("pre");
  ret = a_minnotmax ?
    xfGetDatasetMins( a_dataSetID, a_timesteps, minmaxs ) :
    xfGetDatasetMaxs( a_dataSetID, a_timesteps, minmaxs );
  TIMER("TIME xfGetDatasetMins/Maxs");
  if (ret < 0)
    {
      free(minmaxs);
      ERRORD( "testminmax",
              a_minnotmax ?
              "xfGetDatasetMins" :
              "xfGetDatasetMaxs",
              ret );
      return ret;
    }

  for (ii=0; ii<a_timesteps; ii++)
    {
      if (absd(minmaxs[ii]) > absd(a_expectminmaxs[ii]) * (1+a_tolerance) ||
          absd(minmaxs[ii]) < absd(a_expectminmaxs[ii]) * (1-a_tolerance) )
        {
          ERRORD( "testminmax",
                  a_minnotmax ?
                  "read min does not match parallel write at timestep" :
                  "read max does not match parallel write at timestep",
                  ii );
          PRINTF( "found", minmaxs[ii] );
          PRINTF( "expect", a_expectminmaxs[ii] );
          PRINTF( "difference", a_expectminmaxs[ii]-minmaxs[ii] );
          return -1;
        }
    }

  if (minmaxs) free(minmaxs);

  return !0;
}

/* FUNCTION  testScalar
 * PURPOSE   Runs the scalar parallel XMDF test
 * NOTES
 */
int testScalar()
{
  double *values;               /* pointer to data value buffer */
  int *indices;                 /* index of each data item (parallel issue) */
  int subnum;                   /* number of values written by each process */
  char *filename = "dbl_test.h5";
  char *groupname = "myData";
  int timesteps = 2;
  double average[2];
  int rank;
  int offset;
  double writemins[2];
  double writemaxs[2];

  /* create file */
  {
    xid ret = fileinit(filename, groupname);
    if (0 > ret)
      return ret;
  }

  /* init subnum */
  subnum = 1000000;
  PRINTD( "subnum", subnum );

  values = malloc( sizeof(double) * subnum); /* allocate data buffer */
  indices = malloc( sizeof(int) * subnum );

  /* init indices and values */
  {
    MPI_Comm_rank( MPI_COMM_WORLD, &rank ); /* get process number */

    /* Returns in offset the sum of this process's subnum and
     * all subnum's from processes of smaller rank, thus
     * rank 0 gets subnum*1, rank 1 gets subnum*2, etc */
    MPI_Scan( &subnum, &offset, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
    offset -= subnum;                /* start at 0, not at subnum */

  }

  /* reopen and write file with scalar */
  {
    int ii;
    xid ret;
    xid dataSetID;
    int totalnum;
    closeObjPtr closeobj = closeCreate();

    dataSetID = setup( closeobj, filename, groupname, subnum, &totalnum );
    if (dataSetID < 0) return dataSetID;

    PRINTF( "totalnum", totalnum );

    /* write timesteps */
    for (ii=0; ii<timesteps; ii++)
      {
        /* init values and indices */
        {
          int jj;
          int pick = 777;

          average[ii] = 0;

          for (jj = 0; jj<subnum; jj++)
            {
              /* set index and value for each entry */
              indices[jj] = (offset + jj - 3) % totalnum + 1;
              values[jj] = indices[jj] + (rank * .01) + (ii * .001);

              /* calculate test average, mins, maxs, and test pick */
              average[ii] += values[jj]/subnum;
              if (jj==0 || writemins[ii]>values[jj]) writemins[ii] = values[jj];
              if (jj==0 || writemaxs[ii]<values[jj]) writemaxs[ii] = values[jj];
              if( (indices[jj]-1)%subnum == 0 ) pick = jj;
            }

          PRINTD( "offset", offset );
          PRINTD( "pick", pick );
          PRINTD( "indices[pick]-1", indices[pick]-1 );
          PRINTF( "values[pick]", values[pick] );
          average[ii] *= values[pick]; 
        }

        TIMER("pre");
        xfiPrintDataset( __FILE__, __LINE__, "test", dataSetID);
        PRINTD( "totalnum", totalnum );
        PRINTD( "subnum", subnum );
        PRINTF( "values", values[0] );
        PRINTD( "indices", indices[0] );
        ret = xfWriteScalarTimestepPar( dataSetID,
                                        0.0, /* time */
                                        totalnum,
                                        subnum,
                                        values,
                                        indices,
                                        0 ); /* stride */
        TIMER("TIME xfWriteScalarTimestepPar");
        /* xid          a_Id        xid of the group containing the dataset
         * double       a_Time      time of the current time step
         * int          a_NumTotal  total number of values in the dataset
         * int          a_NumValues number of values to be written by this processor
         * const double a_Values    values to be written to the dataset
         * const int    a_WhichVals indices corresponding to the location of the
         *                          values in the dataset
         * int          a_Stride    number of values (or bytes) to be skipped in the
         *                          dataset after each value
         */

        if (closePush( closeobj, NULL, ret ) < 0)
          {
            ERRORD( "testScalar", "xfWriteScalarTimestepPar", ret );
            return ret;
          }
      }

    closePush( closeobj, NULL, -1 ); /* close all */
  }

  /* reopen and read file */
  {
    int ii;
    xid ret;
    xid dataSetID;
    int totalnum;
    closeObjPtr closeobj = closeCreate();

    dataSetID = setup( closeobj, filename, groupname, subnum, &totalnum );
    if (dataSetID < 0) return dataSetID;

    for (ii=0; ii<timesteps; ii++)
      {
        TIMER("pre");
        PRINTD("dataSetID",(int)dataSetID);
        xfiPrintDataset( __FILE__, __LINE__, "read", dataSetID );
        {
          xid       DatasetId = H5Dopen(dataSetID, "Values");
          printf( "%s:%d: H5Dopen %d\n", __FILE__, __LINE__, (int)DatasetId );
          H5Dclose(DatasetId);
        }
        ret = xfReadScalarValuesTimestepPar( dataSetID,
                                             ii+1, /* timestep */
                                             offset+1, /* start */
                                             subnum,
                                             values,
                                             0 );
        {
          xid       DatasetId = H5Dopen(dataSetID, "Values");
          PRINTD( "H5Dopen", DatasetId );
          H5Dclose(DatasetId);
        }
        TIMER("TIME xfReadScalarValuesTimestepPar");
        /* xid   a_Id          xid of the group containing the dataset
         * int   a_Time        time index of the current time step (1 based)
         * int   a_StartIndex  indices corresponding to the starting location of the
         *                     values in the dataset to be read
         * int   a_NumValues   number of values to be read
         * float a_Values      values to be read from the dataset
         * int   a_Stride      number of values (or bytes) to be skipped in the
         *                     dataset after each value
         */
        if (closePush( closeobj, NULL, ret ) < 0)
          {
            printf( "%s:%d: ERROR in %s: %d %d %d %d %d %d %d\n",
                    __FILE__, __LINE__,
                    "testScalar", ret,
                    (int)dataSetID,
                    ii+1, /* timestep */
                    offset+1, /* start */
                    subnum,
                    values==0,
                    0 );
            return ret;
          }

        /* compare results */
        {
          double avg = 0;
          int jj;

          for (jj = 0; jj<subnum; jj++)
            {
              avg += values[jj]/subnum;
            }
          avg *= values[0];     /* verify ordering (crc would be better) */

          if (avg > average[ii]*(1+.00001) ||
              avg < average[ii]*(1-.00001) )
            {
              ERRORD( "testScalar", "read does not match write for timestep", ii+1 );
              PRINTF( "found", avg );
              PRINTF( "expected", average[ii] );
              PRINTF( "values[0]", values[0] );
            }
          else
            PRINTD( "PASS: parallel read matches parallel write, step ", ii );
        }
      }

    closePush( closeobj, NULL, -1 ); /* close all */
  }

  /* reopen and read file sequentially via processor 0 */
  MPI_Barrier(MPI_COMM_WORLD);
  if (!rank) {
    int ii;
    xid ret;
    xid dataSetID;
    xid fileID;
    closeObjPtr closeobj = closeCreate();

    /* open file for write */
    ret = xfOpenFile(filename, &fileID, XTRUE);
    if (0 > closePush( closeobj, xfCloseFile, ret < 0 ? ret : fileID ))
      {
        PRINTM( "xfOpenFilePar returned error" );
        return ret;
      }

    /* open group */
    ret = xfOpenGroup(fileID, groupname, &dataSetID);
    if (ret >= 0) ret = dataSetID;
    if (0 > closePush( closeobj, xfCloseGroup, ret ))
      {
        PRINTM( "xfOpenGroup returned error" );
        return ret;
      }

    ret = testminmax( timesteps, dataSetID, XTRUE, 0.001, writemins );
    if (ret < 0)
      {
        ERRORD( "testScalar", "testminmax (min)", ret );
        return ret;
      }

    ret = testminmax( timesteps, dataSetID, XFALSE, 0.001, writemaxs );
    if (ret < 0)
      {
        ERRORD( "testScalar", "testminmax (max)", ret );
        return ret;
      }

    for (ii=0; ii<timesteps; ii++)
      {
        TIMER("pre");
        ret = xfReadScalarValuesTimestepDoublePortion( dataSetID,
                                                       ii+1, /* time */
                                                       1,
                                                       subnum,
                                                       values );
        TIMER("TIME xfReadScalarValuesTimestepDoublePortion");
        /* xid   a_Id          xid of the group containing the dataset
         * int   a_Time        time index of the current time step (1 based)
         * int   a_StartIndex  indices corresponding to the starting location of the
         *                     values in the dataset to be read
         * int   a_NumValues   number of values to be read
         * float a_Values      values to be read from the dataset
         */
        if (closePush( closeobj, NULL, ret ) < 0)
          {
            PRINTD( "ERROR in testScalar:  xfReadScalarTimestepDoublePortion returned", ret );
            return ret;
          }

        /* compare results */
        {
          double avg = 0;
          int jj;

          for (jj = 0; jj<subnum; jj++)
            {
              avg += values[jj]/subnum;
            }
          avg *= values[0];     /* verify ordering (crc would be better) */

          if (avg > average[ii]*(1+.00001) ||
              avg < average[ii]*(1-.00001) )
            {
              ERRORD( "testScalar", "read does not match write for timestep", ii+1 );
              PRINTF( "found", avg );
              PRINTF( "expected", average[ii] );
            }
          else
            PRINTM( "PASS: sequential read matches parallel write" );
        }
      }

    closePush( closeobj, NULL, -1 ); /* close all */
  }

  return !0;
}

/* FUNCTION  testVector
 * PURPOSE   Runs the vector parallel XMDF test
 * NOTES
 */
int testVector()
{
  double *values;               /* pointer to data value buffer */
  int valuesRank = 2;           /* number of dimensions of the data */
  int *indices;                 /* index of each data item (parallel issue) */
  int subnum;                   /* number of values written by each process */
  char *filename = "dbl_test.h5"; /* OS file name */
  char *groupname = "myData";   /* hdf5 group name */

  /* create file */
  {
    xid ret = fileinit(filename, groupname);
    if (0 > ret)
      return ret;
  }

  /* init subnum */
  subnum = 100;
  PRINTD( "subnum", subnum );

  values = malloc( sizeof(*values) * subnum * valuesRank ); /* allocate data buffer */
  indices = malloc( sizeof(*indices) * subnum );

  /* init indices and values */
  {
    int rank;
    int offset;

    MPI_Comm_rank( MPI_COMM_WORLD, &rank ); /* get process number */

    /* Returns in offset the sum of this process's subnum and
     * all subnum's from processes of smaller rank, thus
     * rank 0 gets subnum*1, rank 1 gets subnum*2, etc */
    MPI_Scan( &subnum, &offset, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
    offset -= subnum;                /* start at 0, not at subnum */

    {
      int ii;
      int jj;

      for (ii = 0; ii<subnum; ii++)
        {
          indices[ii] = offset + ii + 1; /* set index of each data value */

          for (jj = 0; jj<valuesRank; jj++)
            {
              /* init data in buffer */
              values[ii*valuesRank + jj] = (ii*10 + jj)*10 + rank;
            }
        }
    }
  }

  /* reopen and write file with vector */
  {
    xid ret;
    xid dataSetID;
    int totalnum;
    closeObjPtr closeobj = closeCreate();

    dataSetID = setup( closeobj, filename, groupname, subnum, &totalnum );
    if (dataSetID < 0) return dataSetID;

    ret = xfWriteVectorTimestepPar( dataSetID,
                                    0.0,
                                    totalnum,
                                    subnum,
                                    valuesRank,
                                    values,
                                    indices,
                                    0 );
    /* xid          a_Id        xid of the group containing the dataset
     * double       a_Time      time of the current time step
     * int          a_NumTotal  total number of values in the dataset
     * int          a_NumProc   number of values to be written by this processor
     * int          a_NumComponents number of components in each vector value
     * const double*a_Values    values to be written to the dataset
     * const int   *a_WhichVals indices corresponding to the location of the
     *                          values in the dataset
     * int          a_Stride    number of values (or bytes) to be skipped in the
     *                          dataset after each value
     */

    /* xfWriteVectorTimestepPar does the following:
     *
     * return error if number of process values is less than 1
     * copy data to a new array if stride is set ??? expand
     * flag whether using null value ??? expand
     * calculate min and max vector lengths
     * open datasets for data, time, min, max
     * write time after waiting for other processes, or append
     * write data after waiting for other processes, or append ??? expand
     * . create dataset
     * . select elements
     * write min/max after waiting for other processes, or append
     */

    if (closePush( closeobj, NULL, ret ) < 0)
      {
        PRINTD( "ERROR in testVector: xfWriteVectorTimestepPar returned", ret );
        return ret;
      }

    printf( "%s:%d: wrote values[5][0-1] = %f %f\n", __FILE__, __LINE__,
            values[5*valuesRank + 0],
            values[5*valuesRank + 1] );

    closePush( closeobj, NULL, -1 ); /* close all */
  }

  /* reopen and read file with vector */
  {
    xid ret;
    xid dataSetID;
    int totalnum;
    closeObjPtr closeobj = closeCreate();

    dataSetID = setup( closeobj, filename, groupname, subnum, &totalnum );
    if (dataSetID < 0) return dataSetID;

    ret = xfReadScalarValuesTimestepPar( dataSetID,
                                         1,
                                         1,
                                         subnum,
                                         values,
                                         0 );
    /* xid   a_Id          xid of the group containing the dataset
     * int   a_Time        time index of the current time step (1 based)
     * int   a_StartIndex  indices corresponding to the starting location of the
     *                     values in the dataset to be read
     * int   a_NumValues   number of values to be read
     * float a_Values      values to be read from the dataset
     * int   a_Stride      number of values (or bytes) to be skipped in the
     *                     dataset after each value
     */
    if (closePush( closeobj, NULL, ret ) < 0)
      {
        PRINTD( "ERROR in testVector: xfReadVectorTimestepPar returned", ret );
        return ret;
      }

    printf( "%s:%d: wrote values[5][0-1] = %f %f\n", __FILE__, __LINE__,
            values[5*valuesRank + 0],
            values[5*valuesRank + 1] );

    closePush( closeobj, NULL, -1 ); /* close all */
  }

  return !0;
}

void timeposix()
{
  int fd;
  char *fname = "data.posix";
  int numentries = 3000000;
  int memsize = sizeof(double) * numentries;
  int uniqueloc = numentries/3;
  double *ptr = (double*)malloc(memsize);
  int ret;
  int offset;
  double uniquedata;
  
  uniquedata = deltatime();

  /* WRITE */

  fd=open(fname, O_WRONLY|O_CREAT|O_APPEND, 0760); /* open file for writing */
  if (fd < 0 )
    {
      printf("%s:%d: ERROR: Failed to open file %s\n", __FILE__,__LINE__, fname);
      return;
    }
  
  memset(ptr, 0, memsize);      /* clear buffer */
  ptr[uniqueloc] = uniquedata;     /* set unique value */
  TIMER("PreIsoWrite");
  ret = write(fd, ptr, memsize); /* write buffer */
  offset = lseek(fd, 0, SEEK_CUR);       /* get offset into file */
  close(fd);                    /* close file */
  TIMER("POSTiSOwRITE");
  if (ret != memsize)           /* check write failure */
    {
      printf("%s:%d: ERROR: Failed to write file %s\n", __FILE__,__LINE__, fname);
      return;
    }

  /* READ */

  memset(ptr, 0, memsize);      /* clear buffer */
  if (ptr[uniqueloc] == memsize) /* verify cleared buffer */
    {
      printf("%s:%d: ERROR: buffer was not cleared\n", __FILE__,__LINE__);
      return;
    }

  fd=open(fname, O_RDONLY, 0760); /* open file for reading */
  if (fd < 0 )
    {
      printf("%s:%d: ERROR: Failed to open readable file %s\n", __FILE__,__LINE__, fname);
      return;
    }
  
  ret = lseek(fd, offset-memsize, SEEK_CUR); /* offset to where data was written */
  if (ret != offset-memsize)
    {
      printf("%s:%d: ERROR: lseek return value is off by %d\n", __FILE__,__LINE__, offset - memsize - ret);
      close(fd);
      return;
    }

  TIMER("PreIsoRead");
  ret = read(fd, ptr, memsize);
  TIMER("POSTiSOrEAD");
  close(fd);
  if (ret != memsize)
    {
      printf("%s:%d: ERROR: Failed to read file %s\n", __FILE__,__LINE__, fname);
      return;
    }

  if (ptr[uniqueloc] != uniquedata) /* verify cleared buffer */
    {
      printf("%s:%d: ERROR: data incorrect.\n", __FILE__,__LINE__);
      return;
    }

  printf("%s:%d: PASS PASS PASS.\n", __FILE__,__LINE__);
}

/* FUNCTION  main
 * PURPOSE   runs the test cases for parallel XMDF
 * NOTES
 */
int main(int argc, char **argv)
{
  int fail = !0;
  int ret;

  MPI_Init(&argc, &argv);       /* init multiprocessor interface */
  TIMER(NULL);                  /* initialize timer */

  MPI_Barrier(MPI_COMM_WORLD);
  timeposix();
  MPI_Barrier(MPI_COMM_WORLD);

  ret = testScalar(); if (ret < 0) fail = ret;
  TIMER("final time");          /* initialize timer */
  PRINTD( ret < 0 ? "FAIL" : "ALL PASS", fail );

  xfClosePar();                 /* must be closed before MPI */
  MPI_Finalize();
  return fail;
}

---