!
! Matrix multiplication
!
Program op
  implicit none
  real :: start, finish
  integer, parameter :: K=1000
  integer, parameter :: L=500
  integer, parameter :: M=2000
  integer, parameter :: N=10
  real*8 :: A1(K,L), A2(L,M), A3(M,N)
  real*8, dimension(1:K,1:N) :: A123
  real f1, f2, time1, time2, time3

  CALL RANDOM_NUMBER(A1)
  CALL RANDOM_NUMBER(A2)
  CALL RANDOM_NUMBER(A3)
  
  call cpu_time(start)
  a123 = matmul(matmul(A1,A2), A3)
  call cpu_time(finish)
  Time1 = finish-start
  print '(/,"Multiplication (A1*A2)*A3")'
  print '("Time1 = ",f6.3," seconds.")', Time1

  call cpu_time(start)
  a123 = matmul(A1, matmul(A2, A3))
  call cpu_time(finish)
  Time2 = finish-start
  print '(/,"Multiplication A1*(A2*A3)")'
  print '("Time2 = ",f6.3," seconds.")', Time2

  ! optimal multiplication of three matrices
  call cpu_time(start)
  CALL OPT_MAT_MUL(A1, A2, A3, A123, K, L, M, N)
  call cpu_time(finish)
  Time3 = finish-start
  print '(/,"Optimal Multiplication")'
  print '(  "Time3 = ",f6.3," seconds.")',  Time3

  print '(/,"Time1/Time3 = ", f6.1)', Time1/Time3
  print '(  "Time2/Time3 = ", f6.1)', Time2/Time3

end Program op


SUBROUTINE OPT_MAT_MUL(A, B, C, D, K, L, M, N)
  implicit none
  INTEGER :: K, L, M, N
  REAL*8  :: A(K,L), B(L,M), C(M,N), D(K,N)
  ! --------------------------------------------
  !    Optimized matrix multiplication
  !       A(:,:)*B(:,:)*C(:,:) -> D()
  ! --------------------------------------------
  IF ( K*L*M + K*M*N .GE. L*M*N + K*L*N) THEN
     D = MATMUL(A, MATMUL(B, C))
  ELSE 
     D = MATMUL(MATMUL(A, B), C)
  END IF
  RETURN
END SUBROUTINE OPT_MAT_MUL

