对不同形状的数组使用单个虚拟参数

Question

假设我们有几个不同的连续数组

real :: a1(2),a2(2),a3(2)
real :: b1(2,2), b2(2,2),b3(2,2)
real :: c1(3,3,3), c2(3,3,3),c3(3,3,3)

和一个子程序

subroutine mtpy(x,y,z)
  real,contiguous, intent(in) :: x(:), y(:)
  real, intent(out) :: z(size(x))

  z=x*y
end subroutine mtpy

如何在以下一系列调用中使用mtpy：

call mtpy(a1,a2,a3)
call mtpy(b1,b2,b3)
call mtpy(c1,c2,c3)

显然，这将导致编译器错误，因为实际参数和虚拟参数的形状不匹配。在这样的情况下，我曾经声明几个特定的​​过程，每个过程都处理一个特定的形状，然后使用一个接口将它们包装起来。然而，这是相当乏味的（想象一下，有大量简单的元素函数和将多维数组（最多三个维度）视为单维数组的纯过程，然后为每个数组提供诸如 sub_1d、sub_2d、sub_3d 等实例尽管他们所有人实际上都在做同样的工作）。

我想，一个部分解决方案是使用 RESHAPE

call mtpy(reshape(b1,[4]),reshape(b2,[4]),bb)

但是，我可以确定编译器（我主要对 gfortran 和 ifort 感兴趣）不会开始创建 1d 临时对象来保存重构的 b1 和 b2 数组吗？

现在，我也知道可以声明一个数组指针，例如

real, pointer, contiguous :: p1(:),p2(:),p3(:)

并进行如下指针赋值如

p1(1:size(c1))=>c1

但是，这种方法有一个缺点，我需要将原始数组声明为目标。这不会影响编译器能够执行的优化吗？

我想，另一种解决方案是使用假定大小的数组，但我注意到 Metcalf 等人称其用法为“已弃用”，而且我不确定对优化的影响。

那么，有没有一种简单的方法可以将多维 fortran 数组视为一维数组（在子例程或函数中），并且不会对该数组强加不必要的假设（例如 TARGET）？如果我可以使用 RESHAPE 而不必担心创建临时对象（我只处理连续数组），我会这样做。有什么建议吗？

Answer 1

未来的 Fortran 2018 标准将提供假定秩数组（这将允许接收任何秩的数组）和一个 select rank 构造，它可以很容易地允许使用一个假定秩数组来解决这种情况（参见，例如, The new features of Fortran 2018 从第 16 页开始），并且使用多个假定秩数组更困难。

假定大小的数组虽然不流行或不推荐，但它们是当前标准 (Fortran 2008) 以及下一个标准草案 (Fortran 2018) 的有效且不会过时的功能，因此可以在需要时使用它们。由于很多 Fortran 77 代码都依赖于此，而且其中大部分代码已有数十年的历史，因此我预计它会在大多数编译器中得到显着优化。

但是，您不需要使用假定大小的数组，您可以使用显式形状数组（具有显式维度的数组），并且只要数组的实际参数有足够的元素，代码就会有效，因为，根据 2008 年标准第 12.5.2.11 节第 4 段，

如果虚拟参数是显式形状或假定大小的数组，则表示元素序列并对应于数组的虚拟参数的实际参数是与虚拟参数相关联的序列。实参的秩和形状不必与哑元的秩和形状一致，但哑元中的元素个数不得超过实参的元素序列中的元素个数。如果 dummy 参数是假定大小，则 dummy 参数中的元素数正好是元素序列中的元素数。

所以你可以

call mtpy(a1,a2,a3,size(a3))
call mtpy(b1,b2,b3,size(b3))
call mtpy(c1,c2,c3,size(c3))
...
subroutine mtpy(x,y,z,n)
  integer, intent(in) :: n
  real, intent(in) :: x(n), y(n)
  real, intent(out) :: z(n)
  z=x*y
end subroutine mtpy

Answer 2

因为我也不确定reshape() 是否会在连续情况下创建一个临时数组，所以我尝试打印原始地址并通过c_loc() 传递数组。然后，即使对于小的一维数组，gfortran-8 和 ifort-16 中的reshape() 似乎也会创建临时对象（因为第一个元素的地址不同）。因此，假设即使对于简单的情况也会创建临时对象似乎更安全（有关更多信息，请参阅下面 francescalus 的 cmets。）

module test
    use iso_c_binding, only: c_loc
    implicit none

    interface linear
        module procedure linear_r2d, linear_r3d
    endinterface
contains

subroutine calc_ver1( a )  !! assumed-shape dummy array
    real, contiguous, target :: a(:)
    print *, "addr = ", c_loc( a(1) )
    print *, "vals = ", a
endsubroutine

subroutine calc_ver2( a, n )  !! explicit-shape dummy array
    integer      :: n
    real, target :: a( n )
    print *, "addr = ", c_loc( a(1) )
    print *, "vals = ", a
endsubroutine

function linear_r2d( a ) result( ptr )  !! returns a 1-d pointer from 2-d array
    real, contiguous, target  :: a(:,:)
    real, contiguous, pointer :: ptr(:)
    ptr( 1 : size(a) ) => a
endfunction

function linear_r3d( a ) result( ptr )  !! returns a 1-d pointer from 3-d array
    real, contiguous, target  :: a(:,:,:)
    real, contiguous, pointer :: ptr(:)
    ptr( 1 : size(a) ) => a
endfunction

endmodule

program main
    use test
    implicit none
    integer i
    real, target :: a(2), b(2,2), c(2,2,2)

    a = [1,2]
    b = reshape( [( 2*i, i=1,4 )], [2,2]   )
    c = reshape( [( 3*i, i=1,8 )], [2,2,2] )

    print *, "addr(a) = ", c_loc( a(1) )
    print *, "addr(b) = ", c_loc( b(1,1) )
    print *, "addr(c) = ", c_loc( c(1,1,1) )

    print *, "[ use assumed-shape dummy ]"
    call calc_ver1( a )
    ! call calc_ver1( b )  ! rank mismatch
    ! call calc_ver1( c )  ! rank mismatch

    print *, "--- with reshape() ---"
    call calc_ver1( reshape( b, [size(b)] ) )
    call calc_ver1( reshape( c, [size(c)] ) )

    print *, "--- with linear() ---"
    call calc_ver1( linear( b ) )
    call calc_ver1( linear( c ) )

    print *
    print *, "[ use explicit-shape dummy ]"
    call calc_ver2( a, size(a) )
    call calc_ver2( b, size(b) )
    call calc_ver2( c, size(c) )
end

ifort-16 在 Linux 上的结果：

 addr(a) =                7040528
 addr(b) =                7040544
 addr(c) =                7040560
 [ use assumed-shape dummy ]
 addr =                7040528
 vals =    1.000000       2.000000    
 --- with reshape() ---
 addr =        140736361693536
 vals =    2.000000       4.000000       6.000000       8.000000    
 addr =        140736361693560
 vals =    3.000000       6.000000       9.000000       12.00000       15.00000       18.00000       21.00000       24.00000  
 --- with linear() ---
 addr =                7040544
 vals =    2.000000       4.000000       6.000000       8.000000    
 addr =                7040560
 vals =    3.000000       6.000000       9.000000       12.00000       15.00000       18.00000       21.00000       24.00000    

 [ use explicit-shape dummy ]
 addr =                7040528
 vals =    1.000000       2.000000    
 addr =                7040544
 vals =    2.000000       4.000000       6.000000       8.000000    
 addr =                7040560
 vals =    3.000000       6.000000       9.000000       12.00000       15.00000       18.00000       21.00000       24.00000    

gfortran-8 在 OSX10.11 上的结果：

 addr(a) =       140734555734776
 addr(b) =       140734555734752
 addr(c) =       140734555734720
 [ use assumed-shape dummy ]
 addr =       140734555734776
 vals =    1.00000000       2.00000000    
 --- with reshape() ---
 addr =       140734555734672
 vals =    2.00000000       4.00000000       6.00000000       8.00000000    
 addr =       140734555733984
 vals =    3.00000000       6.00000000       9.00000000       12.0000000       15.0000000       18.0000000       21.0000000       24.0000000    
 --- with linear() ---
 addr =       140734555734752
 vals =    2.00000000       4.00000000       6.00000000       8.00000000    
 addr =       140734555734720
 vals =    3.00000000       6.00000000       9.00000000       12.0000000       15.0000000       18.0000000       21.0000000       24.0000000    

 [ use explicit-shape dummy ]
 addr =       140734555734776
 vals =    1.00000000       2.00000000    
 addr =       140734555734752
 vals =    2.00000000       4.00000000       6.00000000       8.00000000    
 addr =       140734555734720
 vals =    3.00000000       6.00000000       9.00000000       12.0000000       15.0000000       18.0000000       21.0000000       24.0000000

而且我还认为，根据情况，显式形状的虚拟数组很有用，并且问题中的代码似乎正是这种情况。 （因为实际参数是连续的，所以没有临时创建数组。）如果不需要calc_ver2() 中的大小参数n，我们可以使用返回一维数组指针的函数（参见上面的linear() ），但考虑到calc_ver2() 的简单性，我想这可能有点矫枉过正。（顺便说一句，我在代码的各个位置附加了target，这仅仅是因为c_loc() 需要它）。