Fortran 重塑 - N 维转置

Posted 2023-02-21

【中文标题】Fortran 重塑 - N 维转置

【英文标题】：Fortran reshape - N-dimensional transpose

【发布时间】：2016-09-23 08:46:12

【问题描述】：

我正在尝试在 Fortran 中编写一些需要重新排序 n 维数组的代码。我认为 reshape 内在函数结合 order 参数应该允许这样做，但是我遇到了困难。

考虑以下最小示例

program test
    implicit none
    real, dimension(:,:,:,:,:), allocatable :: matA, matB
    integer, parameter :: n1=3, n2=5, n3=7, n4=11, n5=13
    integer :: i1, i2, i3, i4, i5

    allocate(matA(n1,n2,n3,n4,n5)) !Source array
    allocate(matB(n3,n2,n4,n1,n5)) !Reshaped array

    !Populate matA
    do i5=1, n5
       do i4=1, n4
          do i3=1, n3
             do i2=1, n2
                do i1=1, n1
                   matA(i1,i2,i3,i4,i5) = i1+i2*10+i3*100+i4*10000+i5*1000000
                enddo
             enddo
          enddo
       enddo
    enddo

    print*,"Ad1 : ",matA(:,1,1,1,1),shape(matA)
    matB = reshape(matA, shape(matB), order = [3,2,4,1,5])
    print*,"Bd4 : ",matB(1,1,1,:,1),shape(matB) !Leading dimension of A is the fourth dimension of B
end program test

我希望这会导致

Ad1 : 1010111.00       1010112.00       1010113.00               3           5           7          11          13

Bd4 : 1010111.00       1010112.00       1010113.00               7           5          11           3          13

但我发现：

Ad1 : 1010111.00       1010112.00       1010113.00               3           5           7          11          13

Bd4 : 1010111.00       1010442.00       1020123.00               7           5          11           3          13

我已经用gfortran（4.8.3 和 4.9）和ifort（11.0）尝试了这个并找到了相同的结果，所以我很可能只是误解了 reshape 的工作原理。

任何人都可以阐明我的问题所在以及如何实现我的目标吗？

【问题讨论】：

您可能想要matB 的第三列，所以也许您只是没有在正确的点使用逆排列？ [我认为这是你的问题，但这当然会改变你所期望的答案的其他方面。]

【参考方案1】：

当在reshape 中指定order= 时，以置换下标顺序获取的结果的元素对应于源数组的元素。这可能并不完全清楚。 Fortran 2008 标准将其声明为（忽略关于pad= 的部分）

结果的元素，按照置换的下标顺序 ORDER (1), ..., ORDER (n)，是 SOURCE 的正常数组元素顺序的元素..

这意味着从您的 order=[3,2,4,1,5] 示例中可以映射到

matA(1,1,1,1,1), matA(2,1,1,1,1), matA(3,1,1,1,1), matA(1,2,1,1,1), ...

的

matB(1,1,1,1,1), matB(1,1,2,1,1), matB(1,1,3,1,1), matB(1,1,4,1,1), ...

matB 的第三个索引中偏移量变化最快，对应于matA 的第一个索引中变化最快。 matB 中的下一个最快变化是维度 2，然后是维度 4，依此类推。

所以，是元素 matB(1,1,1:3,1,1) 对应于 matA(:,1,1,1,1)。

我已经明确了matB 切片的范围，因为您对matB 的形状有疑问：您希望matB 的形状与@ 给出的排列相反987654335@ 说明符。

你可以把你的例子写成

  implicit none
  integer, parameter :: n1=3, n2=5, n3=7, n4=11, n5=13
  integer matA(n1,n2,n3,n4,n5)
  integer matB(n4,n2,n1,n3,n5)  ! Inverse of permutation (3 2 4 1 5)
  integer i1, i2, i3, i4, i5

  forall (i1=1:n1, i2=1:n2, i3=1:n3, i4=1:n4, i5=1:n5) &
          matA(i1,i2,i3,i4,i5)=i1+i2*10+i3*100+i4*10000+i5*1000000

  print*,"Ad1 : ",matA(:,1,1,1,1),shape(matA)
  matB = reshape(matA, shape(matB), order = [3,2,4,1,5])
  print*,"Bd3 : ",matB(1,1,:,1,1),shape(matB)

end

或者，如果你想要的是matB 的形状，那么它就是想要反转的顺序排列：

  matB = reshape(matA, shape(matB), order = [4,2,1,3,5])

---

乍一看，查看与源尺寸相关的顺序可能很自然。但是，以下可能会澄清：无论源的形状如何，重塑的结果都是相同的（使用的是自然顺序的数组元素）； order= 值的大小等于 shape= 值的大小。对于其中的第一个，如果源是 [1,2,3,4,5,6]（回想一下我们如何构造 rank-2 数组），那么 order= 如果用于来源。

【讨论】：

非常感谢，它澄清了order 的工作。我已经阅读了规范，但我认为这是阅读我期望看到的内容而不是它实际所说的内容。这意味着如果我们有N=reshape(M,shape(N),order=[a,b,c])，那么M 的第一个维度将成为N 等的ath 维度。

很自然地期望顺序是源的尺寸。但是，也许以下内容可能有助于改变这种直觉（或提醒）：无论源的形状如何，重塑的结果都是相同的（使用的是自然顺序的数组元素）； order= 值的大小等于 shape= 值的大小。对于其中的第一个，如果源是 [1,2,3,4,5,6]（召回二维数组的构造），那么如果在源上使用 order= 将永远不会有任何效果（它必须是 [1]） .

这确实有助于我的理解，再次感谢您。

【参考方案2】：

因为我也觉得order 对多维数组的行为很不直观，所以我在下面做了一些代码比较以使情况更加清楚（除了已经完整的@francescalus 的答案）。首先，在一个简单的例子中，reshape() 有和没有order 给出以下结果：

mat = reshape( [1,2,3,4,5,6,7,8], [2,4] )

=> [ 1  3  5  7  ;
     2  4  6  8  ]

mat = reshape( [1,2,3,4,5,6,7,8], [2,4], order=[2,1] )

=> [ 1  2  3  4  ;
     5  6  7  8  ]

这个例子表明，没有order，元素以通常的列主要方式填充，而使用order=[2,1]，第二维运行得更快，因此元素按行填充。这里的关键是 order 指定 LHS 的哪个维度（而不是源数组）运行得更快（正如上面的答案所强调的）。

现在我们将相同的机制应用于更高维的情况。一、没有order的5维数组的reshape()

matB = reshape( matA, [n3,n2,n4,n1,n5] )

对应于显式循环

k = 0
do i5 = 1, n5   !! (5)-th dimension of LHS
do i1 = 1, n1   !! (4)
do i4 = 1, n4   !! (3)
do i2 = 1, n2   !! (2)
do i3 = 1, n3   !! (1)-st dimension of LHS
    k = k + 1
    matB( i3, i2, i4, i1, i5 ) = matA_seq( k )
enddo;enddo;enddo;enddo;enddo

其中matA_seq 是matA 的顺序视图

real, pointer :: matA_seq(:)
matA_seq( 1 : n1*n2*n3*n4*n5 ) => matA(:,:,:,:,:)

现在将order=[3,2,4,1,5] 附加到reshape()，

matB = reshape( matA, [n3,n2,n4,n1,n5], order = [3,2,4,1,5] )

然后改变 DO 循环的顺序，使得

k = 0
do i5 = 1, n5   !! (5)-th dim of LHS
do i3 = 1, n3   !! (1)
do i1 = 1, n1   !! (4)
do i2 = 1, n2   !! (2)
do i4 = 1, n4   !! (3)-rd dim of LHS
    k = k + 1
    matB( i3, i2, i4, i1, i5 ) = matA_seq( k )
enddo;enddo;enddo;enddo;enddo

这意味着matB（因此i4）的第三维运行速度最快（对应于上述答案中的第二行）。但是OP想要的是

k = 0
do i5 = 1, n5   !! (5)-th dim of LHS
do i4 = 1, n4   !! (3)
do i3 = 1, n3   !! (1)
do i2 = 1, n2   !! (2)
do i1 = 1, n1   !! (4)-th dim of LHS
    k = k + 1
    matB( i3, i2, i4, i1, i5 ) = matA_seq( k )
enddo;enddo;enddo;enddo;enddo

对应

matB = reshape( matA, [n3,n2,n4,n1,n5], order = [4,2,1,3,5] )

即法国人回答的最后一行。

希望这个比较能进一步阐明情况......

【讨论】：

谢谢，我同意这不是很直观，尽管我确实发现您使用 matA 的“扁平化”视图确实有助于进一步阐明 order 的使用。

谢谢。你的解释很有用。最后，在 order 选项中，我们根据原始顺序选择目标数组的位置。