PROGRAM MatrMult 
   IMPLICIT NONE

   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)
   
   DO j = 1, N
      DO i = 1, L
         C(i,j) = 0.0
         DO k = 1, M
            C(i,j) = C(i,j) + A(i,k) * B(k,j)
         END DO
      END DO
   END DO

   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)
   
   DO j = 1, N
      DO i = 1, L
         C(i,j) = DDOT(M,At(1,i),1,B(1,j),1)
      END DO
   END DO

   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