module testm
! File: http://ftp.cac.psu.edu/pub/ger/fortran/hdk/OptTest.f90
! by highegg as posted at comp.lang.fortran on 9 Jamiaru 2007.
implicit none

contains

  pure function interpolate(x,xmax,xmin,y) result(yi)
    real,intent(in):: x,xmax,xmin,y(0:)
    real:: yi
    integer:: n,i
    real:: sx
    n = ubound(y,1)
    sx = n*(x-xmin)/(xmax-xmin)
    i = max(min(ceiling(sx),n),1)
    yi = (i-sx)*y(i-1) + (sx-i+1)*y(i)
  end function

end module

program testp
! For y an array of dimension,n, this program illustrates
! that n computatons of 1/y are not optimized out of the
! i-loop. Case 1: 1/y removed from i-loop manually
!         Case 2: 1/y not manually removed from i-loop.
!
  use testm
  real,dimension(:),allocatable:: x,y,y1,z
  integer:: n,i,iter
  real:: xmax,xmin,t_start,t_end
  n = 4000

  allocate(y(n),z(n),x(n),y1(n))

! Case 1
  call random_number(x)
  xmin=minval(x)
  xmax=maxval(x)
  forall(i=1:n) y(i) = sin(xmin+(xmax-xmin)/n*i)
  x = xmin + (xmax-xmin)*x
  z = 0
  call cpu_time(t_start)
  do iter=1,100
    x(1) = x(1) + .001*z(min(iter,n))
    y1 = 1/y  ! manual optimization.
    do i=1,n
      z(i) = interpolate(x(i),xmin,xmax,y1)
    end do
  end do
  call cpu_time(t_end)
  print *,'time 1: ',t_end-t_start

! Case 2
  call random_number(x)
  xmin=minval(x)
  xmax=maxval(x)
  forall(i=1:n) y(i) = sin(xmin+(xmax-xmin)/n*i)
  x = xmin + (xmax-xmin)*x
  z = 0
  call cpu_time(t_start)
  do iter=1,100
    x(1) = x(1) + .001*z(min(iter,n))
    do i=1,n
      z(i) = interpolate(x(i),xmin,xmax,1/y) ! no manual optimization.
    end do
  end do
  call cpu_time(t_end)
  print *,'time 2: ',t_end-t_start

! If compiler does not optimize 1/y out of the i-loop, then
! time 2 (for Case 2) will be many times bigger than time 1 (Case 1).

  deallocate(y,z,x,y1)
end program