Fortran - 可分配派生类型的可分配数组

Posted 2023-02-21

【中文标题】Fortran - 可分配派生类型的可分配数组

【英文标题】：Fortran - Allocatable array of allocatable derived type

【发布时间】：2020-03-17 10:54:32

【问题描述】：

所以我过去 3-4 天一直在寻找，但找不到这个问题的答案。 它与特定派生类型的可分配数组有关。这都是计算流体动力学求解器的一部分。但是，实际应用并不重要。让我提供一些（快速）编程上下文。

假设我们有一个简单的模块，它定义了一个固定大小的派生类型，主程序分配了一个包含多种类型的数组：

module types

   integer, parameter  :: equations = 10

   type array_t
      double precision :: variable(equations) ! variables to solve
      double precision :: gradient(equations,3) ! gradient of variables in x,y,z direction
      double precision :: limiter(equations) ! limiter of variables
   end type

end module

program test

   use types

   implicit none

   type(array_t), allocatable :: array(:)
   integer :: elements

   elements = 100
   allocate(array(elements))

end program

当然，这段代码 sn-p 可以使用每个编译器进行编译。由于 array_t 的大小是固定的并且在编译时是已知的，因此我们只需在一行中分配结构 array（定义结构内的 array_t 重复次数）。

当涉及到内存位置时，变量将按如下方式存储：

array(1)%variable(1) ! element 1
array(1)%variable(2)
...
...
array(1)%gradient(1,1) ! the rest of this 2D array will be written column-major in fortran
array(1)%gradient(2,1)
array(1)%gradient(3,1)
...
...
array(1)%limiter(1)
array(1)%limiter(2)
...
...
array(2)%variable(1) ! element 2
array(2)%variable(2)
...
...

在本例中，我们设置参数equations=10。在求解器中，我们始终将此大小保留为最大值 (10)：所有派生类型都具有求解器可能需要的最大维度。不幸的是，这意味着我们分配的内存可能比我们实际需要的更多——一些模拟可能只需要 5 或 6 个方程而不是 10 个。此外，派生类型维度保持固定在大小 10 的事实使求解器在求解较少方程时变慢 - 未使用的内存位置将减少内存带宽 -。

我想做的是利用具有 allocatable 属性的派生类型。这样，我可以只使用所需数量的方程（即 array_t 的维度）来分配结构，这些方程将在运行时（而不是在编译时）定义，并且会根据模拟参数而改变。

看看下面的代码sn-p：

module types

   integer, save:: equations

   type array_t
      double precision, allocatable :: variable(:) ! variables to solve
      double precision, allocatable :: gradient(:,:) ! gradient
      double precision, allocatable :: limiter(:) ! limiter of variables
   end type

end module

program test

   use types

   implicit none

   type(array_t), allocatable :: array(:)
   integer :: i,elements

   equations = 10
   elements = 100
   allocate(array(elements))
   do i=1,elements
      allocate(array(i)%variable(equations))
      allocate(array(i)%gradient(equations,3))
      allocate(array(i)%limiter(equations))
   enddo

end program

这是，到目前为止，唯一的方法我设法使它工作。求解器运行并收敛，这意味着语法不仅可编译而且等效于使用固定大小。

但是，即使对于相同数量的方程，这种方法的求解器速度也会明显变慢。

这意味着存在内存错位。根据测量的运行时间，变量的存储方式似乎与使用固定大小时不同。

在第二种方法中，变量如何存储在全局内存中？我想实现与第一种方法相同的模式。我感觉像是分配结构的第一行

allocate(array(elements))

不知道要分配什么，因此要么分配一大块内存（以适应稍后出现的可分配类型），要么仅将索引 array(1) 分配给 array(elements) 而没有其他内容（这意味着结构的实际内容稍后存储在循环内）。

如何让第二种方法像第一种方法一样存储变量？

编辑 #1

由于参数化派生类型得到了一些关注，我认为发布一些额外的细节会很有用。

参数化派生类型将在数组分配在主程序内部的情况下工作（如我发布的示例代码）。

然而，我的“真实世界”案例更像是这样的：

(file_modules.f90)
module types

   integer, parameter  :: equations = 10

   type array_t
      double precision :: variable(equations) ! variables to solve
      double precision :: gradient(equations,3) ! gradient pf variables
      double precision :: limiter(equations) ! limiter of variables
   end type

end module

module flow_solution
   use types
   type (array_t), allocatable, save :: cell_solution(:)
end module

(file_main.f90)
program test

   use flow_solution

   implicit none
   integer :: elements

   elements = 100
   allocate(cell_solution(elements))

end program

这些（如您所料）通过 Makefile 单独编译和链接。 如果我使用参数化派生类型，则无法编译模块文件，因为在编译时不知道类型的大小“n”。

编辑 #2

有人建议我提供具有参数化派生类型的工作和非工作代码示例。

工作示例：

module types

   integer, parameter  :: equations=10

   type array_t(n)
      integer, len     :: n
      double precision :: variable(n) ! variables to solve
      double precision :: gradient(n,n) ! gradient
      double precision :: limiter(n) ! limiter of variables
   end type

end module

module flow_solution
    use types
    type(array_t(equations)), allocatable, save :: flowsol(:)
end module

program test

   use flow_solution

   implicit none

   integer :: i,elements

   elements = 100 
   allocate(flowsol(elements))

end program

非工作示例：

module types

   integer, save       :: equations

   type array_t(n)
      integer, len     :: n
      double precision :: variable(n) ! variables to solve
      double precision :: gradient(n,n) ! gradient
      double precision :: limiter(n) ! limiter of variables
   end type

end module

module flow_solution
    use types
    type(array_t(equations)), allocatable, save :: flowsol(:)
end module

program test

   use flow_solution

   implicit none

   integer :: i,elements

   equations = 10
   elements = 100 
   allocate(flowsol(elements))

end program

编译器（ifort）错误：

test.f90(16): error #6754: An automatic object must not appear in a SAVE statement or be declared with the SAVE attribute.   [FLOWSOL]
    type(array_t(equations)), allocatable, save :: flowsol(:)
---------------------------------------------------^
test.f90(16): error #6841: An automatic object must not appear in the specification part of a module.   [FLOWSOL]
    type(array_t(equations)), allocatable, save :: flowsol(:)
---------------------------------------------------^
compilation aborted for test.f90 (code 1)

我应该以不同的方式声明/分配数组吗？

【问题讨论】：

您可能正在寻找 Fortran 的参数化派生类型。我没有时间写答案，但是关于 SO 上的这个主题有一些 Q 和 As，还有很多其他资源在网上。

给定类似type(array_t(:)), allocatable :: cell_solution(:) 的相应分配语句allocate(array_t(5) :: cell_solution(100)) 似乎是合适的。这就是 array_t 是 length 参数化类型的地方（此注释中未显示）。这种长度参数化的类型在这里似乎适用。

感谢您的意见！你的建议是可编译的。但是，我必须在求解器中实现它，看看它是否正常工作。问题是我在使用可分配类型方法时遇到了内存错位问题。我希望这种方法以正确的方式存储变量。

如果您关心派生类型在内存中的布局，您可能需要调查sequence 语句。

【参考方案1】：

麻烦的线是

type(array_t(equations)), allocatable, save :: flowsol(:)

在这里，您希望将长度参数化类型作为保存对象。作为编译器对象，保存的对象不能是自动对象。此外，您不能在模块范围内拥有自动对象。

现在，为什么flowsol 是一个自动对象？它是自动的，因为它的类型是 array_t(equations)：equations 不是常量。相反，您希望延迟类型参数（就像数组的形状一样）：

type(array_t(:)), allocatable, save :: flowsol(:)

当你来分配这样一个对象时你需要提供长度类型参数值：

allocate(array_t(equations) :: flowsol(elements))

如您所见，当equations 是（命名的）常量时，事情会更快乐。在这种情况下，不会声明自动对象。

【讨论】：

不幸的是，这个选项也很慢。乍一看，它似乎比可分配方法还要慢，但是我没有时间比较两者。我可以肯定的是，它肯定比原始方法（具有固定方程=10 作为参数的方法）慢得多。由于这个结构是一个相对较大的求解器的一部分，我认为我不可能发布一个展示我所指的减速的代码示例。

对于上下文，我唯一改变的是我将“方程式”变量更改为“整数，保存”而不是“整数，参数”。相应的类型以长度 n 参数化，求解器初始化中的分配如您所建议的那样。请注意，我没有更改任何其他结构，也没有在其他任何地方使用这种类型。即使这样，求解器的速度也几乎慢了 100%。

大幅放缓是出乎意料的，但正如您所说，在 SO 问题的范围内没有什么可看的。也许作为一种测试，您可以将类型 kind 参数化（甚至可能为 10），看看是否仍然存在减速。