	subroutine dgemm_compare(nt,NMAX,a,b,c,tm)
        implicit none
        integer NMAX
        integer nt
        real*8 a(NMAX,NMAX),b(NMAX,NMAX),c(NMAX,NMAX)
        real*8 err
        real*8 myrandom
        real*8 tzero,tnow
        real*8 tlapack
        real*8 csecond
        real*8 tmin,tmax,tave,titer
        integer i,j,info,lda,ldb,ldc,nrhs,n,k,ij
        integer iseed
        integer niter
        real*8 tm
        real*8 ONE
      
!  nt = 1 =>   Simple Code
!  nt = 2 =>   MatMUL
!  nt = 3 =>   DGEMM
!  nt <1 || nt > 3 => Error 
        tlapack=0.0
        n=NMAX
        tm=0.0
        ONE=1.00000000000000000
         
!  initialize random numbers generator
        call init_randoms_generator(iseed)

        tmax=-1.0
        tmin=1.00000000e37
!  fill arrays A,B,C with random numbers (in range -100 - 100)
        do i = 1,NMAX
           do j = 1,NMAX
              a(j,i)=myrandom(iseed)
              b(j,i)=myrandom(iseed)
              c(j,i)=0.0
           enddo
        enddo

        nrhs = 1
        lda = NMAX
        ldb = NMAX
        ldc = NMAX

        if (nt. lt. 1 .or. nt.gt.3) then
           tm=-1.0
           return 
        endif

        if ( nt .eq. 1 ) then
        tzero=csecond()
        do i = 1, NMAX
           do j = 1, NMAX
                c(i,j)= 0.00000000000000000
                do ij = 1, NMAX
                  c(i,j)=c(i,j)+a(i,ij)*b(ij,j)
                enddo
           enddo
        enddo
        tnow=csecond()
        tm=tm+tnow-tzero
        return
        endif

        if ( nt.eq. 2) then
        tzero=csecond()
        c=matmul(a,b) 
        tnow=csecond()
        tm=tm+tnow-tzero
        return
        endif

        if (nt.eq.3) then
        tzero=csecond()
        call dgemm('N','N',NMAX,NMAX,NMAX,ONE,a,lda,b,ldb,ONE,c,ldc)
        tnow=csecond()
        tm=tm+tnow-tzero
        return
        endif

        end

