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?

Answer 1

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...