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

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)

的一维数组

  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 中复制这段代码，但我不明白读取命令如何将数据分配给 4D 矩阵。谁能给点建议？

Answer 1

首先，Fortran equivalence 和内存布局。我们需要先看看内存布局。为简单起见，我将描述 2D 情况。对于任意维度，概括应该不难理解。

Fortran 数组始终是连续的内存块¹。最左边的索引变化最快（称为列优先顺序）。所以：

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

^{¹几乎总是。在某些情况下，数组切片和 F90+ 指针可能会使该声明不正确。即便如此，那些 "arrays" 由分配给程序的连续内存块支持...}

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

memory_index = i*(N-1) + j
a(memory_index)

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

最终，当您有一个二维数组时，编译器会将您的语句 a(i, j) 转换为 a(i*(N-1) + j)。它可能也做一些边界检查 ^{（例如确保 i 在 1 和 N 之间，但你通常必须使用编译器标志打开它，因为它会稍微影响性能，而且当他们的程序变慢时，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 表示法读取 pg 来获得它。 "flatten" 数组为了将数据读入其中是一个技巧。

在 Python 中，您可以做类似的事情，但这实际上取决于您打算如何使用数组。

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

Translating fortran to python: equivalence command used to generate a matrix from text file

python

arrays

fortran

matrix

equivalence