我正在尝试复制旧的fortran代码，以便可以在python中使用它。最困难的部分是导入以相当复杂的形状存储的数据。
我有一个文本文件，其中的14500个数据点存储在5列中，并用逗号分隔：
9.832e-06, 2.121e-16, 3.862e-21, 2.677e-23, 1.099e-22,
5.314e-23, 1.366e-23, 6.504e-23, 3.936e-23, 1.175e-22,
1.033e-23, 5.262e-23, 3.457e-23, 9.673e-23, 1.903e-22,
3.153e-10, 2.572e-21, 5.350e-23, 4.522e-22, 1.468e-22,
2.195e-23, 2.493e-22, 1.075e-22, 3.525e-22, 1.541e-23,
1.935e-22, 9.220e-23, ..., ..., ...,
fortran代码声明了两个变量来存储这些数据：pg和ping。第一个是4D矩阵，第二个是具有数据点数组长度（14500）的1D数组
parameter(NSIZ=29)
parameter(NTg=5)
parameter(NDg=5)
parameter(NTAUg=20)
real pg(NSIZ,NDg,NTg,NTAUg)
real ping(NSIZ*NDg*NTg*NTAUg)
到目前为止一切顺利，但是现在我有一个等效命令：
equivalence(pg,ping)
据我了解，它将pg矩阵指向ping数组。最后一部分是一个循环，该循环从数据文件中读取行并将其分配给ping数组：
NCOLs=5
NROWS=NTg*NDg*NTAUg*NSIZ / NCOLs
irow=1
do i=1,NROWS
read(nat,*) (ping((irow-1)*NCOLs+k),k=1,NCOLs)
print *, ping(irow)
irow=irow+1
enddo
我想在python中复制此代码，但不了解read命令如何将数据分配给4D矩阵。有人可以提供任何建议吗？
首先，Fortranequivalence和内存布局。我们需要先看一下内存布局。为了简单起见，我将描述2D情况。对于任意维数，泛化不应该太难。
fortran数组始终是内存1的连续块。左索引最远变化最快（称为主列排序）。所以：
a(i, j)
a(i+1, j)
1几乎总是这样。在某些情况下，数组切片和F90 +指针可能会使该语句不正确。即使这样，那些“数组”仍由分配给程序的连续内存块支持...
在内存中相邻。现在，假设我们有一个用length声明的数组a(N, M)。从内存的角度来看，这a(N*M)与索引映射为的数组没有什么不同：
memory_index = i*(N-1) + j
a(memory_index)
（这N-1是因为fortran在索引中默认为从1开始）。
最终，当您拥有2D数组时，编译器会将您的语句转换a(i, j)为a(i*(N-1) + j)。它还可能会进行一些边界检查（例如，确保i介于1和之间N，但是您通常必须使用编译器标志将其打开，因为它会稍微影响性能，并且当程序速度变慢时，Fortran人群真的不喜欢它； -）。
好吧，我们在哪里？我所说的几乎所有内容都是对于编译器而言，N维数组实际上只是具有某些索引重新映射的1维数组。
现在我们可以开始处理了equivalence。C语言中的类似构造就是指针。您的Equivalence语句表示中的第一个元素与中的第一个元素位于ping相同的内存位置pg。这使用户可以读取ping平面数组中的数据，但将其填充到n维数组中（因为最终它们共享相同的存储位置）。值得注意的是，现代的Fortran具有指针。它们并不完全像C指针，但我认为大多数Fortran爱好者都建议您在新代码中使用它们，而不要使用等效代码。
现在，我将如何在python中读取此数据？我可能会做一些非常简单的事情：
def yield_values(filename):
with open(filename) as f_input:
for line in f_input:
for val in line.split(','):
yield float(val)
然后，您可以将其打包到一个numpy数组中：
import numpy as np
arr = np.array(list(yield_values('data.dat')))
最后，您可以调整数组的形状：
# Account for python's Row-major ordering.
arr = arr.reshape((NTAUg, NTg, NDg, NSIZ))
如果您更喜欢fortran的列主要订购方式：
arr = arr.reshape((NSIZ, NDg, NTg, NTAUg), order='F')
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句