program test_speed
      implicit none

      integer :: i, j
      integer, parameter :: N_max = 6000
      integer, parameter :: DOUBLE = 8        ! compiler depended value for double precision


      real(KIND=DOUBLE) :: rtime(2)
      real(KIND=DOUBLE) :: t_sum, t_prod, t_div, t_pow, t_tot
      real(KIND=DOUBLE) :: sum(1:N_max,1:N_max), eps_sum
      real(KIND=DOUBLE) :: prod(1:N_max,1:N_max), eps_prod
      real(KIND=DOUBLE) :: div(1:N_max,1:N_max), eps_div
      real(KIND=DOUBLE) :: pow(1:N_max,1:N_max),  eps_pow


      eps_sum = 1.0D-07 
      eps_prod = 1.0D+00  + 1.0D-09
      eps_div = 1.0D+00  - 1.0D-07
      eps_pow = 1.0D+00  + 1.0D-20

      sum  = 1.0D+00
      prod = 1.0D+00
      div  = 1.0D+00
      pow  = 1.0D+00  + 1.0D-04

      t_sum  = 0.0D+00
      t_prod = 0.0D+00
      t_div  = 0.0D+00
      t_pow  = 0.0D+00
      t_tot  = 0.0D+00


!     determine the kind numbers accociated with single and double precision

      write(*,*) 'on this compiler the single precision kind is', KIND(0.0)
      write(*,*) 'on this compiler the double precision kind is', KIND(0.0D+00)

      write(*,*) 
      write(*,*) 'DO LOOPS:'
      write(*,*) 

      CALL CPU_TIME(RTIME(1))
          
      IF (RTIME(1) < 0) THEN
        PRINT *,'CPU_TIME NOT WORKING, QUIT'
        PRINT *,'RTIME(1) =',RTIME(1)
        STOP 'UNSUCCESSFUL'
      ENDIF


      CALL CPU_TIME(RTIME(1))

      Do j=1,N_max
        Do i=1,N_max
          IF(sum(i,j) /= 0.0 ) THEN
            sum(i,j) = sum(i,j) + eps_sum
          ENDIF
        ENDDO
      ENDDO

      CALL CPU_TIME(RTIME(2))
      PRINT *, 'ELAPSED CPUTIME TIME, SUM = ',RTIME(2) - RTIME(1)
      t_sum = RTIME(2) - RTIME(1)


      CALL CPU_TIME(RTIME(1)) 

      Do j=1,N_max
        Do i=1,N_max
          IF(prod(i,j) /= 0.0 ) THEN
            prod(i,j) = prod(i,j)*eps_prod
          ENDIF
        ENDDO
      ENDDO

      CALL CPU_TIME(RTIME(2))
      PRINT *, 'ELAPSED CPUTIME TIME, PRODUCT = ',RTIME(2) - RTIME(1)
      t_prod = RTIME(2) - RTIME(1)


      CALL CPU_TIME(RTIME(1))

      Do j=1,N_max
        Do i=1,N_max
          IF(div(i,j) /= 0.0 ) THEN
            div(i,j) = div(i,j)/eps_div
          ENDIF
        ENDDO
      ENDDO
         
      CALL CPU_TIME(RTIME(2))
      PRINT *, 'ELAPSED CPUTIME TIME, DIVISION = ',RTIME(2) - RTIME(1)
      t_div = RTIME(2) - RTIME(1)


      CALL CPU_TIME(RTIME(1))
        
      Do j=1,N_max
        Do i=1,N_max
          IF(pow(i,j) /= 0.0 ) THEN
            pow(i,j) = pow(i,j)**eps_pow 
          ENDIF
        ENDDO
      ENDDO

      CALL CPU_TIME(RTIME(2))
      PRINT *, 'ELAPSED CPUTIME TIME, POWER = ',RTIME(2) - RTIME(1)
      t_pow = RTIME(2) - RTIME(1)


      t_tot = t_sum + t_prod + t_div + t_pow
         
      write(*,*) 'Total CPU time , do loops = ', t_tot


      END