! --------------------------------------------------------
! matrix-matrix multiply function with our own REAL
! --------------------------------------------------------
    MODULE MyRealModule
      TYPE MyReal
        REAL x			! <---- Our own REAL type
      END TYPE MyReal
    END MODULE MyRealModule   

   function MatMatMult(A, B) result (C)
    use MyRealModule
    implicit none

    TYPE(MyReal), dimension(:,:), INTENT(IN) :: A
    TYPE(MyReal), dimension(:,:), INTENT(IN) :: B
    TYPE(MyReal), 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(:,:).x = 0.0

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


  program Main
   use MyRealModule
   implicit none

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

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

   DO i=1, 3
     DO j = 1, 3
       A(i,j).x = 0
     END DO
     A(i,i).x = 1
   END DO

   DO i=1, 3
     DO j = 1, 3
       CALL random_number(B(i,j).x)
     END DO
   END DO

   C = A * B

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

  end program

