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(*,*)
      write(*,*) 'FOR ALL LOOPS:'
      write(*,*)


      CALL CPU_TIME(RTIME(1))

      FORALL (i=1:N_max, j=1:N_max, sum(i,j) /= 0.0 )
          sum(i,j) = sum(i,j) + eps_sum
      END FORALL

      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))

      FORALL (i=1:N_max, j=1:N_max, prod(i,j) /= 0.0 )
        prod(i,j) = prod(i,j)*eps_prod
      END FORALL

      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))

      FORALL (i=1:N_max, j=1:N_max, div(i,j) /= 0.0 )
          div(i,j) = div(i,j)/eps_div
      END FORALL


      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))

      FORALL (i=1:N_max, j=1:N_max, pow(i,j) /= 0.0 )
          pow(i,j) = pow(i,j)**eps_pow
      END FORALL

      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, FOR ALL loops = ', t_tot



      END