如何构造一个包含 9 个较小矩阵的矩阵

Question

我有九个矩阵，其维度为 (N x 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)

然后我想构造一个更大的矩阵（3N x 3N），包括这九个矩阵：

A = [A1 A2 A3
     A4 A5 A6
     A7 A8 A9]

在fortran中，我可以使用命令行吗

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

Answer 1

虽然 Fortran 在数组操作方面很有帮助，但块矩阵的创建并不像您从示例中想要的那样优雅（并且来自某些其他语言）。

可以使用数组构造函数来创建所需的矩阵，就像使用scalar elements 所做的那样。也就是说，RESHAPE([A1, A2, A3, ..., A9],[3*N,3*N]) 会给你一个3*Nx3*N 矩阵。只是它不会是你想要的。

与其他问题/答案一样，数组构造函数[...] 考虑数组元素顺序以创建长度为9*N**2 的秩为1 的数组，然后将reshaped 转换为方阵。在这些其他示例在这里使用标量元素的地方，您在构造函数中有数据数组。构造函数的元素本身是按数组元素顺序获取的，相当于

[A1(1,1), A1(2,1), ..., A1(1,2), A1(2,2), ..., A2(1,1), ... ]

这是不需要的。

所以，构造函数应该是这样的

[A1(:,1), A4(:,1), A7(:,1), A1(:,2), ..., A6(:,3), A9(:,3)]

有效，但笨拙。

可能还有其他技巧可以让构造函数更“优雅”，但作为 Vladimir F cmets，直接分配给各个块可能会更好：

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

Answer 2

只是为了好玩，你也可以使用 do-loops 作为大 A 矩阵

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

它以行主要方式填充 A 矩阵，因此小矩阵也以这种方式出现。如果真的有必要，这也可以写成单行

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] ))

原来是@francescalus 答案中第二个数组构造函数的转置（单行形式）

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] )

---

更进一步，我们可以像在其他语言中一样定义hcat 和vcat 例程（注意这里需要显式接口）：

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

那么我们可以写

A = vcat( hcat( A1, A2, A3 ), hcat( A4, A5, A6 ), hcat( A7, A8, A9 ) )

这有点类似于问题中所需的形式：

A = [ A1 A2 A3 ; A4 A5 A6 ; A7 A8 A9 ]

Answer 3

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