!
! we test the performance of a very lean type(MyMatrix) for matrix multiplication
!
!#-----------------------------------------------------------------------------------------------------
!ifort -heap-arrays -mkl benchmarkMyMatrix.f90
!#-----------------------------------------------------------------------------------------------------
!#my:       N      MatmulMyMatrix    MatmulArray2
!#-----------------------------------------------------------------------------------------------------
!my        30       0.61004E-05       0.38389E-05
!my        60       0.21718E-04       0.13885E-04
!my       120       0.10824E-03       0.78837E-04
!my       190       0.36937E-03       0.29097E-03
!my       250       0.78433E-03       0.66674E-03
!my       500       0.62739E-02       0.54070E-02
!#-----------------------------------------------------------------------------------------------------
!ifort -O3 -heap-arrays -xAVX  -parallel benchmarkMyMatrix.f90 -mkl=parallel 
!export  OMP_NUM_THREADS=8
!export  MKL_NUM_THREADS=8
!#-----------------------------------------------------------------------------------------------------
!#my:       N      MatmulMyMatrix    MatmulArray2        ratio       1/ratio
!#-----------------------------------------------------------------------------------------------------
!my        30       0.63192E-05       0.38982E-05        1.62        0.62
!my        60       0.21967E-04       0.13776E-04        1.59        0.63
!my       120       0.10932E-03       0.79098E-04        1.38        0.72
!my       190       0.36823E-03       0.30644E-03        1.20        0.83
!my       250       0.78471E-03       0.65853E-03        1.19        0.84
!my       500       0.64193E-02       0.53616E-02        1.20        0.84
!my      1000       0.45260E-01       0.41173E-01        1.10        0.91
!my      2000       0.34830           0.32136            1.08        0.92
!my      4000       2.5440            2.4220             1.05        0.95
!------------------------------------------------------------------------------------------------------
!my      30         0.59557E-05       0.38672E-05        1.54        0.65 
!my      60         0.21323E-04       0.13763E-04        1.55        0.65 
!my      120        0.10740E-03       0.78995E-04        1.36        0.74 
!my      190        0.36715E-03       0.29190E-03        1.26        0.80 
!my      250        0.78403E-03       0.64898E-03        1.21        0.83 
!my      500        0.63756E-02       0.53554E-02        1.19        0.84 
!my      1000       0.46292E-01       0.42395E-01        1.09        0.92 
!my      2000       0.35856           0.33732            1.06        0.94 
!my      4000       2.5467            2.4140             1.05        0.95 
!------------------------------------------------------------------------------------------------------
! Post to intel forum 200531:
!------------------------------------------------------------------------------------------------------
!I define a very simple type for a matrix:
! type               MyMatrix
!  integer                               :: n=0
!  complex(dp), allocatable              :: v(:,:)
! end type           MyMatrix
!Then, operations turn out to have a large overhead  for small N x N matrices. For example, matrix multiplication (which I am mostly interested in)
!
!type(MyMatrix) ::  m1, m2, m3
!
!complex(dp), allocatable ::    a1(:,:), a2(:,:), a3(:,:)
!
!m3 = m1 * m2;  m3 = MyMatmul(m1,m2); a3 = MyMatmul(a2,a3)
!
!for N=60 gives about 1.6 performance penalty. For the more complex matrix type that I have in my real program, containing also procedure components, the penalty can be up to 2.5! The comparison is done using an interface that is simply passing m1 % v and m2 % v to BLAS's zgemm . The operator(*) is overloaded  to MyMatmul (same performance as the function call).
!
! The penalty reduces for large matrices, reaching 10% for N=1000 and 5% for N=4000. In my real program, the penalty remains 10% up to N=16000. Other operations (like scalar multiplication/division and or matrix addition) get a performance penalty up to 10!
!
!The performance is compared using wall time.
!
!Is this unavoidable?
!
!-------------------------------------------------------------------------------------------------------------
!
!I am attaching the simplified source code that focuses on this particular problem
!
!Benchmarks on  i7-4790 CPU @ 3.60GHz, Ubuntu 18.04, 
!
!Fortran Intel(R) 64 Compiler for applications running on Intel(R) 64, Version 19.0.5.281 Build 20190815
!
!ifort -O3 -heap-arrays -xAVX  -parallel benchmarkMyMatrix.f90 -mkl=parallel 
!
!export  OMP_NUM_THREADS=8
!export  MKL_NUM_THREADS=8
!
!and on the ARIS supercomputer (Intel Xeon E5-2680v2)
!
!#    compiler_version: Intel(R) Fortran Intel(R) 64 Compiler for applications running on Intel(R) 64, Version 18.0.0.128 Build 20170811
!
!#    compiler_options: -I/apps/compilers/intel/18.0.0/mkl/include -c -O3 -heap-arrays -xCORE-AVX-I -parallel -mkl=parallel
!
!export  OMP_NUM_THREADS=20
!export  MKL_NUM_THREADS=20
!
!Indicative wall times for each operation: (times are in wall time seconds per matrix multiplication averaged over several thousand random matrices - less statistics for N=4000).
!
!N     MatmulMyMatrix  MatmulArray2  ratio
!30    0.63192E-05     0.38982E-05   1.62
!60    0.21967E-04     0.13776E-04   1.59
!120   0.10932E-03     0.79098E-04   1.38
!190   0.36823E-03     0.30644E-03   1.20
!250   0.78471E-03     0.65853E-03   1.19
!500   0.64193E-02     0.53616E-02   1.20
!1000  0.45260E-01     0.41173E-01   1.10
!2000  0.34830         0.32136       1.08
!4000  2.5440          2.4220        1.05
!
!-----------------------------------------------------------------------------------------------------------------------------------
module              matrix_benchmark
 implicit none
 integer, parameter                     :: dp  = kind(0.0_8)
 complex(dp), parameter                 :: ONE = cmplx(1.0_dp,0.0_dp,kind=dp), ZERO = cmplx(0.0_dp,0.0_dp,kind=dp)
 integer                                :: t_count,t_count_rate,t_count_max, f_mout
 integer, parameter                     :: len_label = 25 !length of string labeling tests
!-----------------------------------------------------------------------------------------------------------------------------------
!A mock Matrix type, to check bottleneck:
 type               MyMatrix
  integer                               :: n=0
  complex(dp), allocatable              :: v(:,:)
 end type           MyMatrix
 interface          MyMatrix
  module procedure                      :: mymatrix_construct
 end interface      MyMatrix
!-----------------------------------------------------------------------------------------------------------------------------------
 interface          operator(*)
  module procedure                      :: mymatrix_mult_mymatrix
 end interface      operator(*)
 interface          lmatmul
  module procedure                      :: mymatrix_mult_mymatrix, array_matmul_array
 end interface      lmatmul
!-----------------------------------------------------------------------------------------------------------------------------------
contains
!-----------------------------------------------------------------------------------------------------------------------------------
 subroutine         main_benchmark
  real(dp)                              :: dt
  integer , allocatable                 :: Ns (:)
  real(dp), allocatable                 :: dts(:)
  real(dp)                              :: t_tot=0.0_dp

  call init_benchmark

  Ns     = [30, 60, 120, 190, 250, 500, 1000, 2000, 4000]
  dts    = [10, 10,  10,  10,  10,  20,  100,  100,  100]
! dts    = [ 1,  1,   1,   1,   1,   1]
  call benchmark_mymat(Ns,dts)

 end subroutine     main_benchmark
!-----------------------------------------------------------------------------------------------------------------------------------
 subroutine         init_benchmark
  real(dp)                              :: t0
  character(len_label)                  :: tlabel 
  character(100)                        :: user,host,mach,tdate

  f_mout = 6
  call matrix_random_init

 end subroutine     init_benchmark
!-----------------------------------------------------------------------------------------------------------------------------------
 subroutine         benchmark_mymat(Ns,dts)
  integer , intent(in)                  :: Ns (:)
  real(dp), intent(in)                  :: dts(:)
  integer                               :: i,j, N, iN, itest, ntests
  real(dp)                              :: dt
  integer , allocatable                 :: nops   (:,:)
  real(dp), allocatable                 :: dt_test(:,:)
  !---------------------------------------------------------------------------------------------------------------------------------
  ntests = 2
  !---------------------------------------------------------------------------------------------------------------------------------
  allocate(nops(ntests,size(Ns)),dt_test(ntests,size(Ns)))
  do iN = 1, size(Ns) 
   N = Ns(iN) ; dt = dts(iN); call test_MatmulMyMatrix   (N,dt,nops(1,iN),dt_test(1,iN)); 
  end do

  do iN = 1, size(Ns)
   N = Ns(iN) ; dt = dts(iN); call test_LMatmulArray2    (N,dt,nops(2,iN),dt_test(2,iN)); 
  end do

  write(f_mout,'(A)') '#-----------------------------------------------------------------------------------------------------'
  write(f_mout,'(A)') '#my:       N      MatmulMyMatrix    MatmulArray2'                  !&
  write(f_mout,'(A)') '#-----------------------------------------------------------------------------------------------------'

  do iN = 1, size(Ns)
   write(f_mout,'(A,I8,1000G18.5)') 'my  ',Ns(iN),(dt_test(itest,iN)/nops(itest,iN), itest=1,ntests)
  end do
  !---------------------------------------------------------------------------------------------------------------------------------
 end subroutine     benchmark_mymat
!-----------------------------------------------------------------------------------------------------------------------------------
!Change nsuffle elemets of the array a:
 subroutine         shuffle_array2(a,nsuffle)
  integer, intent(in)                   :: nsuffle
  complex(dp), intent(in out)           :: a(:,:)
  integer                               :: N, i,j,isuffle
  real   (dp)                           :: r

  N           = size(a,1)
  do  isuffle = 1, nsuffle
   call random_number(r); i = int(N * r) + 1
   call random_number(r); j = int(N * r) + 1
   call c_random_number(a(i,j))
  end do

 end subroutine     shuffle_array2
!-----------------------------------------------------------------------------------------------------------------------------------
 subroutine         test_MatmulMyMatrix(N,delt,nops,dt_test)
  integer, intent(in)                   :: N
  real(dp),intent(in)                   :: delt
  real(dp),intent(out)                  :: dt_test
  integer ,intent(out)                  :: nops
  type(MyMatrix)                        :: m1,m2, m3
  real(dp)                              :: ti, tf, dt
  character(len_label)                  :: tlabel

  tlabel = '#t_MatmulMyMatrix:';

  m1   = MyMatrix(N); m2 = MyMatrix(N); m3 = MyMatrix(N)
  call my_random_number(m1 % v)
  call my_random_number(m2 % v)

  ti   = wtime() ; tf = ti; dt = 0.0_dp; nops = 0
  do while( dt < delt)
   call shuffle_array2(m1 % v,5);call shuffle_array2(m2 % v,5); 
   m3  = m1 * m2
   tf  = wtime(); dt = tf - ti; nops = nops + 1
  end do
  dt_test = dt
  write(f_mout,'(A,I8,I20,G28.16)') tlabel,N, nops, dt

 end subroutine     test_MatmulMyMatrix
!-----------------------------------------------------------------------------------------------------------------------------------
 subroutine         test_LMatmulArray2(N,delt,nops,dt_test)
  integer    ,intent(in)                :: N
  real   (dp),intent(in)                :: delt
  real   (dp),intent(out)               :: dt_test
  integer    ,intent(out)               :: nops
  complex(dp)                           :: m1(N,N),m2(N,N), m3(N,N)
  real   (dp)                           :: ti, tf, dt
  character(len_label)                  :: tlabel

  tlabel = '#t_LMatmulArray2:';

  call my_random_number(m1); call my_random_number(m2); 

  ti   = wtime() ; tf = ti; dt = 0.0_dp; nops = 0
  do while( dt < delt)
   call shuffle_array2(m1,5);call shuffle_array2(m2,5);
   m3  = LMATMUL(m1, m2)
   tf  = wtime(); dt = tf - ti; nops = nops + 1
  end do
  dt_test = dt

  write(f_mout,'(A,I8,I20,G28.16)') tlabel,N, nops, dt

 end subroutine     test_LMatmulArray2
!-----------------------------------------------------------------------------------------------------------------------------------
 function           mymatrix_construct(n) result(MAT)
  integer, intent(in)                   :: n
  type(MyMatrix)                        :: MAT

  allocate( MAT % v(n,n) )

  MAT % n = n
  MAT % v = ZERO

 end function       mymatrix_construct
!-----------------------------------------------------------------------------------------------------------------------------------
 function           mymatrix_mult_mymatrix (MATA,MATB) result (MATC)
  class(MyMatrix), intent(in)           :: MATA, MATB
  type (MyMatrix)                       :: MATC

  MATC % n = MATA % n

  allocate( MATC % v( MATA % n , MATA % n) )

  MATC % v = LMATMUL( MATA % v , MATB % v)

 end function       mymatrix_mult_mymatrix
!-----------------------------------------------------------------------------------------------------------------------------------
 function           wtime() result (t)
  real(dp)                              :: t
  integer                               :: count, count_rate, count_max

  call system_clock(count,count_rate,count_max)
  t = dble(count)/dble(count_rate)

 end function       wtime
!-----------------------------------------------------------------------------------------------------------------------------------
 function           array_matmul_array(A,B) result (C)
  complex  (dp), dimension(:,:), intent(in)       :: A,B
  complex  (dp), dimension(size(A,1),size(B,2))   :: C
  complex  (dp)                                   :: alpha, beta
  character(1)                                    :: opA, opB
  integer                                         :: n

  n=size(A,1)
  alpha = ONE ; beta = ZERO
  opA   = 'N' ; opB  = 'N'  ;

  call zgemm(opA ,opB ,n,n,n,alpha,A,n,B,n,beta,C,n)

 end function       array_matmul_array
!-----------------------------------------------------------------------------------------------------------------------------------
 subroutine         my_random_number(m)
  complex(dp), intent(out)                       :: m(:,:)
  real   (dp)                                    :: r(size(m,1),size(m,2)),i(size(m,1),size(m,2))

  call random_number(r) ; call random_number(i)
  m = cmplx(r,i,kind=dp)

 end subroutine     my_random_number
!-----------------------------------------------------------------------------------------------------------------------------------
 subroutine         c_random_number(m)
  complex(dp), intent(out)                       :: m
  real   (dp)                                    :: r,i

  call random_number(r) ; call random_number(i)
  m = cmplx(r,i,kind=dp)

 end subroutine     c_random_number
!-----------------------------------------------------------------------------------------------------------------------------------
 subroutine         matrix_random_init                                                                                              !:.:.:
  integer                                         :: un,nseeds
  integer, allocatable                            :: seeds(:)
  
  call random_seed(size=nseeds)
  ALLOCATE(   seeds    (nseeds))
  open(newunit=un,file="/dev/urandom", access="stream", form="unformatted")
  read(un)    seeds
  seeds=ABS ( seeds)
  call random_seed(put = seeds)

 end subroutine     matrix_random_init
!-----------------------------------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------------------------------
end module          matrix_benchmark
!-----------------------------------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------------------------------
program             matrix_benchmark_program
 use                matrix_benchmark

 call main_benchmark

end  program        matrix_benchmark_program
!-----------------------------------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------------------------------

