h5file.cpp

麻省理工的计算光子晶体的程序

/* Copyright (C) 2006 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
 */

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

#include "meep.hpp"

#define CHECK(condition, message) do { \
     if (!(condition))  { \
          abort("error on line %d of " __FILE__ ": " \
                  message "\n", __LINE__); \
     } \
} while (0)

#include "config.h"

#ifdef HAVE_HDF5
#  include <hdf5.h>

/* HDF5 changed this datatype in their interfaces starting in version 1.6.4 */
#  if H5_VERS_MAJOR > 1 \
   || (H5_VERS_MAJOR == 1 && H5_VERS_MINOR > 6) \
   || (H5_VERS_MAJOR == 1 && H5_VERS_MINOR == 6 && H5_VERS_RELEASE > 3)
     typedef hsize_t start_t;
#  else
     typedef hssize_t start_t;
#  endif

#else
typedef int hid_t;
#endif


#define HID(x) (*((hid_t *) (x)))

/*****************************************************************************/
/* 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

static int h5io_critical_section_tag = 0;

/*****************************************************************************/
/* 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); \
}

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

#ifdef HAVE_HDF5
static bool dataset_exists(hid_t id, const char *name)
{
     hid_t data_id;
     SUPPRESS_HDF5_ERRORS(data_id = H5Dopen(id, name));
     if (data_id >= 0)
          H5Dclose(data_id);
     return (data_id >= 0);
}
#endif

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

namespace meep {

// lazy file creation & locking
void *h5file::get_id() {
  if (HID(id) < 0) {
    if (parallel)
      all_wait();

#ifdef HAVE_HDF5
    hid_t access_props = H5Pcreate (H5P_FILE_ACCESS);
#  if defined(HAVE_MPI) && defined(HAVE_H5PSET_FAPL_MPIO)
    if (parallel)
      H5Pset_fapl_mpio(access_props, MPI_COMM_WORLD, MPI_INFO_NULL);
#  else
    if (parallel)
      begin_critical_section(h5io_critical_section_tag);
#  endif

    if (mode != WRITE || IF_EXCLUSIVE(parallel && !am_master(), 0))
      HID(id) = H5Fopen(filename,
			mode == READONLY ? H5F_ACC_RDONLY : H5F_ACC_RDWR,
			access_props);
    else
      HID(id) = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, access_props);
    
    H5Pclose(access_props);
#endif
  }
  return id;
}

// hackery: in some circumstances, for the exclusive-access mode
// we must close the id (i.e. the file) in order to prevent deadlock.
void h5file::prevent_deadlock() {
  IF_EXCLUSIVE(if (parallel) close_id(), 0);
}

void h5file::close_id() {
  unset_cur();
#ifdef HAVE_HDF5
  if (HID(id) >= 0)
    H5Fclose(HID(id));
  HID(id) = -1;
#endif
  IF_EXCLUSIVE(if (parallel)
	       end_critical_section(h5io_critical_section_tag++),
	       (void) 0);
  if (mode == WRITE) mode = READWRITE; // don't re-create on re-open
}

/* note: if parallel is true, then *all* processes must call this,
   and all processes will use I/O. */
h5file::h5file(const char *filename_, access_mode m, bool parallel_) {
  cur_dataname = NULL;
  id = (void*) (new hid_t);
  cur_id = (void*) (new hid_t);
  HID(id) = -1;
  HID(cur_id) = -1;
  extending = 0;
  filename = new char[strlen(filename_) + 1];
  strcpy(filename, filename_);
  mode = m;
  parallel = parallel_;
}

h5file::~h5file() {
  close_id();
  if (cur_dataname) free(cur_dataname); // allocated with realloc
  for (h5file::extending_s *cur = extending; cur; ) {
    h5file::extending_s *next = cur->next;
    delete cur->dataname;
    delete cur;
    cur = next;
  }
  delete[] filename;
}

bool h5file::ok() {
  return (HID(get_id()) >= 0);
}

void h5file::remove() {
  close_id();
  if (mode == READWRITE) mode = WRITE; // now need to re-create file
  for (h5file::extending_s *cur = extending; cur; ) {
    h5file::extending_s *next = cur->next;
    delete cur->dataname;
    delete cur;
    cur = next;
  }
  extending = 0;
  
  IF_EXCLUSIVE(if (parallel) all_wait(), 0);
  if (am_master() && std::remove(filename))
    abort("error removing file %s", filename);
}

h5file::extending_s *h5file::get_extending(const char *dataname) const {
  for (extending_s *cur = extending; cur; cur = cur->next)
    if (!strcmp(dataname, cur->dataname))
      return cur;
  return NULL;
}

bool h5file::is_cur(const char *dataname) {
  return cur_dataname && !strcmp(cur_dataname, dataname);
}

void h5file::unset_cur() {
#ifdef HAVE_HDF5
  if (HID(cur_id) >= 0)
    H5Dclose(HID(cur_id));
#endif
  HID(cur_id) = -1;
  if (cur_dataname)
    cur_dataname[0] = 0;
}

void h5file::set_cur(const char *dataname, void *data_id) {
#ifdef HAVE_HDF5
  if (HID(cur_id) >= 0 && HID(cur_id) != HID(data_id))
    H5Dclose(HID(cur_id));
#endif
  HID(cur_id) = HID(data_id);
  if (!is_cur(dataname)) {
    if (!cur_dataname || strlen(dataname) < strlen(cur_dataname))
      cur_dataname = (char *) realloc(cur_dataname, strlen(dataname) + 1);
    strcpy(cur_dataname, dataname);
  }
}

void h5file::read_size(const char *dataname, int *rank, int *dims, int maxrank)
{
#ifdef HAVE_HDF5
  if (parallel || am_master()) {
    hid_t file_id = HID(get_id()), space_id, data_id;
    
    CHECK(file_id >= 0, "error opening HDF5 input file");
    
    if (is_cur(dataname))
      data_id = HID(cur_id);
    else {
      CHECK(dataset_exists(file_id, dataname),
	    "missing dataset in HDF5 file");
      data_id = H5Dopen(file_id, dataname);
      set_cur(dataname, &data_id);
    }
    space_id = H5Dget_space(data_id);
    
    *rank = H5Sget_simple_extent_ndims(space_id);
    CHECK(*rank <= maxrank, "input array rank is too big");
    
    hsize_t *dims_copy = new hsize_t[*rank];
    hsize_t *maxdims = new hsize_t[*rank];
    H5Sget_simple_extent_dims(space_id, dims_copy, maxdims);
    for (int i = 0; i < *rank; ++i) dims[i] = dims_copy[i];
    delete[] maxdims;
    delete[] dims_copy;
    H5Sclose(space_id);
  }

  if (!parallel) {
    *rank = broadcast(0, *rank);
    broadcast(0, dims, *rank);
    
    if (*rank == 1 && dims[0] == 1)
      *rank = 0;
  }
#endif
}

double *h5file::read(const char *dataname,
		     int *rank, int *dims, int maxrank)
{
#ifdef HAVE_HDF5
  double *data = 0;
  if (parallel || am_master()) {
    int i, N;
    hid_t file_id = HID(get_id()), space_id, data_id;
    
    CHECK(file_id >= 0, "error opening HDF5 input file");

    if (is_cur(dataname))
      data_id = HID(cur_id);
    else {
      CHECK(dataset_exists(file_id, dataname),
	    "missing dataset in HDF5 file");
      data_id = H5Dopen(file_id, dataname);
    }
    space_id = H5Dget_space(data_id);
    
    *rank = H5Sget_simple_extent_ndims(space_id);
    CHECK(*rank <= maxrank, "input array rank is too big");
    
    hsize_t *dims_copy = new hsize_t[*rank];
    hsize_t *maxdims = new hsize_t[*rank];
    H5Sget_simple_extent_dims(space_id, dims_copy, maxdims);
    delete[] maxdims;
    for (N = 1, i = 0; i < *rank; ++i)
      N *= (dims[i] = dims_copy[i]);
    delete[] dims_copy;
    H5Sclose(space_id);
    
    data = new double[N];
    H5Dread(data_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
	    (void *) data);
    
    if (!is_cur(dataname))
      H5Dclose(data_id);
  }
  
  if (!parallel) {
    *rank = broadcast(0, *rank);
    broadcast(0, dims, *rank);
    int N = 1;
    for (int i = 0; i < *rank; ++i)
      N *= dims[i];
    if (!am_master())
      data = new double[N];
    broadcast(0, data, N);
  }

  if (*rank == 1 && dims[0] == 1)
    *rank = 0;
  
  return data;
#else
  return NULL;
#endif
}

char *h5file::read(const char *dataname)
{
#ifdef HAVE_HDF5
  char *data = 0;
  int len = 0;
  if (parallel || am_master()) {
    hid_t file_id = HID(get_id()), space_id, data_id, type_id;
    
    CHECK(file_id >= 0, "error opening HDF5 input file");
    
    if (is_cur(dataname))
      unset_cur();

    CHECK(dataset_exists(file_id, dataname),
	  "missing dataset in HDF5 file");
    
    data_id = H5Dopen(file_id, dataname);
    space_id = H5Dget_space(data_id);
    type_id = H5Dget_type(data_id);
    
    CHECK(H5Sget_simple_extent_npoints(space_id) == 1,
	  "expected single string in HDF5 file, but didn't get one");
    
    len = H5Tget_size(type_id);
    H5Tclose(type_id);
    type_id = H5Tcopy(H5T_C_S1);
    H5Tset_size(type_id, len);
    
    data = new char[len];
    H5Dread(data_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT,
	    (void *) data);
    
    H5Tclose(type_id);
    H5Sclose(space_id);	  
    H5Dclose(data_id);
  }
  
  if (!parallel) {
    len = broadcast(0, len);
    if (!am_master())
      data = new char[len];
    broadcast(0, data, len);
  }
  
  return data;
#else
  return NULL;
#endif
}

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

/* Delete a dataset, if it exists.  In parallel mode, should be called
   by all processors. */
void h5file::remove_data(const char *dataname)
{
#ifdef HAVE_HDF5
  hid_t file_id = HID(get_id());

  if (is_cur(dataname))
    unset_cur();

  if (get_extending(dataname)) { // delete dataname from extending list
    extending_s *prev = 0, *cur = extending;