Fortran重塑 - N维转置

Question

我正在尝试编写一些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 的工作原理存在误解。请问有谁能够帮忙指出我的错误并告诉我如何实现我的目标吗？

- d_1999

Answer 1

因为我也觉得对于多维数组，order的行为相当不直观，所以我在下面做了一些代码比较，以使情况更加清晰（除了已经完整的@francescalus'的答案）。首先，在一个简单的情况下，有和没有order的reshape()分别如下：

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的五维数组进行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] )

希望这个比较能够进一步澄清情况...即francescalus回答的最后一行。

- roygvib

---

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

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

Answer 2

当在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的形状与order=指定器给出的排列的倒数相同。

您可以将示例编写为：

  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]（回想一下我们如何构建秩为 2 的数组），那么如果在源上使用 order=，它永远不会产生任何影响（必须是 [1]）。

- francescalus

---

非常感谢这个解释，它阐明了 order 的工作原理。我已经阅读了规范，但我认为这是一种读取我期望看到的内容而不是实际内容的情况。这意味着如果我们有N=reshape(M,shape(N),order=[a,b,c])，则 M 的第一个维度成为 N 的 ath 维度，依此类推。 - d_1999

可以自然地期望顺序与源的维度相同。然而，也许以下内容有助于改变这种直觉（或提醒）：重新塑形的结果与源的形状无关（使用的是数组元素的自然顺序）；order=值的大小等于shape=值的大小。对于第一个问题，如果源为[1、2、3、4、5、6]（请回想一下二维数组的构造），那么如果在源上使用order=，它永远不会起作用，因为它必须为[1]。 - francescalus

确实有助于我的理解，再次感谢你。 - d_1999