! --------------------------------------------------------
! matrix-matrix multiply function
! --------------------------------------------------------
   function MatMatMult(A, B) result (C)
    implicit none

    real, dimension(:,:), INTENT(IN) :: A
    real, dimension(:,:), INTENT(IN) :: B
    real, dimension( size(A,1), size(B, 2) ) :: C
    integer :: M, N, X, i, j, k

    M = size(A,1)
    N = size(B,2)
    X = size(A,2)   ! Which is also size(B(:,1))

    C = 0.0       !! clear whole matrix

    DO i = 1, M
       DO j = 1, N
	  DO k = 1, X
             C(i,j) = C(i,j) + A(i,k) * B(k,j)
          END DO
       END DO
    END DO
   end function


  program Main
   implicit none

! -------------------------------------------------------------------------
!     "operator Matrix*Matrix"
! -------------------------------------------------------------------------
   INTERFACE operator(*)
    function MatMatMult(A, B) result (C)     !! <--- Problem A*B exists in F90
     real, dimension(:,:), intent(in) :: A
     real, dimension(:,:), intent(in) :: B
     real, dimension( size(A,1), size(B,2) ) :: C
    end function
   END INTERFACE 

   real, dimension( 3, 3 ) :: A, B, C
   integer  i, j

   DATA A /1, 0, 0, 0, 1, 0, 0, 0, 1 /
   CALL random_number(B)

   C = A * B

   DO i = 1, 3
      write (6, "(9F7.3)" ) A(i,:), B(i,:), C(i,:)
   END DO

  end program

