C#######################################################################
C
C program matmul to test compiler optimizations and autoparallelization
C
C multiplies a(ardim,mdim)*b(mdim,ardim)
C compile with -DARDIM=<leading a dimension> -DMDIM=<leading b dimension>
C (default for MDIM is ARDIM)
C compile with -DSIZEARGS to provide dimensions as runtime arguments
C compile with -DDGEMM to add call to blas routine dgemm
C
C suggested parameters for ifort:
C ifort -O3 -parallel -mkl -par_threshold90 -DSIZEARGS -DDGEMM \
C   -Wl,-rpath,<mkl lib dir> -o matmulvar.mkl matmul.F
C
C Copyright: Axel Philipp, POEA            last change: 24. April 2014
C
C#######################################################################

#ifndef ARDIM
#define ARDIM 2000
#endif
#ifndef MDIM
#define MDIM ARDIM
#endif

      implicit real*8(a-h,o-z)
      integer ardim,mdim
      real rdtime, tarray(2),dtime
      character*1 trans

#ifdef SIZEARGS

      real*8, allocatable :: a(:,:),b(:,:),c(:,:)
      integer cnt, len, status
      character arg*30
c      
      ardim = ARDIM
      mdim = MDIM
c      
      cnt = command_argument_count ()
c      write (*,*) 'number of command arguments = ', cnt
      if ( cnt .ge. 1 ) then
        call get_command_argument (1, arg, len, status)
        if (status .ne. 0) then
	  write (*,*) 'get_command_argument 1 failed: status = ', status
	  stop
	end if
	read (arg,*) ardim
	mdim = ardim
c	write (*,*) 'command arg 1 = ', ardim	
      end if
      if ( cnt .eq. 2 ) then
        call get_command_argument (2, arg, len, status)
        if (status .ne. 0) then
	  write (*,*) 'get_command_argument 2 failed: status = ', status
	  stop
	end if
	read (arg,*) mdim
c	write (*,*) 'command arg 2 = ', mdim	
      end if
      allocate (a(ardim,mdim),b(mdim,ardim),c(ardim,ardim))

#else
      parameter (ardim = ARDIM)
      parameter (mdim = MDIM)
      dimension a(ardim,mdim),b(mdim,ardim),c(ardim,ardim)
#endif

      write (*,*) 'Running with array sizes', ardim, 'by', mdim
c
      rdtime=dtime(tarray)
      call clockx(stime)
    1 format('  dtime:',f10.3,'   real time:',f10.3,2x,a6)
c
      do j=1,mdim
        dj=dfloat(j)
      do i=1,ardim
        di=dfloat(i)
        a(i,j)=di/dj
      enddo
      enddo
c
c
      do j=1,ardim
        dj=dfloat(j)
      do i=1,mdim
        di=dfloat(i)
        b(i,j)=dj/di
      enddo
      enddo
c
      do j=1,ardim
      do i=1,ardim
        c(i,j)=0.0d0
      enddo
      enddo
c
      rdtime=dtime(tarray)
      call clockx(etime)
      elaptime=(etime-stime)/1.d6
      stime=etime
      write (6,1) rdtime,elaptime,'init'
c      write (6,1) tarray
c
      do i=1,ardim
      do j=1,ardim
      do k=1,mdim
        c(i,j)=c(i,j)+a(i,k)*b(k,j)
      enddo
      enddo
      enddo
c
      rdtime=dtime(tarray)
      call clockx(etime)
      elaptime=(etime-stime)/1.d6
      stime=etime
      write (6,1) rdtime,elaptime,' ijk'
c      write (6,1) tarray
c
      c=0.0d0
      do i=1,ardim
      do k=1,mdim
      do j=1,ardim
        c(i,j)=c(i,j)+a(i,k)*b(k,j)
      enddo
      enddo
      enddo
c
      rdtime=dtime(tarray)
      call clockx(etime)
      elaptime=(etime-stime)/1.d6
      stime=etime
      write (6,1) rdtime,elaptime,' ikj'
c      write (6,1) tarray
c
      c=0.0d0
      do j=1,ardim
      do i=1,ardim
      do k=1,mdim
        c(i,j)=c(i,j)+a(i,k)*b(k,j)
      enddo
      enddo
      enddo
c
      rdtime=dtime(tarray)
      call clockx(etime)
      elaptime=(etime-stime)/1.d6
      stime=etime
      write (6,1) rdtime,elaptime,' jik'
c      write (6,1) tarray
c
      c=0.0d0
      do j=1,ardim
      do k=1,mdim
      do i=1,ardim
        c(i,j)=c(i,j)+a(i,k)*b(k,j)
      enddo
      enddo
      enddo
c
      rdtime=dtime(tarray)
      call clockx(etime)
      elaptime=(etime-stime)/1.d6
      stime=etime
      write (6,1) rdtime,elaptime,' jki'
c      write (6,1) tarray
c
      c=0.0d0
      do k=1,mdim
      do i=1,ardim
      do j=1,ardim
        c(i,j)=c(i,j)+a(i,k)*b(k,j)
      enddo
      enddo
      enddo
c
      rdtime=dtime(tarray)
      call clockx(etime)
      elaptime=(etime-stime)/1.d6
      stime=etime
      write (6,1) rdtime,elaptime,' kij'
c      write (6,1) tarray
c
      c=0.0d0
      do k=1,mdim
      do j=1,ardim
      do i=1,ardim
        c(i,j)=c(i,j)+a(i,k)*b(k,j)
      enddo
      enddo
      enddo
c
      rdtime=dtime(tarray)
      call clockx(etime)
      elaptime=(etime-stime)/1.d6
      stime=etime
      write (6,1) rdtime,elaptime,' kji'
c      write (6,1) tarray
c
      res=0.0
      do j=1,ardim
      do i=1,ardim
        res=res+c(i,j)
      enddo
      enddo
      write (6,3) res
c
      rdtime=dtime(tarray)
c
      c=matmul(a,b)
c
      rdtime=dtime(tarray)
      call clockx(etime)
      elaptime=(etime-stime)/1.d6
      stime=etime
      write (6,1) rdtime,elaptime,'matmul'
c      write (6,1) tarray
c
#ifdef DGEMM
      trans='n'
      one=1.0d0
      zero=0.0d0
      call dgemm(trans,trans,ardim,ardim,mdim,one,a,ardim,b,mdim,
     +  zero,c,ardim)

      rdtime=dtime(tarray)
      call clockx(etime)
      elaptime=(etime-stime)/1.d6
      stime=etime
      write (6,1) rdtime, elaptime,' dgemm'
c
#endif
      res=0.0
      do j=1,ardim
      do i=1,ardim
        res=res+c(i,j)
      enddo
      enddo
c
      write (6,3) res
    3 format(' Sum of elements:',f30.3)
c
      stop
      end
