! ----------------------------------------------------------------------
! Testing matrix multiplication
! BLAS sgemm/dgemm, f90 matmul, f77 do loops
! Defined operators .BLASgemm., .matmul., .mmloops.
! ----------------------------------------------------------------------
! Command-line compilation:
! ----------------------------------------------------------------------
! Linux
! gfortran -O3 testMM.f90 -lblas (default BLAS library)
! ifort -Ofast -heap-arrays testMM.f90 -mkl (Intel MKL library)
! nvfortran -fast testMM.f90 -lblas (Nvidia HPC SDK library)
! ----------------------------------------------------------------------
! Windows
! (set LIBRARY_PATH=C:\Program Files\PGI\win64\19.10\lib;%LIBRARY_PATH%)
! gfortran -O3 testMM.f90 -llibblas (PGI library)
! ifort -Ofast -heap-arrays -o a testMM.f90 -Qmkl (Intel MKL library)
! (call setvars.bat or ifortvars.bat)
! gfortran -O3 testMM.f90 "%MKLROOT%"\lib\intel64\mkl_rt.lib (Intel MKL library)
! ----------------------------------------------------------------------
! external matmul
! Linux: gfortran -O3 -fexternal-blas testMM.f90 -lblas
!        ifort -Ofast -heap-arrays -qopt-matmul testMM.f90 -mkl
! Win:   gfortran -O3 -fexternal-blas testMM.f90 -llibblas
!        ifort -Ofast -heap-arrays -Qopt-matmul -o a testMM.f90 -Qmkl
! ----------------------------------------------------------------------
! more flags
! gfortran -Ofast -fmax-stack-var-size=0 ...
! gfortran -fcheck=all ..., gfortran -std=f2018 ...
! ifort -C ..., ifort -check all ..., ifort -std18 ...
! nvfortran -C ..., nvfortran -Mbounds ..., nvfortran -Mstandard ...
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Global constants
! ----------------------------------------------------------------------
MODULE mConst

! single and double precision, size of arrays
integer,parameter :: sp=4,dp=8,nmax=5000,mmax=nmax,kmax=nmax
! matrix initialization: 1 constant elements, 2 random numbers
integer,parameter :: iselect=1

! tests to run
logical,parameter,private :: T=.TRUE.,F=.FALSE.
logical,parameter :: runBLAS=T
logical,parameter :: runMatmul=T
logical,parameter :: runLoops=T

END MODULE

! ----------------------------------------------------------------------
! Wallclock timer
! Usage: call tic(); call task(); print *,toc()
! ----------------------------------------------------------------------
MODULE mTicToc

implicit none
private
public tic,toc
integer,parameter :: wp=8
real(wp) :: time

CONTAINS

subroutine tic()
integer count,count_rate
call system_clock(count,count_rate)
time=real(count,wp)/real(count_rate,wp)
end subroutine

real(wp) function toc()
integer count,count_rate
call system_clock(count,count_rate)
toc=real(count,wp)/real(count_rate,wp)-time
end function

END MODULE

! ----------------------------------------------------------------------
! Matrix multiplication via the BLAS library
! Overloaded defined operator .BLASgemm.
! Usage: c4=a4.BLASgemm.b4; c8=a8.BLASgemm.b8
!        for real(4) and real(8) matrices
! ----------------------------------------------------------------------
MODULE mBLAS

implicit none
private
public operator(.BLASgemm.)
integer,parameter :: sp=4,dp=8

interface operator (.BLASgemm.)
module procedure BLASgemm_sp,BLASgemm_dp
end interface

CONTAINS

function BLASgemm_sp(a,b) RESULT (c)
real(sp),intent(in) :: a(:,:),b(:,:)
real(sp) :: c(size(a,1),size(b,2))
real(sp) alpha,beta
integer m,k,n,lda,ldb,ldc
character(1) transa,transb
transa='N' ! no-op 'N', transpose 'T', Hermite conjugate 'H'
transb='N'
lda=size(a,1)
ldb=size(b,1)
ldc=size(c,1)
m=lda
k=ldb
n=size(b,2)
if (size(a,2)/=size(b,1)) print *,'gemm: size(a,2)/=size(b,1)' 
alpha=1.
beta=0.
call sgemm(transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc)
end function

function BLASgemm_dp(a,b) RESULT (c)
real(dp),intent(in) :: a(:,:),b(:,:)
real(dp) :: c(size(a,1),size(b,2))
real(dp) alpha,beta
integer m,k,n,lda,ldb,ldc
character(1) transa,transb
transa='N' ! no-op 'N', transpose 'T', Hermite conjugate 'H'
transb='N'
lda=size(a,1)
ldb=size(b,1)
ldc=size(c,1)
m=lda
k=ldb
n=size(b,2)
if (size(a,2)/=size(b,1)) print *,'gemm: size(a,2)/=size(b,1)' 
alpha=1.
beta=0.
call dgemm(transa,transb,m,n,k,alpha,a,lda,b,ldb,beta,c,ldc)
end function

END MODULE

! ----------------------------------------------------------------------
! Matrix multiplication via the matmul function
! Overloaded defined operator .matmul.
! Usage: c4=a4.matmul.b4; c8=a8.matmul.b8
!        for real(4) and real(8) matrices
! ----------------------------------------------------------------------
MODULE mMatmul

implicit none
private
public operator(.matmul.)
integer,parameter :: sp=4,dp=8

interface operator (.matmul.)
module procedure matmul_sp,matmul_dp
end interface

CONTAINS

function matmul_sp(a,b) RESULT (c)
real(sp),intent(in) :: a(:,:),b(:,:)
real(sp) :: c(size(a,1),size(b,2))
c=matmul(a,b)
end function

function matmul_dp(a,b) RESULT (c)
real(dp),intent(in) :: a(:,:),b(:,:)
real(dp) :: c(size(a,1),size(b,2))
c=matmul(a,b)
end function

END MODULE

! ----------------------------------------------------------------------
! Matrix multiplication via indexed do loops
! Overloaded defined operator .mmloops.
! Usage: c4=a4.mmloops.b4; c8=a8.mmloops.b8
!        for real(4) and real(8) matrices
! ----------------------------------------------------------------------
MODULE mLoops

implicit none
private
public operator(.mmloops.)
integer,parameter :: sp=4,dp=8

interface operator (.mmloops.)
module procedure mmloops_sp,mmloops_dp
end interface

CONTAINS

! C_{MxN} = A_{MxK} * B_{KxN}
function mmloops_sp(a,b) RESULT (c)
real(sp),intent(in) :: a(:,:),b(:,:)
real(sp) :: c(size(a,1),size(b,2))
real(sp) :: t
integer :: mmax,kmax,nmax,m,k,n
mmax=size(a,1); kmax=size(a,2); nmax=size(b,2)
do n=1,nmax
do m=1,mmax
  t=0.
  do k=1,kmax
    t=t+a(m,k)*b(k,n)
  enddo
  c(m,n)=t
enddo
enddo
end function

! C_{MxN} = A_{MxK} * B_{KxN}
function mmloops_dp(a,b) RESULT (c)
real(dp),intent(in) :: a(:,:),b(:,:)
real(dp) :: c(size(a,1),size(b,2))
real(dp) :: t
integer :: mmax,kmax,nmax,m,k,n
mmax=size(a,1); kmax=size(a,2); nmax=size(b,2)
do n=1,nmax
do m=1,mmax
  t=0.
  do k=1,kmax
    t=t+a(m,k)*b(k,n)
  enddo
  c(m,n)=t
enddo
enddo
end function

END MODULE

! ----------------------------------------------------------------------
  PROGRAM TestMM
! ----------------------------------------------------------------------
use mConst
use mTicToc
use mBLAS
use mMatmul
use mLoops
implicit none
real(dp),allocatable :: a(:,:),b(:,:),c(:,:)

ALLOCATE (a(mmax,kmax),b(kmax,nmax),c(mmax,nmax))
select case (iselect)
case (1)
  a=1.
  b=1.
case (2)
  call random_number(a)
  call random_number(b)
end select

if (runBLAS) then
  print *,'BLAS'
  call tic()
  c=real(a,sp).BLASgemm.real(b,sp)
  print *,'dt (SP) = ',toc()
  print '(1x,a,10e20.12)','c =',c(1:1,1:1),sum(c)

  call tic()
  c=a.BLASgemm.b
  print *,'dt (DP) = ',toc()
  print '(1x,a,10e20.12)','c =',c(1:1,1:1),sum(c)
endif

if (runMatmul) then
  print *,'MATMUL'
  call tic()
  c=a.matmul.b
  print *,'dt (DP) = ',toc()
  print '(1x,a,10e20.12)','c =',c(1:1,1:1),sum(c)
endif

if (runLoops) then
  print *,'LOOPS'
  call tic()
  c=a.mmloops.b
  print *,'dt (DP) = ',toc()
  print '(1x,a,10e20.12)','c =',c(1:1,1:1),sum(c)
endif

END PROGRAM