matrixio.c

MIT开发出来的计算光子晶体的软件

/* Copyright (C) 1999, 2000, 2001, 2002, Massachusetts Institute of Technology.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/* matrixio.c: This file layers a "matrixio" abstraction on top of the
   HDF5 binary i/o interface.  This abstraction should make HDF5 much
   easier to use for our purposes, and could also conceivably allow
   us to replace HDF5 with some other file format (e.g. HDF4). */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include "../config.h"

#include <check.h>
#include <scalar.h>
#include <mpiglue.h>
#include <mpi_utils.h>

#include "matrixio.h"

/*****************************************************************************/

/* If we have the H5Pset_fapl_mpio function (which is available if HDF5 was
   compiled for MPI), then we can perform collective file i/o operations
   (e.g. all processes call H5Fcreate at the same time to create one file).

   If we don't, however, then we deal with it by having one process
   work with the file at a time: "exclusive" access.  The following
   macro helps us select different bits of code depending upon whether
   this is the case. */

#ifdef HAVE_H5PSET_MPI /* old name for this routine */
#  define H5Pset_fapl_mpio H5Pset_mpi
#  ifndef HAVE_H5PSET_FAPL_MPIO
#    define HAVE_H5PSET_FAPL_MPIO 1
#  endif
#endif

#ifdef HAVE_H5PSET_FAPL_MPIO
#  define IF_EXCLUSIVE(yes,no) no
#else
#  define IF_EXCLUSIVE(yes,no) yes
#endif

/*****************************************************************************/

/* Normally, HDF5 prints out all sorts of error messages, e.g. if a dataset
   can't be found, in addition to returning an error code.  The following
   macro can be wrapped around code to temporarily suppress error messages. */

#define SUPPRESS_HDF5_ERRORS(statements) { \
     H5E_auto_t xxxxx_err_func; \
     void *xxxxx_err_func_data; \
     H5Eget_auto(&xxxxx_err_func, &xxxxx_err_func_data); \
     H5Eset_auto(NULL, NULL); \
     { statements; } \
     H5Eset_auto(xxxxx_err_func, xxxxx_err_func_data); \
}

/*****************************************************************************/

/* Wrappers to write/read an attribute attached to id.  HDF5 attributes
   can *not* be attached to files, in which case we'll write/read it
   as an ordinary dataset.  Ugh. */

static void write_attr(matrixio_id id, matrixio_id type_id,
		       matrixio_id space_id, const char *name, void *val)
{
#if defined(HAVE_HDF5)
     hid_t attr_id;

     if (!mpi_is_master())
	  return; /* only one process should add attributes */
     
     if (H5I_FILE == H5Iget_type(id)) {
          attr_id = H5Dcreate(id, name, type_id, space_id, H5P_DEFAULT);
          CHECK(id >= 0, "error creating HDF attr");
          H5Dwrite(attr_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, val);
          H5Dclose(attr_id);
     }
     else {
	  attr_id = H5Acreate(id, name, type_id, space_id, H5P_DEFAULT);
	  CHECK(id >= 0, "error creating HDF attr");
	  H5Awrite(attr_id, type_id, val);
	  H5Aclose(attr_id);
     }
#endif
}

static matrixio_id open_attr(matrixio_id id, matrixio_id *type_id,
			     matrixio_id *space_id, const char *name)
{
#if defined(HAVE_HDF5)
     hid_t attr_id;

     if (H5I_FILE == H5Iget_type(id)) {
          SUPPRESS_HDF5_ERRORS(attr_id = H5Dopen(id, name));
	  if (attr_id >= 0) {
	       *type_id = H5Dget_type(attr_id);
	       *space_id = H5Dget_space(attr_id);
	  }
     }
     else {
          SUPPRESS_HDF5_ERRORS(attr_id = H5Aopen_name(id, name));
	  if (attr_id >= 0) {
	       *type_id = H5Aget_type(attr_id);
	       *space_id = H5Aget_space(attr_id);
	  }
     }

     return attr_id;
#else
     return -1;
#endif
}

static void read_attr(matrixio_id id, matrixio_id attr_id,
		      matrixio_id mem_type_id, void *val)
{
#if defined(HAVE_HDF5)
     if (H5I_FILE == H5Iget_type(id)) {
	  H5Dread(attr_id, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, val);
     }
     else {
	  H5Aread(attr_id, mem_type_id, val);
     }
#endif
}

static void close_attr(matrixio_id id, matrixio_id attr_id)
{
#if defined(HAVE_HDF5)
     if (H5I_FILE == H5Iget_type(id)) {
          H5Dclose(attr_id);
     }
     else {
          H5Aclose(attr_id);
     }
#endif
}

/*****************************************************************************/

void matrixio_write_string_attr(matrixio_id id, const char *name,
				const char *val)
{
#if defined(HAVE_HDF5)
     hid_t type_id;
     hid_t space_id;

     if (!val || !name || !name[0] || !val[0])
	  return; /* don't try to create empty attributes */
     
     type_id = H5Tcopy(H5T_C_S1);
     H5Tset_size(type_id, strlen(val) + 1);
     space_id = H5Screate(H5S_SCALAR);
     write_attr(id, type_id, space_id, name, (void*) val);
     H5Sclose(space_id);
     H5Tclose(type_id);
#endif
}

void matrixio_write_data_attr(matrixio_id id, const char *name,
			      const real *val, int rank, const int *dims)
{
#if defined(HAVE_HDF5)
     hid_t type_id;
     hid_t space_id;
     hsize_t *space_dims;
     int i;

     if (!val || !name || !name[0] || rank < 0 || !dims)
	  return; /* don't try to create empty attributes */
     
#ifdef SCALAR_SINGLE_PREC
     type_id = H5T_NATIVE_FLOAT;
#else
     type_id = H5T_NATIVE_DOUBLE;
#endif

     if (rank > 0) {
	  CHK_MALLOC(space_dims, hsize_t, rank);
	  for (i = 0; i < rank; ++i)
	       space_dims[i] = dims[i];
	  space_id = H5Screate_simple(rank, space_dims, NULL);
	  free(space_dims);
     }
     else {
	  space_id = H5Screate(H5S_SCALAR);
     }

     write_attr(id, type_id, space_id, name, (void*) val);
     H5Sclose(space_id);
#endif
}

char *matrixio_read_string_attr(matrixio_id id, const char *name)
{
#if defined(HAVE_HDF5)
     hid_t attr_id;
     hid_t type_id;
     hid_t space_id;
     int len;
     char *s = NULL;

     if (!name || !name[0])
	  return NULL; /* don't try to read empty-named attributes */
     
     attr_id = open_attr(id, &type_id, &space_id, name);
     if (attr_id < 0)
	  return NULL;

     if (H5Sget_simple_extent_npoints(space_id) == 1) {
	  len = H5Tget_size(type_id);
	  H5Tclose(type_id);
	  
	  type_id = H5Tcopy(H5T_C_S1);
	  H5Tset_size(type_id, len);
	  
	  CHK_MALLOC(s, char, len);
	  read_attr(id, attr_id, type_id, (void*) s);
     }

     H5Tclose(type_id);
     H5Sclose(space_id);
     close_attr(id, attr_id);

     return s;
#else
     return NULL;
#endif
}

real *matrixio_read_data_attr(matrixio_id id, const char *name,
			      int *rank, int max_rank, int *dims)
{
#if defined(HAVE_HDF5)
     hid_t attr_id, type_id, mem_type_id, space_id;
     real *d = NULL;

     if (!name || !name[0] || max_rank < 0 || !dims)
	  return NULL; /* don't try to create empty attributes */
     
#ifdef SCALAR_SINGLE_PREC
     mem_type_id = H5T_NATIVE_FLOAT;
#else
     mem_type_id = H5T_NATIVE_DOUBLE;
#endif

     attr_id = open_attr(id, &type_id, &space_id, name);
     if (attr_id < 0)
          return NULL;

     *rank = H5Sget_simple_extent_ndims(space_id);
     if (*rank <= max_rank) {
	  if (*rank > 0) {
	       int i;
	       hsize_t *space_dims, *maxdims;
	       CHK_MALLOC(space_dims, hsize_t, *rank);
	       CHK_MALLOC(maxdims, hsize_t, *rank);
	       H5Sget_simple_extent_dims(space_id, space_dims, maxdims);
	       for (i = 0; i < *rank; ++i)
		    dims[i] = space_dims[i];
	       free(maxdims);
	       free(space_dims);
	  }
	  CHK_MALLOC(d, real, H5Sget_simple_extent_npoints(space_id));
          read_attr(id, attr_id, mem_type_id, (void*) d);
     }

     H5Tclose(type_id);
     H5Sclose(space_id);
     close_attr(id, attr_id);

     return d;
#else
     return NULL;
#endif
}

/*****************************************************************************/

#define FNAME_SUFFIX ".h5"  /* standard HDF5 filename suffix */

static char *add_fname_suffix(const char *fname)
{
     int oldlen = strlen(fname);
     int suflen = strlen(FNAME_SUFFIX);
     char *new_fname;

     CHECK(fname, "null filename!");

     CHK_MALLOC(new_fname, char, oldlen + suflen + 1);

     strcpy(new_fname, fname);

     /* only add suffix if it is not already there: */
     if (strstr(new_fname, FNAME_SUFFIX) != new_fname + oldlen - suflen)
	  strcat(new_fname, FNAME_SUFFIX);

     return new_fname;
}

/*****************************************************************************/

static int matrixio_critical_section_tag = 0;

matrixio_id matrixio_create(const char *fname)
{
#if defined(HAVE_HDF5)
     char *new_fname;
     matrixio_id id;
     hid_t access_props;

     access_props = H5Pcreate (H5P_FILE_ACCESS);
     
#  if defined(HAVE_MPI) && defined(HAVE_H5PSET_FAPL_MPIO)
     H5Pset_fapl_mpio(access_props, MPI_COMM_WORLD, MPI_INFO_NULL);
#  endif

     new_fname = add_fname_suffix(fname);

#  ifdef HAVE_H5PSET_FAPL_MPIO
     id = H5Fcreate(new_fname, H5F_ACC_TRUNC, H5P_DEFAULT, access_props);
#  else
     mpi_begin_critical_section(matrixio_critical_section_tag);
     if (mpi_is_master())
	  id = H5Fcreate(new_fname, H5F_ACC_TRUNC, H5P_DEFAULT, access_props);
     else
	  id = H5Fopen(new_fname, H5F_ACC_RDWR, access_props);
#  endif

     CHECK(id >= 0, "error creating HDF output file");

     free(new_fname);

     H5Pclose(access_props);

     return id;
#else
     mpi_one_fprintf(stderr,
		     "matrixio: cannot output \"%s\" (compiled without HDF)\n",
		     fname);
     return 0;
#endif
}

matrixio_id matrixio_open(const char *fname, int read_only)
{
#if defined(HAVE_HDF5)
     char *new_fname;
     matrixio_id id;
     hid_t access_props;

     access_props = H5Pcreate (H5P_FILE_ACCESS);
     
#  if defined(HAVE_MPI) && defined(HAVE_H5PSET_FAPL_MPIO)
     H5Pset_fapl_mpio(access_props, MPI_COMM_WORLD, MPI_INFO_NULL);
#  endif

     new_fname = add_fname_suffix(fname);

     IF_EXCLUSIVE(mpi_begin_critical_section(matrixio_critical_section_tag),0);

     if (read_only)
	  id = H5Fopen(new_fname, H5F_ACC_RDONLY, access_props);
     else
	  id = H5Fopen(new_fname, H5F_ACC_RDWR, access_props);
     CHECK(id >= 0, "error opening HDF input file");

     free(new_fname);

     H5Pclose(access_props);

     return id;
#else
     CHECK(0, "no matrixio implementation is linked");
     return 0;
#endif
}

void matrixio_close(matrixio_id id)
{
#if defined(HAVE_HDF5)
     CHECK(H5Fclose(id) >= 0, "error closing HDF file");
     IF_EXCLUSIVE(mpi_end_critical_section(matrixio_critical_section_tag++),0);
#endif
}

/*****************************************************************************/

matrixio_id matrixio_create_sub(matrixio_id id, 
				const char *name, const char *description)
{
#if defined(HAVE_HDF5)
     matrixio_id sub_id;

     /* when running a parallel job, only the master process creates the
	group.  It flushes the group to disk and then the other processes
	open the group.  Is this the right thing to do, or is the
        H5Gcreate function parallel-aware? */

     if (mpi_is_master()) {
	  sub_id = H5Gcreate(id, name, 0 /* ==> default size */ );
	  matrixio_write_string_attr(sub_id, "description", description);
	  
	  H5Fflush(sub_id, H5F_SCOPE_GLOBAL);

	  IF_EXCLUSIVE(0,MPI_Barrier(MPI_COMM_WORLD));
     }
     else {
	  IF_EXCLUSIVE(0,MPI_Barrier(MPI_COMM_WORLD));

	  sub_id = H5Gopen(id, name);
     }

     return sub_id;