PROGRAM PROVA_MMULT
   IMPLICIT NONE
   INTEGER, PARAMETER :: l=2, m=10, n=4
   REAL :: a(l,m), b(m,n), c(l,n)
   INTEGER :: i, j
   
   DO i = 1, l
      DO j = 1, m
         a(i,j) = SIN(DBLE(i+j))
      END DO
   END DO
      
   DO i = 1, m
      DO j = 1, n
         b(i,j) = DBLE(i+j)/DBLE(i)
      END DO
   END DO
      
   CALL MMULT(l,m,n,c,a,b)
   WRITE(*,*) " MMULT(a,b) = "
   CALL PMATR2(l,n,c)
   
   c = MATMUL(a,b)
   WRITE(*,*) "      = "
   CALL PMATR2(l,n,c)
   
   STOP
CONTAINS 

   FUNCTION VDOT(n,v,w)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n
      REAL, INTENT(IN), DIMENSION(n) :: v, w
      REAL :: c, VDOT
      INTEGER :: i
      
      c = 0.0
      DO i = 1, 100
         c = c + v(i) * w(i)
      END DO
      
      VDOT = c
      
      RETURN
   END FUNCTION VDOT


   SUBROUTINE MMULT(nl,nm,nn,v,w,z)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: nl, nm, nn
      REAL, INTENT(IN) :: w(nl,nm), z(nm,nn)
      REAL, INTENT(OUT), DIMENSION(nl,nn) :: v
      INTEGER :: i, j
      
      DO j = 1, nn
         DO i = 1, nl
           v(i,j) = DOT_PRODUCT(w(i,:),z(:,j))
         END DO
      END DO
      
      RETURN
   END SUBROUTINE MMULT   
   
   SUBROUTINE PMATR2(l,m,z)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: l, m
      REAL, INTENT(IN), DIMENSION(l,m) :: z
      INTEGER :: i, j
      
         DO i = 1, l
            WRITE(*,*) (z(i,j), j=1,m)
         END DO
      WRITE(*,'(/)')
      RETURN
   END SUBROUTINE PMATR2

END PROGRAM PROVA_MMULT