如何在 C 中将动态分配的 3D 数组写入 hdf5 文件？

Posted 2023-02-21

【中文标题】如何在 C 中将动态分配的 3D 数组写入 hdf5 文件？

【英文标题】：How to write dynamically allocated 3D array into hdf5 file in C?

【发布时间】：2022-01-14 23:09:37

【问题描述】：

我有一个动态分配的 3D 数组，实现为 指向数组指针数组的指针（至少这是我对我正在做的事情的解释），并希望将该数据存储在 hdf5 文件中。虽然文件中存储了某些内容，但它不是原始数据。

这是我的代码（这里省略了错误检查的内容）：

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

double ***arr3D_d( size_t dim1, size_t dim2, size_t dim3 ) 
    size_t  ii, jj;
    double  ***arr;

    arr = calloc( (size_t)dim1, sizeof(double**) );
    for ( ii=0 ; ii<dim1 ; ++ii ) 
        arr[ii] = calloc( (size_t)(dim2*dim3), sizeof(double*) );
        for ( jj=0 ; jj<dim2 ; ++jj ) 
            arr[ii][jj] = calloc( (size_t)(dim3), sizeof(double) );
        
    
    return arr;


int main( int argc, char *argv[] ) 
    size_t  ii, jj, kk,
            dim1, dim2, dim3;
    double  ***arr3D;

    // hdf5 related variables
    hid_t   file_id, dataset_id, dataspace_id;
    hsize_t dims[3];
    herr_t  status;

    dim1    = 2;
    dim2    = 3;
    dim3    = 4;
    arr3D   = arr3D_d( dim1, dim2, dim3 );

    for (ii=0 ; ii<dim1 ; ++ii)
        for (jj=0 ; jj<dim2 ; ++jj)
            for (kk=0 ; kk<dim3 ; ++kk)
                arr3D[ii][jj][kk]   = ii + jj + kk;

    for (ii=0 ; ii<dim1 ; ++ii)
        for (jj=0 ; jj<dim2 ; ++jj)
            for (kk=0 ; kk<dim3 ; ++kk)
                printf( "arr3D[%ld][%ld][%ld] = %f\n", 
                        ii, jj, kk, arr3D[ii][jj][kk] ); 

    // create new file for hdf5 data to be written into
    file_id = H5Fcreate( "data.h5", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT );
    // create simple dataspace for the dataset
    dims[0] = dim1;
    dims[1] = dim2;
    dims[2] = dim3;
    dataspace_id    = H5Screate_simple( 3, dims, NULL );
    // create dataset
    dataset_id      = H5Dcreate( file_id, "dataset", H5T_NATIVE_DOUBLE, dataspace_id, 
        H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT );
    // write the dataset
    status          = H5Dwrite( dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, 
        H5P_DEFAULT, arr3D[0][0] );
    // terminate access and free identifiers
    status          = H5Dclose(dataset_id);
    status          = H5Sclose(dataspace_id);
    status          = H5Fclose(file_id);

    return 0;

当我现在用h5dump 输出数据时，它的内容如下：

HDF5 "data.h5" 
GROUP "/" 
   DATASET "dataset" 
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE  ( 2, 3, 4 ) / ( 2, 3, 4 ) 
      DATA 
      (0,0,0): 0, 1, 2, 3,
      (0,1,0): 0, 2.42092e-322, 1, 2,
      (0,2,0): 3, 4, 0, 2.42092e-322,
      (1,0,0): 2, 3, 4, 5,
      (1,1,0): 0, 5.58294e-322, 4.64561e-310, 4.64561e-310,
      (1,2,0): 4.64561e-310, 0, 0, 0
      
   

这与代码中的arr3D 不对应，它在运行时打印到控制台 - 输出内容为：

arr3D[0][0][0] = 0.000000
arr3D[0][0][1] = 1.000000
arr3D[0][0][2] = 2.000000
arr3D[0][0][3] = 3.000000
arr3D[0][1][0] = 1.000000
arr3D[0][1][1] = 2.000000
arr3D[0][1][2] = 3.000000
arr3D[0][1][3] = 4.000000
arr3D[0][2][0] = 2.000000
arr3D[0][2][1] = 3.000000
arr3D[0][2][2] = 4.000000
arr3D[0][2][3] = 5.000000
arr3D[1][0][0] = 1.000000
arr3D[1][0][1] = 2.000000
arr3D[1][0][2] = 3.000000
arr3D[1][0][3] = 4.000000
arr3D[1][1][0] = 2.000000
arr3D[1][1][1] = 3.000000
arr3D[1][1][2] = 4.000000
arr3D[1][1][3] = 5.000000
arr3D[1][2][0] = 3.000000
arr3D[1][2][1] = 4.000000
arr3D[1][2][2] = 5.000000
arr3D[1][2][3] = 6.000000

如上所述，这不是写入 hdf5 文件的内容。我做错了什么？

【问题讨论】：

为第二个“维度”分配 dim2*dim3 元素是没有意义的。这将使它成为dim1 x (dim2*dim3) x dim3 元素的矩阵。

【参考方案1】：

垃圾值的原因是在 HDF5 文件中存储指针而不是实际的双精度值。为了很好地存储连续分配的数组，最简单的方法是使用指向可变长度数组 (VLA) 的指针分配 3d 数组。只需将所有分配代码替换为以下行：

 double (*arr3D)[dim2][dim3] = calloc(dim1, sizeof *arr3D);

就是这样。

记得拨打free(arr3D)释放它。

与流行的看法相反，将 VLA 添加到 C 的主要原因是为了简化处理多维数组，而不是用于运行时定义大小的对象的堆栈分配。

使用H5Dwrite() 存储数组时，只需将arr3D 作为最后一个参数。

修改后的hdf5文件内容为：

HDF5 "data.h5" 
GROUP "/" 
   DATASET "dataset" 
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE  ( 2, 3, 4 ) / ( 2, 3, 4 ) 
      DATA 
      (0,0,0): 0, 1, 2, 3,
      (0,1,0): 1, 2, 3, 4,
      (0,2,0): 2, 3, 4, 5,
      (1,0,0): 1, 2, 3, 4,
      (1,1,0): 2, 3, 4, 5,
      (1,2,0): 3, 4, 5, 6
      
   

【讨论】：

有趣，我不熟悉您为 3D 数组分配内存的方式。如果我直接将其写入main，它就可以工作。愚蠢的问题，但是当把它放到一个单独的函数中时会是什么样子......？

@Alf，通常你不会。数组的维度绑定到数组的type。 VLA 类型的问题在于，完成此类型需要对大小表达式进行运行时评估（即[dim2]）。并且运行时表达式可以仅在块范围进行评估。因此，不可能在文件范围内定义 VLA 类型，因此函数无法返回 VLA 类型。解决方法是返回void* 并将其转换回double (*)[dim2][dim3]，但它与显式解决方案没有太大区别。

再次感谢。我是否正确，而不是将 VLA 传递给函数时，我需要 int some_funct( size_t dim2, size_t dim3, double *arr3D[dim2][dim3] ) 之类的东西？

@Alf，这实际上很容易。函数参数在块范围内评估。只需使用int foo(size_t dim1, size_t dim2, size_t dim3, double arr3D[dim1][dim2][dim3])