将fortran翻译成python：用于从文本fi生成矩阵的等价命令

Question

我试图复制一个旧的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）

^{pr2}$

到目前为止还不错，但现在我有一个等价的命令：

  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矩阵的。有人能给点建议吗？在

Answer 1

首先，Fortranequivalence和内存布局。我们需要先看看内存布局。为了简单起见，我将描述二维情况。对于任意维，一般化不应该太难理解。在

fortran数组总是一个连续的内存块¹。最左边的索引变化最快（称为列主顺序）。所以：

a(i, j)
a(i+1, j)

几乎总是。在某些情况下，数组切片和F90+指针可能会使该语句不真实。即使这样，这些“数组”也由分配给程序的连续内存块来支持…

在内存中是相邻的。现在，假设我们有一个用长度a(N, M)声明的数组。从内存的角度来看，这与数组a(N*M)没有什么不同，其中索引的映射是：

^{pr2}$

（N-1是因为fortran在索引中默认为基于1）。在

最终，当您有一个2D数组时，编译器将您的语句a(i, j)转换为a(i*(N-1) + j)。它可能也会进行一些边界检查（例如，确保i在1和{}之间，但通常必须用编译器标志打开它，因为它会轻微影响性能，而且Fortran群众确实不喜欢它在程序减速时；-）。在

好吧，我们在哪？对于编译器来说，N维数组实际上只是一个带有索引重新映射的一维数组。

现在我们可以开始处理equivalence。C中类似的构造是指针。等价语句表示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')

Answer 2

由于equivalence语句，ping中存储的任何数据也存储在pg中，因为它们指向内存中的同一位置。因此，当read命令读取一个值并将其存储在ping中时，可以使用4D表示法读取{}。为了将数据读入数组，这是一个“展平”数组的技巧。在

在Python中，可以执行类似的操作，但这实际上取决于如何使用数组。在