我有一个用Fortran编写的基于MPI的程序,它在每个节点( 2D时间序列的部分)生成复杂数据的3D数组。我想使用并行I/O将这些数组写到一个文件中,该文件可以相对容易地在python中打开,以便进行进一步的分析/可视化。理想情况下,我希望解决方案是内存高效的(即避免创建中间临时数组)。
使用NetCDF,我已经成功地改编了一个subroutine,它可以实现实数的3D数组。然而,当涉及到复杂数组时,我遇到了一个绊脚石。
在下面的代码中,我尝试通过创建一个由两个实数组成的复合数据类型,并假设Fortran复杂数据类型的实部和虚部连续存储在3D数组的第一维中,将子例程从实数扩展到复数。
module IO
use NetCDF
use, intrinsic :: iso_fortran_env, only: dp => real64
implicit none
contains
subroutine output_3D(dataname, starts, ends, global_data_dims, &
local_data, MPI_communicator)
character(len=*), intent(in) :: dataname
integer, dimension(3), intent(in) :: starts
integer, dimension(3), intent(in) :: ends
integer, dimension(3), intent(in) :: global_data_dims
complex(dp), intent(in) :: local_data( 1:(ends(1) - starts(1)+ 1), &
1:(ends(2) - starts(2) + 1), &
1:(ends(3) - starts(3) + 1))
integer, dimension(3) :: expanded_starts
integer, intent(in) :: MPI_communicator
integer :: ncid, varid, dimid(3)
integer :: counts(3)
integer :: typeid
expanded_starts(1) = (starts(1))* 2 + 1
expanded_starts = starts(2)
expanded_starts(3) = starts(3)
call check(nf90_create( trim(dataname)//'.cdf', &
IOR(NF90_NETCDF4, NF90_MPIIO), &
ncid, &
comm = MPI_communicator, &
info = MPI_INFO_NULL))
call check(nf90_def_dim(ncid, "x", global_data_dims(1), dimid(1)))
call check(nf90_def_dim(ncid, "y", global_data_dims(2) * 2, dimid(2)))
call check(nf90_def_dim(ncid, "z", global_data_dims(3), dimid(3)))
! define a complex data type consisting of two real(8)
call check(nf90_def_compound(ncid, 16, "COMPLEX", typeid))
call check(nf90_insert_compound(ncid, typeid, "REAL", 0, NF90_DOUBLE))
call check(nf90_insert_compound(ncid, typeid, "IMAG", 8, NF90_DOUBLE))
! define a 3D variable of type "complex"
call check(nf90_def_var(ncid, dataname, typeid, dimid, varid))
! exit define mode
call check(nf90_enddef(ncid))
! Now in NETCDF data mode
! set to use MPI/PnetCDF collective I/O
call check(nf90_var_par_access(ncid, varid, NF90_COLLECTIVE))
counts(1) = (ends(1) - starts(1) + 1) * 2
counts(2) = (ends(2) - starts(2) + 1)
counts(3) = (ends(3) - starts(3) + 1)
call check(nf90_put_var(ncid, &
varid, &
local_data, &
start = expanded_starts,&
count = counts))
! close the file
call check(nf90_close(ncid))
return
end subroutine output_3D
subroutine check(status)
integer, intent ( in) :: status
if(status /= nf90_noerr) then
print *, trim(nf90_strerror(status))
stop 2
end if
end subroutine check
end module IO
program test_write
use IO
use MPI
complex(dp) :: data(2,2,3)
integer :: flock
integer :: rank
integer :: ierr
integer :: i, j, k
call MPI_init(ierr)
call MPI_comm_size(MPI_comm_world, flock, ierr)
call MPI_comm_rank(MPI_comm_world, rank, ierr)
do k = 1, 3
do j = 1, 2
do i = 1, 2
data(i,j,k) = cmplx(i, j, 8)
enddo
enddo
enddo
if (rank == 0) then
call output_3D_hdf5('out', [1,1,1], [2,2,3], [2,2,6], &
data, MPI_comm_world)
else
call output_3D_hdf5('out', [1,1,4], [2,2,6], [2,2,6], &
data, MPI_comm_world)
endif
call MPI_finalize(ierr)
end program test_write
上面的代码在编译时会导致“没有用于nf90_put_var的特定函数”错误。这表明该函数对输入数组的数据类型不满意,因此很明显,在复合数据类型的使用方面,我遗漏了一些东西。
编辑:一个简单的解决方法是将复杂的数组赋给一个真实的指针,如this post中所述。然后可以使用numpy对数组进行整形/重新转换,以得到python中的复杂数组。它有点笨拙,也有点令人不满意--但对于我现在的目的来说,可能已经足够好了。
发布于 2019-07-21 19:19:22
这只是你将在下面看到的部分原因的答案-但它对于评论来说太长了。希望我能够找到丢失的信息并“升级”它,但这是我到目前为止所拥有的。
如果您在https://www.unidata.ucar.edu/software/netcdf/fortran/docs/f90-user-defined-data-types.html#f90-compound-types-introduction的“复合类型简介”下查看NetCDF4文档,您将看到:
若要在复合类型中写入数据,请首先使用
nf90_def_compound
创建类型,多次调用nf90_insert_compound
以添加到复合类型中,然后使用适当的nf90_put_var1
、nf90_put_vara
、nf90_put_vars
或nf90_put_varm call
写入数据。
注意,它根本没有提到nf90_put_var
,而是四个不同的函数。这在某种程度上是有道理的,nf90_put_var
应该很好地重载了所有的内部类型NetCDF支持(它完全是垃圾,它不支持复杂的),所以对于非内部类型,可能有一些类似C的接口,比如void *
,我猜上面提到的四个函数就是这样实现的。
到目前为止还不错,您应该调用nf90_put_var1
、nf90_put_vara
、nf90_put_vars
或nf90_put_varm
中的一个,而不是nf90_put_var
。现在坏消息是-我找不到这4个函数的任何文档。等价的C函数在https://www.unidata.ucar.edu/software/netcdf/docs/group__variables.html上,所以你也许能够从那里找出需要什么,但它不是很好-我至少向Unidata提交了一个错误报告,但对我来说,缺乏对complex的内在支持足以让我在其他地方寻找我的I/O解决方案……
当我在这里的时候,你真的不应该对变量的种类使用显式的数字,我可以向你展示complex(8)
将无法编译的编译器。取而代之的是使用Selected_real_kind
或类似的,或者使用内部模块iso_fortran_env中的常量,或者可能使用iso_c_binding中的常量,以及复数的类型与组成它的实数的类型相同的事实。
https://stackoverflow.com/questions/57130311
复制相似问题