python - 从 .npy 文件中读取数组到 Fortran 90

Question

我正在使用 Python 以二维数组的形式生成一些初始数据，例如“X”，然后使用 Fortran 对它们进行一些计算。最初，当数组大小约为 10,000 x 10,000 时，np.savetxt 在速度方面运行良好。但是一旦我开始增加数组的维度，savetxt 的速度就会显着降低。所以我尝试了 np.save ，这导致保存速度更快，但文件以 .npy 格式保存。我将如何在 Fortran 中读取这样的文件以重建原始数组？根据我的阅读，二进制通常会导致最低的空间消耗和最快的速度。

在 Fortran 90 中，

open(10,file='/home/X.npy')

read(10,*)X

我收到以下错误:Fortran 运行时错误:列表输入的第 1 项中的实数错误

编辑:在 python 中我这样做，

import numpy as np 
x = np.linspace(-50,50,10)
y = np.linspace(-50,50,10) 
X,Y = np.meshgrid(x,y) 
np.save('X',X) 

然后在 Fortran 中我这样做，

program read_unformatted_binary_from_python
    use iso_c_binding
    implicit none

    integer, parameter :: DPR = selected_real_kind(p=15)
    character(len=*), parameter :: filename = 'X.npy'

    integer, parameter :: N = 10
    real(DPR), allocatable, dimension(:, :) :: X

    allocate(X(N, N))


    open(40, file=filename, status='old', access='stream', form='unformatted')
    read(40) X
    close(40)

    write(*,*) X

end program read_unformatted_binary_from_python

输出以 1.8758506894003703E-309 1.1711999948023422E+171 5.2274167985502976E-037 8.4474009688929314E+252 2.651412321034566210E+180E+252 2.651412321034566210E+180 73E+247 2.1620996769994603E+233 7.5805790251297605E-096 3.4756671925988047E-152 6.5549091408576423E-260 -50.000000000000000 - 38.888888888888886 -27.777777777777779 等等

使用 .bin 格式时不会发生此错误，仅当使用 np.save 时会导致 npy。

Answer 1

Vladimir F 是正确的，您希望在 Fortran 中对“原始二进制”文件进行“流式”访问。这是一个 MWE:

python

import numpy as np
A = np.random.rand(10000, 10000)
print(A.sum())
A.tofile('data.bin')

中文

program read_unformatted_binary_from_python
    use iso_c_binding
    implicit none

    integer, parameter :: DPR = selected_real_kind(p=15)
    character(len=*), parameter :: filename = 'data.bin'

    integer, parameter :: N = 10000
    real(DPR), allocatable, dimension(:, :) :: dat

    allocate(dat(N, N))

    open(40, file=filename, status='old', access='stream', form='unformatted')
    read(40) dat
    close(40)

    write(*,*) sum(dat)

end program read_unformatted_binary_from_python

我的 Fortran 示例可能比必要的要长一些，因为我使用许多不同的系统和编译器套件，并且也不喜欢大型静态数组(毕竟我是 Fortran 用户)。

我在 MacBook Pro 上使用 Python 2.7.x、Numpy 13.x 和来自 Homebrew GCC 6.3.0_1 的 gfortran 快速编写了代码，但这应该适用于所有系统。

更新:这里需要特别注意阵列的形状和大小。如果分配的 dat 大于文件中的内容，则流式 read 应该尝试填充整个数组，点击 EOF符号，并发出错误。在 Python 中，np.fromfile() 方法将读取到 EOF 然后返回具有适当长度的一维数组，即使 A 是最初是多维的。这是因为原始二进制文件没有元数据，只是来自 RAM 的连续字节串。

因此，Python 中的以下行生成相同的文件:

A = np.random.rand(10000, 10000)
A.tofile('file.whatever')
A.ravel().tofile('file.whatever')
A.reshape((100, 1000, 1000)).tofile('file.whatever')

并且该文件可以读入并 reshape 为:

B = np.fromfile('file.whatever').reshape(A.shape)
B = np.fromfile('file.whatever').reshape((100, 1000, 100, 10))
# or something like
B = np.fromfile('file.whatever') # just a 1D array
B.resize(A.shape)  # resized in-place

在 Fortran 中，使用流访问很容易读取整个原始文件，而无需提前知道其大小，尽管您显然需要某种用户输入来 reshape 数据:

program read_unformatted_binary_from_python
    use iso_c_binding
    implicit none

    integer, parameter :: DPR = selected_real_kind(p=15)
    character(len=*), parameter :: filename = 'data.bin'
    integer :: N = 10000, bytes, reals, M
    real(DPR), allocatable :: A(:,:), D(:, :), zeros(:)
    real(DPR), allocatable, target :: B(:)
    real(DPR), pointer :: C(:, :)

    allocate(A(N, N))

    open(40, file=filename, status='old', access='stream', form='unformatted')

    read(40) A
    write(*,*) 'sum of A', sum(A)

    inquire(unit=40, size=bytes)
    reals = bytes/8
    allocate(B(reals))

    read(40, pos=1) B
    write(*,*) 'sum of B', sum(B)

    ! now reshape B in-place assuming the user wants the first dimension 
    ! (which would be the outer dimension in Python) to be length 100
    N = 100
    if(mod(reals, N) == 0) then
         M = reals/N
         call C_F_POINTER (C_LOC(B), C, [N, M])
         write(*, *) 'sum of C', sum(C)
         write(*, *) 'shape of C', shape(C)
    else
         write(*,*) 'file size is not divisible by N!, did not point C to B'
    end if

    ! now reshape B out-of-place with first dimension length 9900, and
    ! pad the result so that there is no size mismatch error  
    N = 9900
    M = reals/N
    if(mod(reals, N) > 0) M=M+1

    allocate(D(N, M))
    allocate(zeros(N), source=real(0.0, DPR))
    D = reshape(B, [N, M], pad=zeros)

    write(*,*) 'sum of D', sum(D)
    write(*,*) 'shape of D', shape(D)

    ! obviously you can also work with shape(A) in fortran the same way you
    ! would use A.shape() in Python, if you already knew that A was the
    ! correct shape of the data
    deallocate(D)
    allocate(D, mold=A)
    D = reshape(B, shape(A))
    write(*,*) 'sum of D', sum(D)
    write(*,*) 'shape of D', shape(D)

    ! or, just directly read it in, skipping the whole reshape B part
    read(40, pos=1) D
    write(*,*) 'sum of D', sum(D)

    close(40)

end program read_unformatted_binary_from_python