Logo Search packages:      
Sourcecode: h5utils version File versions  Download package

arrayh5.c

/* Copyright (c) 2005 Massachusetts Institute of Technology
 *
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 * 
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

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

#include <hdf5.h>

#include "arrayh5.h"

#define CHECK(cond, msg) { if (!(cond)) { fprintf(stderr, "arrayh5 error: %s\n", msg); exit(EXIT_FAILURE); } }

#define CHK_MALLOC(p, t, n) CHECK(p = (t *) malloc(sizeof(t) * (n)), "out of memory")

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

arrayh5 arrayh5_create_withdata(int rank, const int *dims, double *data)
{
     arrayh5 a;
     int i;

     CHECK(rank >= 0, "non-positive rank");
     a.rank = rank;
     
     CHK_MALLOC(a.dims, int, rank);

     a.N = 1;
     for (i = 0; i < rank; ++i) {
        a.dims[i] = dims[i];
        a.N *= dims[i];
     }

     if (data)
        a.data = data;
     else {
        CHK_MALLOC(a.data, double, a.N);
     }
     return a;
}

arrayh5 arrayh5_create(int rank, const int *dims)
{
     return arrayh5_create_withdata(rank, dims, NULL);
}

arrayh5 arrayh5_clone(arrayh5 a)
{
     return arrayh5_create(a.rank, a.dims);
}

void arrayh5_destroy(arrayh5 a)
{
     free(a.dims);
     free(a.data);
}

int arrayh5_conformant(arrayh5 a, arrayh5 b)
{
     int i;

     if (a.rank != b.rank)
        return 0;
     for (i = 0; i < a.rank; ++i)
        if (a.dims[i] != b.dims[i])
             return 0;
     return 1;
}

static void rtranspose(int curdim, int rank, const int *dims,
                   int curindex, int curindex_t,
                   const double *data, double *data_t)
{
     int prod_before = 1, prod_after = 1;
     int i;
     
     if (rank == 0) {
        *data_t = *data;
        return;
     }

     for (i = 0; i < curdim; ++i)
        prod_before *= dims[i];
     for (i = curdim + 1; i < rank; ++i)
        prod_after *= dims[i];

     if (curdim == rank - 1) {
        for (i = 0; i < dims[curdim]; ++i)
             data_t[curindex_t + i * prod_before] = data[curindex + i];
     }
     else {
        for (i = 0; i < dims[curdim]; ++i)
             rtranspose(curdim + 1, rank, dims,
                    curindex + i * prod_after,
                    curindex_t + i * prod_before,
                    data, data_t);
     }
}

void arrayh5_transpose(arrayh5 *a)
{
     double *data_t;
     int i;

     CHK_MALLOC(data_t, double, a->N);
     rtranspose(0, a->rank, a->dims, 0, 0, a->data, data_t);
     free(a->data);
     a->data = data_t;

     for (i = 0; i < a->rank - 1 - i; ++i) {
        int dummy = a->dims[i];
        a->dims[i] = a->dims[a->rank - 1 - i];
        a->dims[a->rank - 1 - i] = dummy;
     }
}

void arrayh5_getrange(arrayh5 a, double *min, double *max)
{
     int i;

     CHECK(a.N > 0, "no elements in array");
     *min = *max = a.data[0];
     for (i = 1; i < a.N; ++i) {
        if (a.data[i] < *min)
             *min = a.data[i];
        if (a.data[i] > *max)
             *max = a.data[i];
     }
}

static herr_t find_dataset(hid_t group_id, const char *name, void *d)
{
     char **dname = (char **) d;
     H5G_stat_t info;

     H5Gget_objinfo(group_id, name, 1, &info);
     if (info.type == H5G_DATASET) {
        CHK_MALLOC(*dname, char, strlen(name) + 1);
        strcpy(*dname, name);
        return 1;
     }
     return 0;
}

typedef enum { NO_ERROR = 0, OPEN_FAILED, NO_DATA, READ_FAILED, SLICE_FAILED,
           INVALID_SLICE, INVALID_RANK, OPEN_DATA_FAILED } arrayh5_err;

const char arrayh5_read_strerror[][100] = {
     "no error",
     "error opening HD5 file",
     "couldn't find data set in HDF5 file",
     "error reading data from HDF5",
     "error reading data slice from HDF5",
     "invalid slice of HDF5 data",
     "non-positive rank in HDF file",
     "error opening data set in HDF file",
};

int arrayh5_read(arrayh5 *a, const char *fname, const char *datapath,
             char **dataname,
             int nslicedims, const int *slicedim_, const int *islice_,
             const int *center_slice)
{
     hid_t file_id = -1, data_id = -1, space_id = -1;
     char *dname = NULL;
     int err = NO_ERROR;
     hsize_t i, rank, *dims_copy, *maxdims;
     int *islice = 0, *slicedim = 0;
     int *dims = 0;

     CHECK(a, "NULL array passed to arrayh5_read");
     a->dims = NULL;
     a->data = NULL;

     file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
     if (file_id < 0) {
        err = OPEN_FAILED;
        goto done;
     }
 
     if (datapath && datapath[0]) {
        CHK_MALLOC(dname, char, strlen(datapath) + 1);
        strcpy(dname, datapath);
     }
     else {
        if (H5Giterate(file_id, "/", NULL, find_dataset, &dname) <= 0) {
             err = NO_DATA;
             goto done;
        }
     }

     data_id = H5Dopen(file_id, dname);
     if (data_id < 0) {
        err = OPEN_DATA_FAILED;
        goto done;
     }

     space_id = H5Dget_space(data_id);
     rank = H5Sget_simple_extent_ndims(space_id);
     if (rank <= 0) {
        err = INVALID_RANK;
        goto done;
     }
     
     CHK_MALLOC(dims, int, rank);
     CHK_MALLOC(dims_copy, hsize_t, rank);
     CHK_MALLOC(maxdims, hsize_t, rank);

     H5Sget_simple_extent_dims(space_id, dims_copy, maxdims);
     for (i = 0; i < rank; ++i)
        dims[i] = dims_copy[i];

     free(maxdims);
     free(dims_copy);
     
     for (i = 0; i < nslicedims && slicedim_[i] == NO_SLICE_DIM; ++i)
        ;

     if (i == nslicedims) { /* no slices */
        *a = arrayh5_create(rank, dims);
        
        if (H5Dread(data_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
                  H5P_DEFAULT, (void *) a->data) < 0) {
             err = READ_FAILED;
             goto done;
        }
     }
     else if (nslicedims > 0) {
        int j, rank2 = rank;
        hssize_t *start;
        hsize_t *count;
        hid_t mem_space_id;
        herr_t readerr;

        CHK_MALLOC(slicedim, int, nslicedims);
        CHK_MALLOC(islice, int, nslicedims);

        for (i = j = 0; i < nslicedims; ++i)
             if (slicedim_[i] != NO_SLICE_DIM) {
                if (slicedim_[i] == LAST_SLICE_DIM)
                   slicedim[j] = rank - 1;
                else
                   slicedim[j] = slicedim_[i];
                if (slicedim[j] < 0 || slicedim[j] >= rank) {
                   err = INVALID_SLICE;
                   goto done;
                }
                islice[j] = islice_[i];
                if (center_slice[i])
                   islice[j] += dims[slicedim[j]] / 2;
                if (islice[j] < 0 || islice[j] >= dims[slicedim[j]]) {
                   err = INVALID_SLICE;
                   goto done;
                }
                j++;
             }
        nslicedims = j;
        
        CHK_MALLOC(start, hssize_t, rank);
        CHK_MALLOC(count, hsize_t, rank);
        
        for (i = 0; i < rank; ++i) {
             count[i] = dims[i];
             start[i] = 0;
        }
        for (i = 0; i < nslicedims; ++i) {
             start[slicedim[i]] = islice[i];
             count[slicedim[i]] = 1;
        }
        
        H5Sselect_hyperslab(space_id, H5S_SELECT_SET,
                        start, NULL, count, NULL);

        for (i = j = 0; i < rank; ++i)
             if (count[i] > 1)
                dims[j++] = count[i];
        rank2 = j;

        *a = arrayh5_create(rank2, dims);
        
        mem_space_id = H5Screate_simple(rank, count, NULL);
        H5Sselect_all(mem_space_id);
        
        readerr = H5Dread(data_id, H5T_NATIVE_DOUBLE, 
                      mem_space_id, space_id, 
                      H5P_DEFAULT, (void *) a->data);
        
        H5Sclose(mem_space_id);
        free(count);
        free(start);
        
        if (readerr < 0) {
             err = SLICE_FAILED;
             goto done;
        }
     }
     else {
        err = INVALID_SLICE;
        goto done;
     }

 done:
     if (err != NO_ERROR)
        arrayh5_destroy(*a);
     free(islice);
     free(slicedim);
     free(dims);
     if (space_id >= 0)
        H5Sclose(space_id);
     if (data_id >= 0)
        H5Dclose(data_id);
     if (dataname)
        *dataname = dname;
     else
        free(dname);
     if (file_id >= 0)
        H5Fclose(file_id);

     return err;
}

static int 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);
}

void arrayh5_write(arrayh5 a, char *filename, char *dataname,
               short append_data)
{
     int i;
     hid_t file_id, space_id, type_id, data_id;
     hsize_t *dims_copy;

     if (append_data)
        file_id = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
     else
        file_id = H5Fcreate(filename, H5F_ACC_TRUNC,
                        H5P_DEFAULT, H5P_DEFAULT);
     CHECK(file_id >= 0, "error opening HDF5 output file");     

     if (dataset_exists(file_id, dataname))
        H5Gunlink(file_id, dataname);  /* delete it */

     CHECK(a.rank > 0, "non-positive rank");
     CHK_MALLOC(dims_copy, hsize_t, a.rank);
     for (i = 0; i < a.rank; ++i)
        dims_copy[i] = a.dims[i];
     space_id = H5Screate_simple(a.rank, dims_copy, NULL);
     free(dims_copy);

     type_id = H5T_NATIVE_DOUBLE;
     data_id = H5Dcreate(file_id, dataname, type_id, space_id, H5P_DEFAULT);
     H5Sclose(space_id);

     H5Dwrite(data_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, a.data);

     H5Dclose(data_id);
     H5Fclose(file_id);
}

int arrayh5_read_rank(const char *fname, const char *datapath, int *rank)
{
     hid_t file_id = -1, data_id = -1, space_id = -1;
     char *dname = NULL;
     int err = NO_ERROR;

     file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
     if (file_id < 0) {
        err = OPEN_FAILED;
        goto done;
     }
 
     if (datapath && datapath[0]) {
        CHK_MALLOC(dname, char, strlen(datapath) + 1);
        strcpy(dname, datapath);
     }
     else {
        if (H5Giterate(file_id, "/", NULL, find_dataset, &dname) <= 0) {
             err = NO_DATA;
             goto done;
        }
     }

     data_id = H5Dopen(file_id, dname);
     if (data_id < 0) {
        err = OPEN_DATA_FAILED;
        goto done;
     }

     space_id = H5Dget_space(data_id);
     *rank = H5Sget_simple_extent_ndims(space_id);

 done:
     if (space_id >= 0)
        H5Sclose(space_id);
     if (data_id >= 0)
        H5Dclose(data_id);
     free(dname);
     if (file_id >= 0)
        H5Fclose(file_id);

     return err;
}

Generated by  Doxygen 1.6.0   Back to index