Fortran reshape - N-dimensional transpose

Question

I'm trying to write some code in Fortran which requires the re-ordering of an n-dimensional array. I thought the reshape intrinsic combined with the order argument should allow this, however I'm running into difficulties.

Consider the following minimal example

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

I would expect this to result in

Ad1 : 1010111.00       1010112.00       1010113.00               3           5           7          11          13

Bd4 : 1010111.00       1010112.00       1010113.00               7           5          11           3          13

But instead I find:

Ad1 : 1010111.00       1010112.00       1010113.00               3           5           7          11          13

Bd4 : 1010111.00       1010442.00       1020123.00               7           5          11           3          13

I've tried this with gfortran (4.8.3 and 4.9) and ifort (11.0) and find the same results, so it's likely that I am simply misunderstanding something about how reshape works.

Can anybody shed any light on where I'm going wrong and how I can achieve my goal?

Hydroxylamine

Answer 1

When order= is specified in reshape the elements of the result taken with permuted subscript order correspond to the elements of the source array. That probably isn't entirely clear. The Fortran 2008 standard states this as (ignoring the part about pad=)

The elements of the result, taken in permuted subscript order ORDER (1), ..., ORDER (n), are those of SOURCE in normal array element order ..

What this means is that from your example with order=[3,2,4,1,5] there is the mapping to

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

of

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

with offset changing most rapidly in the third index of matB corresponding to most rapidly varying in the first of matA. The next fastest varying in matB being dimension 2, then 4, and so on.

So, it's the elements matB(1,1,1:3,1,1) that correspond the matA(:,1,1,1,1).

I've been explicit in the extent of that matB slice because you've a problem with the shape of matB: you want the shape of matB to be the inverse of the permutation given by the order= specifier.

You could write your example as

  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

Alternatively, if it's the shape of matB that you want, then it's the order permutation that wants inverting:

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

---

At first glance, it may be natural to view the order relating to the dimensions of the source. However, the following may clarify: the result of the reshaping is the same regardless of the shape of source (what is used are the elements of the array in natural order); the order= value has size equal to that of the shape= value. For the first of these, if the source were, say [1,2,3,4,5,6] (recall how we construct rank-2 arrays), then order= could never have any effect (it would have to be [1]) if it were used on the source.

Deva

Many thanks for this, it has clarified the working of order. I had read the specification but I think it was a case of reading what I expected to see rather than what it actually said. It means if we haveN=reshape(M,shape(N),order=[a,b,c]) then the first dimension of M becomes the ath dimension of N etc. – Hydroxylamine

It is natural to expect that the order is of the dimensions of the source. However, perhaps the following may help change that intuition (or remind): the result of the reshaping is the same regardless of the shape of source (what is used are the elements of the array in natural order); the order= value has size equal to that of the shape= value. For the first of these, if source were, say [1,2,3,4,5,6] (recall construction of 2-D arrays), then order= could never have any effect (it would have to be [1]) if it were used on the source. – Deva

That does indeed help my understanding, thank you again. – Hydroxylamine

Answer 2

Because I also feel the behavior of order for multi-dimensional arrays is quite non-intuitive, I made some code comparison below to make the situation even clear (in addition to the already complete @francescalus' answer). First, in a simple case, reshape() with and without order gives the following:

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  ]

This example shows that without order the elements are filled in the usual column-major way, while with order=[2,1] the 2nd dimension runs faster and so the elements are filled row-wise. The key point here is that the order specifies which dimension of the LHS (rather than the source array) runs faster (as emphasized in the above answer).

Now we apply the same mechanism to higher-dimensional cases. First, reshape() of the 5-dimensional array without order

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

corresponds to the explicit loops

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

where matA_seq is a sequential view of matA

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

Now attaching order=[3,2,4,1,5] to reshape(),

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

then the order of DO-loops is changed such that

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

This means that the 3rd dimension of matB (and thus i4) runs fastest (which corresponds to the second line in the above Answer). But what is desired by OP is

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

which corresponds to

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

i.e., the final line of the francescalus' answer.

Hope this comparison further clarifies the situation...

Theatrics

Thanks I agree it is not very intuitive, although I do find your use of the "flattened" view of matA does help clarify the use of order further. – Hydroxylamine

Thanks. Your explanation is very useful. at the end, in the order option, we pick the place of the target array according to the original sequence. – Connett