Python 读取未格式化的直接访问 Fortran 90 给出不正确的输出

Python reading unformatted direct access Fortran 90 gives incorrect output

这是数据的写入方式(它是一个二维浮点矩阵。我不确定大小)。

open(unit=51,file='rmsd/'//nn_output,form='unformatted',access='direct',status='replace',&
     recl=Npoints*sizeofreal)

!a bunch of code omitted 

    write(51,rec=idx-nstart+1) (real(dist(jdx)),jdx=1,Npoints)

这是我试图阅读文件的方式,灵感来自对这些类似问题的回答:1, 2

   f = open(inputfilename,'rb')
   field = np.fromfile(f,dtype='float64')

但是结果不正确。矩阵的对角线应该为 0(或非常接近),因为它是一个相似矩阵。我尝试了不同的 dtypes 但我仍然无法得到正确的结果。

编辑:这是我得到的输出

array([  1.17610188e+01,   2.45736970e+02,   7.79741823e+02, ...,
         9.52930627e+05,   8.93743127e+05,   7.64186127e+05])

我还没有为重塑而烦恼,因为值的范围应该在 0 到 ~20 之间。

编辑 2:这是文件前 100 行的 link: https://drive.google.com/file/d/0B2Mz7CoRS5g5SmRxTUg5X19saGs/view?usp=sharing

编辑 3:文件顶部的声明

integer,parameter :: real_kind=8
!valuables for MPI
integer :: comm, nproc, myid, ierr

integer,allocatable :: idneigh(:)
real :: tmp
real,allocatable :: traj(:,:)
real(real_kind),allocatable :: dist(:)
real(real_kind),allocatable :: xx(:,:),yy(:,:),rot(:),weight(:)
character(200) :: nn_traj, nn_output, nn_neigh
integer,parameter :: sizeofreal=4,sizeofinteger=4

未格式化的 Fortran 二进制文件对记录长度进行编码以分隔记录。通常,记录长度将在每个记录之前和之后写入,但如果内存服务于处理器相关的详细信息(如果记录未以这种方式分隔,请参见 post 的后半部分)。查看您 post 编辑的文件,如果您将前 4 个字节解释为整数,将其余字节解释为 32 位浮点值,您将得到:

0000000        881505604    7.302916e+00    8.723415e+00    6.914254e+00
0000020     9.826199e+00    7.044637e+00    8.601265e+00    6.629045e+00
0000040     6.103047e+00    9.476192e+00    9.326468e+00    6.535160e+00
0000060     8.904651e+00    4.710213e+00    6.534080e+00    1.156603e+01
0000100     1.046533e+01    9.343380e+00    8.574672e+00    7.498291e+00
0000120     1.071538e+01    7.138038e+00    5.898036e+00    6.182026e+00
0000140     7.037515e+00    6.418780e+00    6.294755e+00    8.327971e+00
0000160     6.796582e+00    7.397069e+00    6.493272e+00    1.126087e+01
0000200     6.467663e+00    7.178994e+00    7.867798e+00    5.921878e+00

如果您寻找过去的记录长度字段,您可以将记录的其余部分读入 python 变量。您需要指定正确的字节数,因为在记录的末尾会有另一个记录长度,这将作为错误值读入您的数组。 881505604 意味着你的 NPoints 是 220376401(如果这不是真的,然后看 post 的后半部分,这可能是数据而不是记录长度)。

您可以通过以下方式在 python 中读取此数据:

f = open('fortran_output', 'rb')
recl = np.fromfile(f, dtype='int32', count=1)
f.seek(4)
field = np.fromfile(f, dtype='float32')

print('Record length=',recl)
print(field)

这会读取记录长度,向前查找 4 个字节并将文件的其余部分读入 float32 数组。如前所述,您需要为读取指定适当的 count= 以不占用结束记录字段。

此程序为您的输入文件输出以下内容:

Record length= [881505604]
[ 7.30291557  8.72341537  6.91425419 ...,  6.4588294   6.53710747
  6.01582813]

数组中的所有 32 位实数都在 0 到 20 之间,正如您建议的那样适合您的数据,所以这看起来像您想要的。


但是,如果您的编译器未在定界未格式化的直接访问文件的记录中对记录长度进行编码,则输出将被解释为 32 位浮点数:

0000000    2.583639e-07       7.3029156        8.723415        6.914254
0000020        9.826199        7.044637        8.601265        6.629045
0000040        6.103047       9.4761915        9.326468       6.5351596
0000060        8.904651        4.710213         6.53408       11.566033
0000100       10.465328         9.34338        8.574672        7.498291
0000120       10.715377        7.138038       5.8980355        6.182026
0000140       7.0375147         6.41878       6.2947555       8.3279705
0000160       6.7965817       7.3970685       6.4932723       11.260868
0000200        6.467663        7.178994        7.867798       5.9218783
0000220        6.710998         5.71757       6.1372333        5.809089

其中第一个值远小于零,考虑到您说矩阵对角线应为 0,这可能是合适的。

要阅读此内容,您可以像尝试一样在 python 中进行读取,但谨慎的做法是使用 count= 仅在超过正确记录长度时才读取正确的记录长度文件中的一条记录。

使用python代码:

f = open('fortran_output', 'rb')
field = np.fromfile(f, dtype='float32')

print(field)

产生输出

[  2.58363912e-07   7.30291557e+00   8.72341537e+00 ...,   6.45882940e+00
   6.53710747e+00   6.01582813e+00]

这与解释为 32 位浮点数的文件输出匹配。


当用各种 gfortran 和 ifort 版本测试时,这种输出的写作细节通常没有记录分隔符,但可能像 post 的前半部分或其他编译器的不同之处.

我还想重申这个解决方案不是 100% 可信度的警告,因为我们没有您写入此文件的数据来验证。但是您确实拥有该信息,并且您应该验证生成的数据是否正确。同样值得注意的是,在使用默认标志进行编译时,ifort 记录长度以 4 字节字为单位,而 gfortran 以 1 字节为单位。这不会导致输出问题,但如果不在 ifort 中进行补偿,您的文件将比需要的大 4 倍。您可以使用 -assume byterecl 和 ifort 来获取以字节为单位的记录长度。