[原创] matrix乘法

Program matrix 
     implicit none
     integer n
     parameter (n=16) ! remember setup dimension first!!!
       real:: a(n,n), b(n,n), c(n,n)
       open(11,status='unknown',file='c.DAT')
       read(11,*) c 
       open(12,status='unknown',file='b.DAT')  
       read(12,*) b
     call matmul(a,b,n,n,c,n,n)
       write(*,*) 'Your new matrix looks like: '
     write(*,*)
       call output (a,n)
     write(*,*)
     stop
     end
       subroutine output (a,n)
       implicit none
     integer n
     real a(n,n)
     integer i,j
       character(len = 100) :: for = '(??(1X, F8.1))'  
       write (for(2:3), '(i2)') n
       do i=1,n
     write (*, fmt=for) (a(i,j), j=1,n)
       open(21,status='unknown',file='a.dat')  ! your result is here
       write(21,fmt=for) (a(i,j),j=1,n)
     enddo
     return
     end subroutine output
       subroutine matmul (a, b, br, bc, c, cr, cc)
     implicit none
     integer br  
     integer bc
     real b(br, bc) 
     integer cr 
     integer cc  
     real c(cr,cc)   
     real a(br,cc) 
     integer i, j, k  
       if (bc .ne. cr) then
     write (*,*) 'size error :('
     stop
     endif
     do i=1, br
     do j=1, cc
     a(i,j)=0
     do k=1, bc
       a(i,j) = a(i,j) +b(j,k) * c(k,i)
       enddo
     enddo
     enddo
     return
     end

Module MatMulMod
  Implicit None
contains

  Function mat_mul(a, b) result( c )
    Real a(:,:) , b(:,:)
    Real c( size(a,dim=1) , size(b,dim=2) )
    Integer i, j, k
    c = 0.0
    Do i = 1, size(a,dim=1)
      Do j = 1, size(b,dim=2)
        Do k = 1, size(a,dim=2)
          c(i, j) = c(i, j) + a(i, k) * b(k, j)
        End Do
      End Do
    End Do
  End Function mat_mul
  
  Subroutine output(a)
    Real a(:,:)
    Integer i , n    
    Character (Len=100) :: for = '(??(1X, es13.6))'
    n = size( a , dim = 1 )
    Write (for(2:3), '(i2.2)') n
    !Open (21, Status='unknown', File='c.dat') ! your result is here
    Do i = 1, n
      Write (*, for) a(i,:)      
      !Write (21, for ) a(i,:)
    End Do
    !Close( 21 )
  End Subroutine output

End Module MatMulMod

Program matrix
  use MatMulMod
  Implicit None
  Integer n
  Parameter (n=2) ! remember setup dimension first!!!
  Real :: a(n, n), b(n, n), c(n, n)
  a = reshape([1,2,3,4],[2,2]) ; b=reshape([2,3,4,5],[2,2])
  !Open (11, Status='unknown', File='a.DAT')
  !Read (11, *) c
  !Open (12, Status='unknown', File='b.DAT')
  !Read (12, *) b
  c = mat_mul(a, b)
  Write (*, *) 'Your new matrix looks like: '
  Call output(c)
  c = matmul(a,b)
  Call output(c)
End Program matrix