program reshape_test
implicit none
integer, parameter :: N = 2
integer, dimension(N,N) :: A1,A2,A3,A4,A5,A6,A7,A8,A9
character(20) fmt
integer A(3*N,3*N)
A1 = reshape([11,21,12,22],[N,N])
A2 = reshape([13,23,14,24],[N,N])
A3 = reshape([15,25,16,26],[N,N])
A4 = reshape([31,41,32,42],[N,N])
A5 = reshape([33,43,34,44],[N,N])
A6 = reshape([35,45,36,46],[N,N])
A7 = reshape([51,61,52,62],[N,N])
A8 = reshape([53,63,54,64],[N,N])
A9 = reshape([55,65,56,66],[N,N])
write(fmt,'(*(g0))') '(a/',N,'(i2:1x))'
write(*,fmt) 'A1 = ',transpose(A1)
write(*,fmt) 'A2 = ',transpose(A2)
write(*,fmt) 'A3 = ',transpose(A3)
write(*,fmt) 'A4 = ',transpose(A4)
write(*,fmt) 'A5 = ',transpose(A5)
write(*,fmt) 'A6 = ',transpose(A6)
write(*,fmt) 'A7 = ',transpose(A7)
write(*,fmt) 'A8 = ',transpose(A8)
write(*,fmt) 'A9 = ',transpose(A9)
A = reshape([reshape([A1,A4,A7,A2,A5,A8,A3,A6,A9],[N,3,N,3],order=[1,3,2,4])],[3*N,3*N])
write(fmt,'(*(g0))') '(a/',3*N,'(i2:1x))'
write(*,fmt) 'A = ',transpose(A)
end program reshape_test
4774
! [Sigh]
program reshape_test
!
! PROBLEM: Insert J*K M X N matrices AI(I)%A into a J X K block matrix A.
! The blocks AI(I)%A may be input in column-major or row-major order, and
! may be transposed or direct.
!
! SOLUTION: The RESHAPE intrinsic may perform the required pivot
! operation, if the command is designed carefully.
! Step 1: Work out the strides the inputs in array element order will
! be seen in output matrix.
! Step 2: The first 3 extents of the intermediate array may then be
! obtained by sorting the strides ascending and taking the ratios
! between adjacent strides. The 4th extent is such as to match the
! size of the intermediate array with the inputs.
! Step 3: The order of the indices is set so that each index steps
! through the inputs at its proper stride.
!
! EXAMPLE: In the column-major, direct case the second element of the
! inputs will be AI(1)%A(2,1) which should map to A(2,1) in the output,
! at offset of 1 from A(1,1), thus the first stride is 1.
! The second level block will be AI(1)%A(:,2) which should map to
! A(1:M,2) in the output, at offset of M*J from A(1:M,1), thus the
! second stride is M*J.
! The third level block will be AI(K+1)%A which should map to
! A(M+1:2*M,1:N), at offset M from A(1:M,1:N), so the third stride
! is M.
! The fourth level block is AI(2:K*(J-1)+2:K)%A which should map to
! A(1:M*J,N+1:2*N) at offset of M*N*J from A(1:M*J,1:N), so the fourth
! stride is M*N*J.
! Now strides = [1,M*J,M,M*N*J]
! Sort ascending: [1,M,M*J,M*N*J]
! Take ratios: [M,J,N]
! (M*N)*(J*K)/(M*J*N) = K, so SHAPE = [M,J,N,K].
! The first stride is 1, and that is stride take by the first index.
! The second stride is M*J, taken by the third index.
! The third stride is M, taken by the second index.
! The fourth stride is M*N*J, taken by the fourth index of the
! intermediate array. Hence ORDER = [1,3,2,4]
! Finally, to specify the input blocks in column major order such
! that they look read as row ajor in the output as specified,
! SOURCE = [((AI(IJ*(K-1)+IK)%A,IJ=1,J),IK=1,K)]
! Now that the elements of A are in the right order, a second RESHAPE
! casts the collection into the desired shape.
!
! EXERCISE: Permuting the elements of an array in bit-reversed order is
! another class of pivot operation. Write a RESHAPE invocation that
! performs this pivot operation. What are its limitations?
!
implicit none
integer M,N ! Each A_I is an M X N matrix
integer J,K ! A is a J X K block matrix
type A_type
character(:), allocatable :: A(:,:) ! A block!
end type A_type
type(A_type), allocatable :: AI(:) ! Input blocks
character(:), allocatable :: A(:,:) ! Output block matrix
character(20) test_string ! To find max length of index
integer LI ! Max length of block index
integer LM ! Max length of row index
integer LN ! Max length of column index
integer I ! Block index
integer IM ! Row of A(I)%A
integer IN ! Column of A(I)%A
character(40) fmt ! Variable FORMAT
integer IJ ! Row of A
integer IK ! Column of A
! Define problem dimensions
M = 2
N = 2
J = 3
K = 3
! Get max index lengths
write(test_string,'(i0)') J*K
LI = len_trim(test_string)
write(test_string,'(i0)') M
LM = len_trim(test_string)
write(test_string,'(i0)') N
LN = len_trim(test_string)
! Create FORMAT for array element label
write(fmt,'(4(a,i0))') "('A',i0.",LI,",'(',i",LM,",',',i",LN,",')')"
! Create blocks
allocate(AI(J*K))
do I = 1, J*K
allocate(character(4+LI+LM+LN)::AI(I)%A(M,N))
do IM = 1,M
do IN = 1,N
write(AI(I)%A(IM,IN),fmt) I,IM,IN
end do
end do
end do
! Solve the 4 versions
write(*,'(a)') 'Column-major, direct'
write(*,'(a)') 'Input order'
allocate(character(4+LI+LM+LN)::A(J*M,K*N))
write(fmt,'(3(a,i0))') "(",J*K,"('A',i0.",LI,":','))"
write(*,fmt) ((K*(IJ-1)+IK,IJ=1,J),IK=1,K)
A = reshape([reshape([((AI(K*(IJ-1)+IK)%A,IJ=1,J),IK=1,K)],[M,J,N,K],order=[1,3,2,4])],[M*J,N*K])
write(fmt,'(3(a,i0))') '(',N*K,'(a',len(A),':1x))'
write(*,fmt) transpose(A)
deallocate(A)
write(*,'()')
write(*,'(a)') 'Row-major, direct'
write(*,'(a)') 'Input order'
allocate(character(4+LI+LM+LN)::A(J*M,K*N))
write(fmt,'(3(a,i0))') "(",J*K,"('A',i0.",LI,":','))"
write(*,fmt) (I,I=1,J*K)
A = reshape([reshape([(AI(I)%A,I=1,J*K)],[M,J,N,K],order=[1,3,4,2])],[M*J,N*K])
write(fmt,'(3(a,i0))') '(',N*K,'(a',len(A),':1x))'
write(*,fmt) transpose(A)
deallocate(A)
write(*,'()')
write(*,'(a)') 'Column-major, transposed'
write(*,'(a)') 'Input order'
allocate(character(4+LI+LM+LN)::A(J*N,K*M))
write(fmt,'(3(a,i0))') "(",J*K,"('A',i0.",LI,":','))"
write(*,fmt) ((K*(IJ-1)+IK,IJ=1,J),IK=1,K)
A = reshape([reshape([((AI(K*(IJ-1)+IK)%A,IJ=1,J),IK=1,K)],[N,J,M,K],order=[3,1,2,4])],[N*J,M*K])
write(fmt,'(3(a,i0))') '(',M*K,'(a',len(A),':1x))'
write(*,fmt) transpose(A)
deallocate(A)
write(*,'()')
write(*,'(a)') 'Row-major, transposed'
write(*,'(a)') 'Input order'
allocate(character(4+LI+LM+LN)::A(J*M,K*N))
write(fmt,'(3(a,i0))') "(",J*K,"('A',i0.",LI,":','))"
write(*,fmt) (I,I=1,J*K)
A = reshape([reshape([(AI(I)%A,I=1,J*K)],[N,J,M,K],order=[3,1,4,2])],[N*J,M*K])
write(fmt,'(3(a,i0))') '(',M*K,'(a',len(A),':1x))'
write(*,fmt) transpose(A)
deallocate(A)
end program reshape_test