Fortran Coder

标题: 两个矩阵的直积 [打印本页]

作者: hongcha    时间: 2014-9-3 11:03
标题: 两个矩阵的直积

[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    时间: 2014-9-3 13:57
如果用数组片段的话,这个命题就更简单了。

[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    时间: 2014-9-3 19:48
fcode 发表于 2014-9-3 13:57
如果用数组片段的话,这个命题就更简单了。

[mw_shl_code=fortran,true]Module KroneckerProduct_Mod

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

简单的我来,难的你上
作者: aliouying    时间: 2014-9-4 09:41
fcode 发表于 2014-9-3 19:53
简单的我来,难的你上

简单的还能凑凑份子,难的我就不会了,还得你来啊~
作者: 山林悍匪    时间: 2014-9-11 23:15
每天学习一点
作者: mangix2010    时间: 2014-10-22 04:10
mark~~不错的例子
作者: xiaoshunxiaohu    时间: 2014-12-31 12:40
很好的讨论啊!




欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2