--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/ls-hdf5.cc	Wed Nov 19 21:23:19 2003 +0000
@@ -0,0 +1,1241 @@
+/*
+
+Copyright (C) 1996, 1997 John W. Eaton
+
+This file is part of Octave.
+
+Octave 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, or (at your option) any
+later version.
+
+Octave 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 Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+// Author: Steven G. Johnson <stevenj@alum.mit.edu>
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#if defined (HAVE_HDF5)
+
+#include <cfloat>
+#include <cstring>
+#include <cctype>
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <memory>
+#include <string>
+
+#include <hdf5.h>
+
+#include "byte-swap.h"
+#include "data-conv.h"
+#include "file-ops.h"
+#include "glob-match.h"
+#include "lo-mappers.h"
+#include "lo-sstream.h"
+#include "mach-info.h"
+#include "oct-env.h"
+#include "oct-time.h"
+#include "quit.h"
+#include "str-vec.h"
+
+#include "Cell.h"
+#include "defun.h"
+#include "error.h"
+#include "gripes.h"
+#include "load-save.h"
+#include "oct-obj.h"
+#include "oct-map.h"
+#include "ov-cell.h"
+#include "pager.h"
+#include "pt-exp.h"
+#include "symtab.h"
+#include "sysdep.h"
+#include "unwind-prot.h"
+#include "utils.h"
+#include "variables.h"
+#include "version.h"
+#include "dMatrix.h"
+
+#include "ls-utils.h"
+#include "ls-hdf5.h"
+
+// this is only used for HDF5 import
+// try to convert s into a valid identifier, replacing invalid chars with "_":
+
+static std::string
+make_valid_identifier (const std::string& nm)
+{
+  std::string retval;
+
+  size_t nm_len = nm.length ();
+
+  if (nm_len > 0)
+    {
+      if (! isalpha (nm[0]))
+	retval += '_';
+
+      for (size_t i = 0; i < nm_len; i++)
+	{
+	  char c = nm[i];
+	  retval += (isalnum (c) || c == '_') ? c : '_';
+	}
+    }
+
+  return retval;
+}
+
+static bool
+ident_is_all_digits (const std::string& id)
+{
+  bool retval = true;
+
+  size_t len = 0;
+
+  for (size_t i = 0; i < len; i++)
+    {
+      if (! isdigit (id[i]))
+	{
+	  retval = false;
+	  break;
+	}
+    }
+
+  return retval;
+}
+
+// Define this to 1 if/when HDF5 supports automatic conversion between
+// integer and floating-point binary data:
+#define HAVE_HDF5_INT2FLOAT_CONVERSIONS 0
+
+// Given two compound types t1 and t2, determine whether they 
+// are compatible for reading/writing.  This function only
+// works for non-nested types composed of simple elements (ints, floats...),
+// which is all we need it for
+
+bool
+hdf5_types_compatible (hid_t t1, hid_t t2)
+{
+  int n;
+  if ((n = H5Tget_nmembers (t1)) != H5Tget_nmembers (t2))
+    return false;
+
+  for (int i = 0; i < n; ++i)
+    {
+      hid_t mt1 = H5Tget_member_type (t1, i);
+      hid_t mt2 = H5Tget_member_type (t2, i);
+
+      if (H5Tget_class (mt1) != H5Tget_class (mt2))
+	return false;
+
+      H5Tclose (mt2);
+      H5Tclose (mt1);
+    }
+
+  return true;
+}
+
+// Import a multidimensional (rank >= 3) dataset whose id is data_id, into tc.
+// This works by calling itself recursively, building up lists of lists
+//  of lists ... of 2d matrices.  rank and dims are the rank and dimensions
+//  of the dataset.  type_id is the datatype to read into.  If it is
+//  H5T_NATIVE_DOUBLE, we are reading a real matrix.  Otherwise, type_id
+//  is assumed to be a complex type for reading a complex matrix.
+//
+//  Upon entry, we should have curdim = rank - 1, start = an array
+//  of length rank = all zeros, and count = an array of length rank =
+//  all ones except for the first two dimensions which equal the corresponding
+//  entries in dims[]. 
+//
+//  Note that we process the dimensions in reverse order, reflecting
+//  the fact that Octave is uses column-major (Fortran-order) data while
+//  HDF5 is row-major.  This means that the HDF5 file is read
+//  non-contiguously, but on the other hand means that for a 3d array
+//  we get a list of xy-plane slices, which seems nice.  We could change
+//  this behavior without much trouble; what is the best thing to do?
+//
+//  Returns a positive value upon success.
+
+static herr_t
+hdf5_import_multidim (hid_t data_id, hid_t space_id, hsize_t rank,
+		      const hsize_t *dims, hsize_t curdim,
+		      hssize_t *start, const hsize_t *count,
+		      hid_t type_id, octave_value &tc)
+{
+  herr_t retval = 1;
+
+  if (rank < 3 || curdim < 1 || curdim >= rank)
+    return -1;
+
+  if (curdim == 1)
+    {
+      // import 2d dataset for 1st 2 dims directly as a matrix
+      int nr, nc;    // rows and columns
+      nc = dims[0];  // octave uses column-major & HDF5 uses row-major
+      nr = dims[1];
+
+      hid_t mem_space_id = H5Screate_simple (2, dims, 0);
+
+      if (mem_space_id < 0)
+	return -1;
+
+      if (H5Sselect_all (mem_space_id) < 0)
+	return -1;
+    
+      if (H5Sselect_hyperslab (space_id, H5S_SELECT_SET,
+			       start, 0, count, 0) < 0)
+	{
+	  H5Sclose (mem_space_id);
+	  return -1;
+	}
+    
+      if (type_id == H5T_NATIVE_DOUBLE)
+	{
+	  // real matrix
+	  Matrix m (nr, nc);
+	  double *re = m.fortran_vec ();
+	  if (H5Dread (data_id, type_id, mem_space_id, space_id,
+		       H5P_DEFAULT, (void *) re) < 0)
+	    retval = -1;  // error
+	  else
+	    tc = m;
+	}
+      else
+	{
+	  // assume that we are using complex numbers
+	  // complex matrix
+	  ComplexMatrix m (nr, nc);
+	  Complex *reim = m.fortran_vec ();
+	  if (H5Dread (data_id, type_id, mem_space_id, space_id,
+		       H5P_DEFAULT, (void *) X_CAST (double *, reim)) < 0)
+	    retval = -1;  // error
+	  else
+	    tc = m;
+	}
+    
+      H5Sclose (mem_space_id);
+
+    }
+  else
+    {
+      octave_value_list lst;
+
+      for (hsize_t i = 0; i < dims[curdim]; ++i)
+	{
+	  octave_value slice;
+	  start[curdim] = i;
+	  retval = hdf5_import_multidim (data_id, space_id, rank,
+					 dims, curdim-1, start, count,
+					 type_id, slice);
+	  if (retval < 0)
+	    break;
+	  lst.append (slice);
+	}
+
+      if (retval > 0)
+	tc = octave_value (lst);
+    }
+
+  return retval;
+}
+
+// Return true if loc_id has the attribute named attr_name, and false
+// otherwise.
+
+bool
+hdf5_check_attr (hid_t loc_id, const char *attr_name)
+{
+  bool retval = false;
+
+  // we have to pull some shenanigans here to make sure
+  // HDF5 doesn't print out all sorts of error messages if we
+  // call H5Aopen for a non-existing attribute
+
+  H5E_auto_t err_func;
+  void *err_func_data;
+
+  // turn off error reporting temporarily, but save the error
+  // reporting function:
+
+  H5Eget_auto (&err_func, &err_func_data);
+  H5Eset_auto (0, 0);
+
+  hid_t attr_id = H5Aopen_name (loc_id, attr_name);
+
+  if (attr_id >= 0)
+    {
+      // successful
+      retval = 1;
+      H5Aclose (attr_id);
+    }
+
+  // restore error reporting:
+  H5Eset_auto (err_func, err_func_data);
+
+  return retval;
+}
+
+// The following two subroutines create HDF5 representations of the way
+// we will store Octave complex and range types (pairs and triplets of
+// floating-point numbers, respectively).  NUM_TYPE is the HDF5 numeric
+// type to use for storage (e.g. H5T_NATIVE_DOUBLE to save as 'double').
+// Note that any necessary conversions are handled automatically by HDF5.
+
+static hid_t
+hdf5_make_complex_type (hid_t num_type)
+{
+  hid_t type_id = H5Tcreate (H5T_COMPOUND, sizeof (double) * 2);
+
+  H5Tinsert (type_id, "real", 0 * sizeof (double), num_type);
+  H5Tinsert (type_id, "imag", 1 * sizeof (double), num_type);
+
+  return type_id;
+}
+
+static hid_t
+hdf5_make_range_type (hid_t num_type)
+{
+  hid_t type_id = H5Tcreate (H5T_COMPOUND, sizeof (double) * 3);
+
+  H5Tinsert (type_id, "base", 0 * sizeof (double), num_type);
+  H5Tinsert (type_id, "limit", 1 * sizeof (double), num_type);
+  H5Tinsert (type_id, "increment", 2 * sizeof (double), num_type);
+
+  return type_id;
+}
+
+// Callback data structure for passing data to hdf5_read_next_data, below.
+
+struct
+hdf5_callback_data
+{
+  hdf5_callback_data (void)
+    : name (), global (false), tc (), doc (),
+      complex_type (hdf5_make_complex_type (H5T_NATIVE_DOUBLE)),
+      range_type (hdf5_make_range_type (H5T_NATIVE_DOUBLE)),
+      import (false) { }
+
+  // the following fields are set by hdf5_read_data on successful return:
+
+  // the name of the variable
+  std::string name;
+
+  // whether it is global
+  bool global;
+
+  // the value of the variable, in Octave form
+  octave_value tc;
+
+  // a documentation string (NULL if none)
+  std::string doc;
+
+  // the following fields are input to hdf5_read_data:
+
+  // HDF5 rep's of complex and range type
+  hid_t complex_type, range_type;
+
+  // whether to try extra hard to import "foreign" data
+  bool import;
+};
+
+// This variable, set in read_hdf5_data(), tells whether we are using
+// a version of HDF5 with a buggy H5Giterate (i.e. which neglects to
+// increment the index parameter to the next unread item).
+static bool have_h5giterate_bug = false;
+
+// This function is designed to be passed to H5Giterate, which calls it
+// on each data item in an HDF5 file.  For the item whose name is NAME in
+// the group GROUP_ID, this function sets dv->tc to an Octave representation
+// of that item.  (dv must be a pointer to hdf5_callback_data.)  (It also
+// sets the other fields of dv).
+//
+// It returns 1 on success (in which case H5Giterate stops and returns),
+// -1 on error, and 0 to tell H5Giterate to continue on to the next item
+// (e.g. if NAME was a data type we don't recognize).
+
+static herr_t
+hdf5_read_next_data (hid_t group_id, const char *name, void *dv)
+{
+  hdf5_callback_data *d = static_cast <hdf5_callback_data *> (dv);
+
+  H5G_stat_t info;
+  herr_t retval = 0;
+  bool ident_valid = valid_identifier (name);
+
+  std::string vname = name;
+
+  // Allow identifiers as all digits so we can load lists saved by
+  // earlier versions of Octave.
+
+  if (! ident_valid && (d->import || ident_is_all_digits (vname)))
+    {
+      // fix the identifier, replacing invalid chars with underscores
+      vname = make_valid_identifier (vname);
+
+      // check again (in case vname was null, empty, or some such thing):
+      ident_valid = valid_identifier (vname); 
+    }
+
+  H5Gget_objinfo (group_id, name, 1, &info);
+
+  if (info.type == H5G_DATASET && ident_valid)
+    {
+      retval = 1;
+
+      hid_t data_id = H5Dopen (group_id, name);
+
+      if (data_id < 0)
+	{
+	  retval = data_id;
+
+	  goto done;
+	}
+
+      hid_t type_id = H5Dget_type (data_id);
+
+      hid_t type_class_id = H5Tget_class (type_id);
+
+#if HAVE_HDF5_INT2FLOAT_CONVERSIONS
+      if (type_class_id == H5T_INTEGER || type_class_id == H5T_FLOAT)
+	{
+#else
+      // hdf5 doesn't (yet) support automatic float/integer conversions
+      if (type_class_id == H5T_FLOAT)
+	{
+#endif
+	  // read real matrix or scalar variable
+
+	  hid_t space_id = H5Dget_space (data_id);
+
+	  hsize_t rank = H5Sget_simple_extent_ndims (space_id);
+
+	  if (rank == 0)
+	    {
+	      // real scalar:
+	      double dtmp;
+	      if (H5Dread (data_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, 
+			   H5P_DEFAULT, (void *) &dtmp) < 0)
+		retval = -1;  // error
+	      else
+		d->tc = dtmp;
+	    }
+	  else if (rank > 0 && rank <= 2)
+	    {
+	      // real matrix
+	      OCTAVE_LOCAL_BUFFER (hsize_t, dims, rank);
+	      OCTAVE_LOCAL_BUFFER (hsize_t, maxdims, rank);
+
+	      H5Sget_simple_extent_dims (space_id, dims, maxdims);
+
+	      int nr, nc;  // rows and columns
+	      // octave uses column-major & HDF5 uses row-major
+	      nc = dims[0];
+	      nr = rank > 1 ? dims[1] : 1;
+	      Matrix m (nr, nc);
+	      double *re = m.fortran_vec ();
+	      if (H5Dread (data_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, 
+			   H5P_DEFAULT, (void *) re) < 0)
+		retval = -1;  // error
+	      else
+		d->tc = m;
+	    }
+	  else if (rank >= 3 && d->import)
+	    {
+	      OCTAVE_LOCAL_BUFFER (hsize_t, dims, rank);
+	      OCTAVE_LOCAL_BUFFER (hsize_t, maxdims, rank);
+
+	      H5Sget_simple_extent_dims (space_id, dims, maxdims);
+
+	      OCTAVE_LOCAL_BUFFER (hssize_t, start, rank);
+	      OCTAVE_LOCAL_BUFFER (hsize_t, count, rank);
+
+	      for (hsize_t i = 0; i < rank; ++i)
+		{
+		  start[i] = 0;
+		  count[i] = 1;
+		}
+	      count[0] = dims[0];
+	      count[1] = dims[1];
+	      retval = hdf5_import_multidim (data_id, space_id,
+					     rank, dims, rank-1,
+					     start, count,
+					     H5T_NATIVE_DOUBLE, d->tc);
+	    }
+	  else
+	    {
+	      warning ("load: can't read %d-dim. hdf5 dataset %s",
+		       rank, name);
+	      retval = 0;  // skip; we can't read 3+ dimensional datasets
+	    }
+
+	  H5Sclose (space_id);
+	}
+      else if (type_class_id == H5T_STRING)
+	{
+	  // read string variable
+	  hid_t space_id = H5Dget_space (data_id);
+	  hsize_t rank = H5Sget_simple_extent_ndims (space_id);
+
+	  if (rank == 0)
+	    {
+	      // a single string:
+	      int slen = H5Tget_size (type_id);
+	      if (slen < 0)
+		retval = -1;  // error
+	      else
+		{
+		  OCTAVE_LOCAL_BUFFER (char, s, slen);
+		  // create datatype for (null-terminated) string
+		  // to read into:
+		  hid_t st_id = H5Tcopy (H5T_C_S1);
+		  H5Tset_size (st_id, slen);
+		  if (H5Dread (data_id, st_id, H5S_ALL, H5S_ALL, 
+			       H5P_DEFAULT, (void *) s) < 0)
+		    {
+		      retval = -1;  // error
+		    }
+		  else
+		    d->tc = s;
+
+		  H5Tclose (st_id);
+		}
+	    }
+	  else if (rank == 1)
+	    {
+	      // string vector
+	      hsize_t elements, maxdim;
+	      H5Sget_simple_extent_dims (space_id, &elements, &maxdim);
+	      int slen = H5Tget_size (type_id);
+	      if (slen < 0)
+		retval = -1;  // error
+	      else
+		{
+		  // hdf5 string arrays store strings of all the
+		  // same physical length (I think), which is
+		  // slightly wasteful, but oh well.
+
+		  OCTAVE_LOCAL_BUFFER (char, s, elements * slen);
+
+		  // create datatype for (null-terminated) string
+		  // to read into:
+		  hid_t st_id = H5Tcopy (H5T_C_S1);
+		  H5Tset_size (st_id, slen);
+
+		  if (H5Dread (data_id, st_id, H5S_ALL, H5S_ALL, 
+			       H5P_DEFAULT, (void *) s) < 0)
+		    retval = -1;  // error
+		  else
+		    {
+		      charMatrix chm (elements, slen - 1);
+		      for (hsize_t i = 0; i < elements; ++i)
+			{
+			  chm.insert (s + i*slen, i, 0);
+			}
+		      d->tc = octave_value (chm, true);
+		    }
+
+		  H5Tclose (st_id);
+		}
+	    }
+	  else
+	    {
+	      warning ("load: can't read %d-dim. hdf5 string vector %s",
+		       rank, name); 
+	      // skip; we can't read higher-dimensional string vectors
+	      retval = 0;
+	    }
+	}
+      else if (type_class_id == H5T_COMPOUND)
+	{
+	  // check for complex or range data:
+
+	  if (hdf5_types_compatible (type_id, d->complex_type))
+	    {
+	      // read complex matrix or scalar variable
+
+	      hid_t space_id = H5Dget_space (data_id);
+	      hsize_t rank = H5Sget_simple_extent_ndims (space_id);
+
+	      if (rank == 0)
+		{
+		  // complex scalar:
+		  Complex ctmp;
+		  if (H5Dread (data_id, d->complex_type, H5S_ALL,
+			       H5S_ALL, H5P_DEFAULT,
+			       (void *) X_CAST (double *, &ctmp)) < 0)
+		    retval = -1;  // error
+		  else
+		    d->tc = ctmp;
+		}
+	      else if (rank > 0 && rank <= 2)
+		{
+		  // complex matrix
+		  OCTAVE_LOCAL_BUFFER (hsize_t, dims, rank);
+		  OCTAVE_LOCAL_BUFFER (hsize_t, maxdims, rank);
+		  H5Sget_simple_extent_dims (space_id, dims, maxdims);
+		  int nr, nc;  // rows and columns
+		  // octave uses column-major & HDF5 uses row-major
+		  nc = dims[0];
+		  nr = rank > 1 ? dims[1] : 1;
+		  ComplexMatrix m (nr, nc);
+		  Complex *reim = m.fortran_vec ();
+		  if (H5Dread (data_id, d->complex_type, H5S_ALL,
+			       H5S_ALL, H5P_DEFAULT,
+			       (void *) X_CAST (double *, reim)) < 0)
+		    retval = -1;  // error
+		  else
+		    d->tc = m;
+		}
+	      else if (rank >= 3 && d->import)
+		{
+		  OCTAVE_LOCAL_BUFFER (hsize_t, dims, rank);
+		  OCTAVE_LOCAL_BUFFER (hsize_t, maxdims, rank);
+		  H5Sget_simple_extent_dims (space_id, dims, maxdims);
+		  OCTAVE_LOCAL_BUFFER (hssize_t, start, rank);
+		  OCTAVE_LOCAL_BUFFER (hsize_t, count, rank);
+		  for (hsize_t i = 0; i < rank; ++i)
+		    {
+		      start[i] = 0;
+		      count[i] = 1;
+		    }
+		  count[0] = dims[0];
+		  count[1] = dims[1];
+		  retval = hdf5_import_multidim (data_id, space_id,
+						 rank, dims, rank-1,
+						 start, count,
+						 d->complex_type,
+						 d->tc);
+		}
+	      else
+		{
+		  warning ("load: can't read %d-dim. hdf5 dataset %s",
+			   rank, name);
+		  // skip; we can't read 3+ dimensional datasets
+		  retval = 0;
+		}
+	      H5Sclose (space_id);
+	    }
+	  else if (hdf5_types_compatible (type_id, d->range_type))
+	    {
+	      // read range variable:
+	      hid_t space_id = H5Dget_space (data_id);
+	      hsize_t rank = H5Sget_simple_extent_ndims (space_id);
+
+	      if (rank == 0)
+		{
+		  double rangevals[3];
+		  if (H5Dread (data_id, d->range_type, H5S_ALL, H5S_ALL, 
+			       H5P_DEFAULT, (void *) rangevals) < 0)
+		    retval = -1;  // error
+		  else
+		    {
+		      Range r (rangevals[0], rangevals[1], rangevals[2]);
+		      d->tc = r;
+		    }
+		}
+	      else
+		{
+		  warning ("load: can't read range array `%s' in hdf5 file",
+			   name);
+		  // skip; we can't read arrays of range variables
+		  retval = 0;
+		}
+
+	      H5Sclose (space_id);
+	    }
+	  else
+	    {
+	      warning ("load: can't read `%s' (unknown compound datatype)",
+		       name);
+	      retval = 0; // unknown datatype; skip.
+	    }
+	}
+      else
+	{
+	  warning ("load: can't read `%s' (unknown datatype)", name);
+	  retval = 0; // unknown datatype; skip
+	}
+
+      H5Tclose (type_id);
+
+      // check for OCTAVE_GLOBAL attribute:
+      d->global = hdf5_check_attr (data_id, "OCTAVE_GLOBAL");
+
+      H5Dclose (data_id);
+    }
+  else if (info.type == H5G_GROUP && ident_valid)
+    {
+      // read in group as a list or a structure
+      retval = 1;
+
+      hid_t subgroup_id = H5Gopen (group_id, name);
+
+      if (subgroup_id < 0)
+	{
+	  retval = subgroup_id;
+	  goto done;
+	}
+
+      // an HDF5 group is treated as an octave structure by
+      // default (since that preserves name information), and an
+      // octave list otherwise.
+
+      bool is_list = hdf5_check_attr (subgroup_id, "OCTAVE_LIST");
+
+      hdf5_callback_data dsub;
+
+      dsub.complex_type = d->complex_type;
+      dsub.range_type = d->range_type;
+      dsub.import = d->import;
+
+      herr_t retval2;
+      octave_value_list lst;
+      Octave_map m;
+      int current_item = 0;
+      while ((retval2 = H5Giterate (group_id, name, &current_item,
+				    hdf5_read_next_data, &dsub)) > 0)
+	{
+	  if (is_list)
+	    lst.append (dsub.tc);
+	  else
+	    {
+	      octave_value ov = dsub.tc;
+
+	      if (ov.is_list ())
+		m [dsub.name] = ov.list_value ();
+	      else
+		m [dsub.name] = ov;
+	    }
+
+	  if (have_h5giterate_bug)
+	    current_item++;  // H5Giterate returned the last index processed
+	}
+
+      if (retval2 < 0)
+	retval = retval2;
+      else
+	{
+	  d->global = hdf5_check_attr (group_id, "OCTAVE_GLOBAL");
+
+	  if (is_list)
+	    d->tc = octave_value (lst);
+	  else
+	    d->tc = m;
+	}
+
+      H5Gclose (subgroup_id);
+    }
+  else if (! ident_valid)
+    {
+      // should we attempt to handle invalid identifiers by converting
+      // bad characters to '_', say?
+      warning ("load: skipping invalid identifier `%s' in hdf5 file",
+	       name);
+    }
+
+ done:
+
+  if (retval < 0)
+    error ("load: error while reading hdf5 item %s", name);
+
+  if (retval > 0)
+    {
+      // get documentation string, if any:
+      int comment_length = H5Gget_comment (group_id, name, 0, 0);
+
+      if (comment_length > 1)
+	{
+	  OCTAVE_LOCAL_BUFFER (char, tdoc, comment_length);
+	  H5Gget_comment (group_id, name, comment_length, tdoc);
+	  d->doc = tdoc;
+	}
+      else if (vname != name)
+	{
+	  // the name was changed by import; store the original name
+	  // as the documentation string:
+	  d->doc = name;
+	}
+
+      // copy name (actually, vname):
+      d->name = vname;
+    }
+
+  return retval;
+}
+
+// Read the next Octave variable from the stream IS, which must really be
+// an hdf5_ifstream.  Return the variable value in tc, its doc string
+// in doc, and whether it is global in global.  The return value is
+// the name of the variable, or NULL if none were found or there was
+// and error.  If import is true, we try extra hard to import "foreign"
+// datasets (not created by Octave), although we usually do a reasonable
+// job anyway.  (c.f. load -import documentation.)
+std::string
+read_hdf5_data (std::istream& is,
+		const std::string& filename, bool& global,
+		octave_value& tc, std::string& doc, bool import)
+{
+  std::string retval;
+
+  doc.resize (0);
+
+  hdf5_ifstream& hs = (hdf5_ifstream&) is;
+  hdf5_callback_data d;
+
+  d.import = import;
+
+  // Versions of HDF5 prior to 1.2.2 had a bug in H5Giterate where it
+  // would return the index of the last item processed instead of the
+  // next item to be processed, forcing us to increment the index manually.
+
+  unsigned int vers_major, vers_minor, vers_release;
+
+  H5get_libversion (&vers_major, &vers_minor, &vers_release);
+
+  // XXX FIXME XXX -- this test looks wrong.
+  have_h5giterate_bug
+    = (vers_major < 1
+       || (vers_major == 1 && (vers_minor < 2
+			       || (vers_minor == 2 && vers_release < 2))));
+
+  herr_t H5Giterate_retval = H5Giterate (hs.file_id, "/", &hs.current_item,
+					 hdf5_read_next_data, &d);
+
+  if (have_h5giterate_bug)
+    {
+      // H5Giterate sets current_item to the last item processed; we want
+      // the index of the next item (for the next call to read_hdf5_data)
+
+      hs.current_item++;
+    }
+
+  if (H5Giterate_retval > 0)
+    {
+      global = d.global;
+      tc = d.tc;
+      doc = d.doc;
+    }
+  else
+    {
+      // an error occurred (H5Giterate_retval < 0) or there are no
+      // more datasets print an error message if retval < 0?
+      // hdf5_read_next_data already printed one, probably.
+    }
+
+  H5Tclose (d.complex_type);
+  H5Tclose (d.range_type);
+
+  if (! d.name.empty ())
+    retval = d.name;
+
+  return retval;
+}
+
+// Add an attribute named attr_name to loc_id (a simple scalar
+// attribute with value 1).  Return value is >= 0 on success.
+static herr_t
+hdf5_add_attr (hid_t loc_id, const char *attr_name)
+{
+  herr_t retval = 0;
+
+  hid_t as_id = H5Screate (H5S_SCALAR);
+
+  if (as_id >= 0)
+    {
+      hid_t a_id = H5Acreate (loc_id, attr_name,
+			      H5T_NATIVE_UCHAR, as_id, H5P_DEFAULT);
+
+      if (a_id >= 0)
+	{
+	  unsigned char attr_val = 1;
+
+	  retval = H5Awrite (a_id, H5T_NATIVE_UCHAR, (void*) &attr_val);
+
+	  H5Aclose (a_id);
+	}
+      else
+	retval = a_id;
+
+      H5Sclose (as_id);
+    }
+  else
+    retval = as_id;
+
+  return retval;
+}
+
+
+// save_type_to_hdf5 is not currently used, since hdf5 doesn't yet support
+// automatic float<->integer conversions:
+
+#if HAVE_HDF5_INT2FLOAT_CONVERSIONS
+
+// return the HDF5 type id corresponding to the Octave save_type
+
+static hid_t
+save_type_to_hdf5 (save_type st)
+{
+  switch (st)
+    {
+    case LS_U_CHAR:
+      return H5T_NATIVE_UCHAR;
+
+    case LS_U_SHORT:
+      return H5T_NATIVE_USHORT;
+
+    case LS_U_INT:
+      return H5T_NATIVE_UINT;
+
+    case LS_CHAR:
+      return H5T_NATIVE_CHAR;
+
+    case LS_SHORT:
+      return H5T_NATIVE_SHORT;
+
+    case LS_INT:
+      return H5T_NATIVE_INT;
+
+    case LS_FLOAT:
+      return H5T_NATIVE_FLOAT;
+
+    case LS_DOUBLE:
+    default:
+      return H5T_NATIVE_DOUBLE;
+    }
+}
+#endif /* HAVE_HDF5_INT2FLOAT_CONVERSIONS */
+
+// Add the data from TC to the HDF5 location loc_id, which could
+// be either a file or a group within a file.  Return true if
+// successful.  This function calls itself recursively for lists
+// (stored as HDF5 groups).
+
+static bool
+add_hdf5_data (hid_t loc_id, const octave_value& tc,
+	       const std::string& name, const std::string& doc,
+	       bool mark_as_global, bool save_as_floats)
+{
+  hsize_t dims[3];
+  hid_t type_id = -1, space_id = -1, data_id = -1;
+  bool data_is_group = 0;
+  bool retval = 0;
+
+  if (tc.is_string ())
+    {
+      int nr = tc.rows ();
+      charMatrix chm = tc.char_matrix_value ();
+      int nc = chm.cols ();
+
+      // create datatype for (null-terminated) string to write from:
+      type_id = H5Tcopy (H5T_C_S1); H5Tset_size (type_id, nc + 1);
+      if (type_id < 0)
+	goto error_cleanup;
+
+      dims[0] = nr;
+      space_id = H5Screate_simple (nr > 0 ? 1 : 0, dims, (hsize_t*) 0);
+      if (space_id < 0)
+	goto error_cleanup;
+
+      data_id = H5Dcreate (loc_id, name.c_str (), 
+			   type_id, space_id, H5P_DEFAULT);
+      if (data_id < 0)
+	goto error_cleanup;
+
+      OCTAVE_LOCAL_BUFFER (char, s, nr * (nc + 1));
+
+      for (int i = 0; i < nr; ++i)
+	{
+	  std::string tstr = chm.row_as_string (i);
+	  strcpy (s + i * (nc+1), tstr.c_str ());
+	}
+
+      if (H5Dwrite (data_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+		    (void*) s) < 0) {
+	goto error_cleanup;
+      }
+    }
+  else if (tc.is_range ())
+    {
+      space_id = H5Screate_simple (0, dims, (hsize_t*) 0);
+      if (space_id < 0)
+	goto error_cleanup;
+
+      type_id = hdf5_make_range_type (H5T_NATIVE_DOUBLE);
+      if (type_id < 0)
+	goto error_cleanup;
+
+      data_id = H5Dcreate (loc_id, name.c_str (), 
+			   type_id, space_id, H5P_DEFAULT);
+      if (data_id < 0)
+	goto error_cleanup;
+    
+      Range r = tc.range_value ();
+      double range_vals[3];
+      range_vals[0] = r.base ();
+      range_vals[1] = r.limit ();
+      range_vals[2] = r.inc ();
+
+      if (H5Dwrite (data_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+		    (void*) range_vals) < 0)
+	goto error_cleanup;
+    }
+  else if (tc.is_real_scalar ())
+    {
+      space_id = H5Screate_simple (0, dims, (hsize_t*) 0);
+      if (space_id < 0) goto error_cleanup;
+
+      data_id = H5Dcreate (loc_id, name.c_str (), 
+			   H5T_NATIVE_DOUBLE, space_id, H5P_DEFAULT);
+      if (data_id < 0)
+	goto error_cleanup;
+
+      double tmp = tc.double_value ();
+      if (H5Dwrite (data_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+		    H5P_DEFAULT, (void*) &tmp) < 0)
+	goto error_cleanup;
+    }
+  else if (tc.is_real_matrix ())
+    {
+      Matrix m = tc.matrix_value ();
+      dims[1] = m.rows ();    // Octave uses column-major, while
+      dims[0] = m.columns (); // HDF5 uses row-major ordering
+
+      space_id = H5Screate_simple (dims[1] > 1 ?2:1, dims, (hsize_t*) 0);
+      if (space_id < 0)
+	goto error_cleanup;
+
+      hid_t save_type_id = H5T_NATIVE_DOUBLE;
+
+      if (save_as_floats)
+	{
+	  if (m.too_large_for_float ())
+	    {
+	      warning ("save: some values too large to save as floats --");
+	      warning ("save: saving as doubles instead");
+	    }
+	  else
+	    save_type_id = H5T_NATIVE_FLOAT;
+	}
+#if HAVE_HDF5_INT2FLOAT_CONVERSIONS
+      // hdf5 currently doesn't support float/integer conversions
+      else
+	{
+	  double max_val, min_val;
+
+	  if (m.all_integers (max_val, min_val))
+	    save_type_id
+	      = save_type_to_hdf5 (get_save_type (max_val, min_val));
+	}
+#endif /* HAVE_HDF5_INT2FLOAT_CONVERSIONS */
+
+      data_id = H5Dcreate (loc_id, name.c_str (), 
+			   save_type_id, space_id, H5P_DEFAULT);
+      if (data_id < 0)
+	goto error_cleanup;
+    
+      double *mtmp = m.fortran_vec ();
+      if (H5Dwrite (data_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+		    H5P_DEFAULT, (void*) mtmp) < 0)
+	goto error_cleanup;
+    }
+  else if (tc.is_complex_scalar ())
+    {
+      space_id = H5Screate_simple (0, dims, (hsize_t*) 0);
+      if (space_id < 0)
+	goto error_cleanup;
+
+      type_id = hdf5_make_complex_type (H5T_NATIVE_DOUBLE);
+      if (type_id < 0)
+	goto error_cleanup;
+
+      data_id = H5Dcreate (loc_id, name.c_str (), 
+			   type_id, space_id, H5P_DEFAULT);
+      if (data_id < 0)
+	goto error_cleanup;
+
+      Complex tmp = tc.complex_value ();
+      if (H5Dwrite (data_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+		    (void*) X_CAST (double*, &tmp)) < 0)
+	goto error_cleanup;
+    }
+  else if (tc.is_complex_matrix ())
+    {
+      ComplexMatrix m = tc.complex_matrix_value ();
+
+      dims[1] = m.rows ();    // Octave uses column-major, while
+      dims[0] = m.columns (); // HDF5 uses row-major ordering
+
+      space_id = H5Screate_simple (dims[1] > 1 ?2:1, dims, (hsize_t*) 0);
+      if (space_id < 0)
+	goto error_cleanup;
+
+      hid_t save_type_id = H5T_NATIVE_DOUBLE;
+
+      if (save_as_floats)
+	{
+	  if (m.too_large_for_float ())
+	    {
+	      warning ("save: some values too large to save as floats --");
+	      warning ("save: saving as doubles instead");
+	    }
+	  else
+	    save_type_id = H5T_NATIVE_FLOAT;
+	}
+#if HAVE_HDF5_INT2FLOAT_CONVERSIONS
+      // hdf5 currently doesn't support float/integer conversions
+      else
+	{
+	  double max_val, min_val;
+
+	  if (m.all_integers (max_val, min_val))
+	    save_type_id
+	      = save_type_to_hdf5 (get_save_type (max_val, min_val));
+	}
+#endif /* HAVE_HDF5_INT2FLOAT_CONVERSIONS */
+
+      type_id = hdf5_make_complex_type (save_type_id);
+      if (type_id < 0) goto error_cleanup;
+
+      data_id = H5Dcreate (loc_id, name.c_str (), 
+			   type_id, space_id, H5P_DEFAULT);
+      if (data_id < 0)
+	goto error_cleanup;
+    
+      hid_t complex_type_id = hdf5_make_complex_type (H5T_NATIVE_DOUBLE);
+      if (complex_type_id < 0)
+	goto error_cleanup;
+
+      Complex *mtmp = m.fortran_vec ();
+      if (H5Dwrite (data_id, complex_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+		    (void*) X_CAST (double *, mtmp)) < 0)
+	{
+	  H5Tclose (complex_type_id);
+	  goto error_cleanup;
+	}
+
+      H5Tclose (complex_type_id);
+    }
+  else if (tc.is_list ())
+    {
+      data_id = H5Gcreate (loc_id, name.c_str (), 0);
+      if (data_id < 0)
+	goto error_cleanup;
+
+      data_is_group = 1;
+
+      // recursively add each element of the list to this group
+      octave_value_list lst = tc.list_value ();
+
+      for (int i = 0; i < lst.length (); ++i)
+	{
+	  // should we use lst.name_tags () to label the elements?
+	  char s[20];
+	  sprintf (s, "_%d", i);
+	  bool retval2 = add_hdf5_data (data_id, lst (i), s, "",
+					false, save_as_floats);
+	  if (! retval2)
+	    goto error_cleanup;
+	}
+
+      // mark with an attribute "OCTAVE_LIST" with value 1
+      // to distinguish from structures (also stored as HDF5 groups):
+      if (hdf5_add_attr (data_id, "OCTAVE_LIST") < 0)
+	goto error_cleanup;
+    }
+  else if (tc.is_map ())
+    {
+      // an Octave structure
+      data_id = H5Gcreate (loc_id, name.c_str (), 0);
+      if (data_id < 0)
+	goto error_cleanup;
+
+      data_is_group = 1;
+
+      // recursively add each element of the structure to this group
+      Octave_map m = tc.map_value ();
+      Octave_map::iterator i = m.begin ();
+      while (i != m.end ())
+	{
+	  // XXX FIXME XXX -- if the length of the structure array is
+	  // 1, should we really create a list object?
+	  bool retval2 = add_hdf5_data (data_id, octave_value (m.contents (i)),
+					m.key (i), "", false, save_as_floats);
+	  if (! retval2)
+	    goto error_cleanup;
+
+	  i++;
+	}
+    }
+  else
+    {
+      gripe_wrong_type_arg ("save", tc, false);
+      goto error_cleanup;
+    }
+
+  // attach doc string as comment:
+  if (doc.length () > 0
+      && H5Gset_comment (loc_id, name.c_str (), doc.c_str ()) < 0)
+    goto error_cleanup;
+
+  retval = 1;
+
+  // if it's global, add an attribute "OCTAVE_GLOBAL" with value 1
+  if (mark_as_global)
+    retval = hdf5_add_attr (data_id, "OCTAVE_GLOBAL") >= 0;
+
+ error_cleanup:
+
+  if (! retval)
+    error ("save: error while writing `%s' to hdf5 file", name.c_str ());
+
+  if (data_id >= 0)
+    {
+      if (data_is_group)
+	H5Gclose (data_id);
+      else
+	H5Dclose (data_id);
+    }
+
+  if (space_id >= 0)
+    H5Sclose (space_id);
+
+  if (type_id >= 0)
+    H5Tclose (type_id);
+
+  return retval;
+}
+
+// Write data from TC in HDF5 (binary) format to the stream OS,
+// which must be an hdf5_ofstream, returning true on success.
+
+bool
+save_hdf5_data (std::ostream& os, const octave_value& tc,
+		const std::string& name, const std::string& doc,
+		bool mark_as_global, bool save_as_floats)
+{
+  hdf5_ofstream& hs = (hdf5_ofstream&) os;
+
+  return add_hdf5_data (hs.file_id, tc, name, doc,
+			mark_as_global, save_as_floats);
+}
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/