两个矩阵的直积

hongcha

[Fortran] 纯文本查看 复制代码

subroutine product_AB(H,A,B,m,n,r,s,index1,index2)
  implicit none
  integer(8)::m,n,r,s,ii,jj,kk,ll,row,column,index1,index2       
  complex(8)::A(m,n),B(r,s),H(index1,index2)
  row=0
  do ii=1,m
    do kk=1,r             
      row=row+1         
      column=0
      do jj=1,n
        do ll=1,s
          column=column+1
          h(row,column)=a(ii,jj)*b(kk,ll)
        enddo
      end do
    enddo
  enddo
end subroutine product_AB

fcode

如果用数组片段的话，这个命题就更简单了。

[Fortran] 纯文本查看 复制代码

Module KroneckerProduct_Mod
  Implicit None
  Integer , parameter , private :: DP = Selected_Real_Kind( 9 )
contains

  Subroutine KroneckerProduct( A , B , H )
    Real( Kind = DP ) , Intent( IN )  :: A(:,:) , B(:,:)
    Real( Kind = DP ) , Intent( OUT ) :: H(:,:)
    Integer :: i , j , m , n , p , q
    m = size( A , dim = 1 )
    n = size( A , dim = 2 )
    p = size( B , dim = 1 )
    q = size( B , dim = 2 )
    Do i = 1 , m
      Do j = 1 , n
        H( p*(i-1)+1 : p*i , q*(j-1)+1 : q*j ) = B * A(i,j)
      End Do
    End Do
  End Subroutine KroneckerProduct

End Module KroneckerProduct_Mod

Program www_fcode_cn
  use KroneckerProduct_Mod
  Implicit None
  Integer , parameter :: DP = Selected_Real_Kind( 9 )
  Integer , parameter :: m=2 , n=3 , p=4, q=5 , index1 = m*p , index2 = n*q
  Real(kind=8) :: A(m,n) , B(p,q) , H(index1,index2)
  integer :: i
  A = reshape( (/1,2,3,4,5,6/) , (/2,3/) )
  B = reshape( (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20/) , (/4,5/) )
  call KroneckerProduct( A , B , H )
  Do i = 1 , index2
    Write(*,'(8(f5.1,1x))') H( :, i)
  End do
End Program www_fcode_cn

aliouying

fcode 发表于 2014-9-3 13:57
如果用数组片段的话，这个命题就更简单了。

[mw_shl_code=fortran,true]Module KroneckerProduct_Mod

你这回复得太及时了，根本不给我们机会啊~

fcode

aliouying 发表于 2014-9-3 19:48
你这回复得太及时了，根本不给我们机会啊~

简单的我来，难的你上

aliouying

fcode 发表于 2014-9-3 19:53
简单的我来，难的你上

简单的还能凑凑份子，难的我就不会了，还得你来啊~

山林悍匪

每天学习一点

mangix2010

mark~~不错的例子

xiaoshunxiaohu

很好的讨论啊！