PROGRAM MatrMult 
   IMPLICIT NONE

   INTERFACE
      FUNCTION DDOT(N,X,INCX,Y,INCY) RESULT(XtY)
         REAL(8) :: XtY  ! Y transpose
         INTEGER, INTENT(IN) :: N, INCX, INCY
         REAL(8), DIMENSION(:), INTENT(IN) :: X,Y
     END FUNCTION DDOT
   END INTERFACE

   REAL(8), DIMENSION(:,:), ALLOCATABLE :: A, At, B, C, rC
   INTEGER :: L, M, N, N0, N1
   REAL ::   cput2, cput1

   INTEGER :: i, j, k, t, st
   REAL(8) :: rmin, rmax
   !
   REAL(8), EXTERNAL :: DDOT
   !
   WRITE(*,"(A)") " Input dimensions of matrices  L, M, N "
   READ*,L, M, N
   WRITE(*,"(A,3(1x,I6))") " Matrices dimensions are: ",L,M,N

   ALLOCATE(A(L,M), B(M,N), C(L,N), STAT=st)
   IF ( st /= 0 ) THEN
   	  WRITE(*,*) "Error allocating matrices "
   	  STOP
   END IF

   A = 0.0; B = 0.0; C = 0.0

   DO i = 1, L
      DO k = 1, M
         A(i,k) = REAL(i-k) / REAL(i+k)
      END DO
   END DO

   DO j = 1, N
      DO k = 1, M
         B(k,j) = 0.01 * REAL(k-j) / REAL(k+j)
      END DO
   END DO

   CALL CPU_TIME(cput1)
   
	!! Compute matrix multiply using DO loops (result in C):
		.	.	.

   CALL CPU_TIME(cput2)

   ! Check for correctness
   ALLOCATE(rC(L,N), STAT=st)
   IF ( st /= 0 ) THEN
      	 WRITE(*,*) "Error allocating matrix rC !!!"
      	 STOP
   END IF
      	   	   
      rC = MATMUL(A,B)

      rmin = MINVAL(rC-C); rmax = MAXVAL(rC-C)

      WRITE(*,*) " rmin, rmax = ",rmin,rmax
    
      WRITE(*,*) "Cpu and elapsed time are: ",(cput2-cput1)
!
	WRITE(*,*) " "
	WRITE(*,*) " Try using BLAS DDOT routine:"
!
!	First of all transpose A matrix
	
   ALLOCATE(At(M,L), STAT=st)
   IF ( st /= 0 ) THEN
   	  WRITE(*,*) "Error allocating At matrix "
   	  STOP
   END IF
   
   DO j = 1, M
      DO i = 1, L
         At(j,i) = A(i,j)
      END DO
   END DO
!	
   CALL CPU_TIME(cput1)
   
!! Call ddot for computing matrix multiply (result in C):
		.	.	.
	 

   CALL CPU_TIME(cput2)
   
   
      rmin = MINVAL(rC-C); rmax = MAXVAL(rC-C)

      WRITE(*,*) " rmin, rmax = ",rmin,rmax
    
      WRITE(*,*) "Cpu and elapsed time are: ",(cput2-cput1)
   STOP
END PROGRAM MatrMult