I have nine matrices whose dimension as (N by N)
A1(i,j),A2(i,j),A3(i,j),A4(i,j),A5(i,j),A6(i,j),A7(i,j),A8(i,j),A9(i,j)
Then I want to construct a larger matrix (3N by 3N) including these nine matrices as:
A = [A1 A2 A3
A4 A5 A6
A7 A8 A9]
In fortran, can I use the command line as
do i=1,FN
do j=1,FML
A(i,j) = [A1(i,j),A2(i,j),A3(i,j);A4(i,j),A5(i,j),A6(i,j);A7(i,j),A8(i,j),A9(i,j)]
end do
end do
Just for fun, you can also make the large A matrix by using do-loops as
do i = 1, N
A( i, : ) = [ A1( i,: ), A2( i,: ), A3( i,: ) ]
A( i + N, : ) = [ A4( i,: ), A5( i,: ), A6( i,: ) ]
A( i + N*2, : ) = [ A7( i,: ), A8( i,: ), A9( i,: ) ]
enddo
which fills the A matrix in row-major way and so the small matrices also appear in that way. If really really necessary, this could also be written as one-liner as
A = transpose( reshape( &
[ ( [ A1( i,: ), A2( i,: ), A3( i,: ) ], i=1,N ), &
( [ A4( i,: ), A5( i,: ), A6( i,: ) ], i=1,N ), &
( [ A7( i,: ), A8( i,: ), A9( i,: ) ], i=1,N ) ], [N*3, N*3] ))
which turns out to be the transpose of the second array constructor in the @francescalus answer (in one-liner form)
A = reshape( &
[ ( [ A1( :,i ), A4( :,i ), A7( :,i ) ], i=1,N ), &
( [ A2( :,i ), A5( :,i ), A8( :,i ) ], i=1,N ), &
( [ A3( :,i ), A6( :,i ), A9( :,i ) ], i=1,N ) ], [N*3, N*3] )
To go one-step further, we may define hcat and vcat routines as in other languages (note here that explicit interface is necessary):
function hcat( A, B, C ) result( X )
integer, dimension(:,:) :: A, B, C
integer :: X( size(A,1), size(A,2)+size(B,2)+size(C,2) )
X = reshape( [ A, B, C ], shape( X ) )
endfunction
function vcat( A, B, C ) result( X )
integer, dimension(:,:) :: A, B, C
integer :: X( size(A,1)+size(B,1)+size(C,1), size(A,2) )
X = transpose( reshape( &
[ transpose(A), transpose(B), transpose(C) ], &
[ size(X,2), size(X,1) ] ) )
endfunction
then we can write
A = vcat( hcat( A1, A2, A3 ), hcat( A4, A5, A6 ), hcat( A7, A8, A9 ) )
which is somewhat more similar to the desired form in the question:
A = [ A1 A2 A3 ; A4 A5 A6 ; A7 A8 A9 ]
Although Fortran is helpful when it comes to array manipulation, creation of block matrices isn't as elegant as you want from your examples (and would come from certain other languages).
It's possible to use an array constructors to create your desired matrix, much as would be done with scalar elements. That is, RESHAPE([A1, A2, A3, ..., A9],[3*N,3*N]) will give you a 3*Nx3*N matrix. It's just that it won't be the one you want.
As with the other question/answers the array constructor [...] considers the array element order to create a rank-1 array of length 9*N**2, which is then reshaped to a square matrix. Where those other examples used scalar elements here you have arrays for the data in the constructor. The elements of the constructor are themselves taken in array element order, which would be equivalent to
[A1(1,1), A1(2,1), ..., A1(1,2), A1(2,2), ..., A2(1,1), ... ]
which is unwanted.
So, a constructor would be something like
[A1(:,1), A4(:,1), A7(:,1), A1(:,2), ..., A6(:,3), A9(:,3)]
which works, but is unwieldy.
There may be other tricks to get something more "elegant" into the constructor, but as Vladimir F comments, it may be much nicer to just do direct assignment to the various blocks:
A(1:N,1:N) = A1
A(1:N,N+1:2*N) = A2
A(1:N,2*N+1:3*N) = A3
A(N+1:2*N,1:N) = A4
....
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
Too long by 4774 characters
! [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
