calculate pi by sectional measurement in multi threads.

pi_from_sectional_measurement.f90

program pi_from_sectional_measurement
  !$ USE OMP_LIB

  implicit none

  !$ integer threads_count, thread_number

  ! Const
  real(8), parameter :: pi25 = 3.141592653589793238462643D0

  ! Input data
  integer intervals_count

  ! Benchmark
  real(8) time_start, time_end

  ! Temporary val
  integer i
  real(8) width, x

  ! Answer
  real(8) :: area = 0.0D0
  !$ real(8) :: total_area = 0.0D0

  open(10, FILE = "input_pi2.dat")
  read(10, *) intervals_count
  close(10)
  print *, "Intervals count: ", intervals_count

  width = 1.0D0 / dfloat(intervals_count)

  !$OMP PARALLEL PRIVATE(threads_count, thread_number)
  !$ threads_count = OMP_GET_NUM_THREADS()
  !$ thread_number = OMP_GET_THREAD_NUM()
  !$ write(6, "(A, I6, A, A, I6)") "Thread number: ", thread_number, ", ", "The number of threads: ", threads_count
  !$OMP END PARALLEL

  !$ if (threads_count == 1) then
    call cpu_time(time_start)
  !$ else
  !$ time_start = OMP_GET_WTIME()
  !$ end if

  !$OMP PARALLEL DO REDUCTION(+: total_area) PRIVATE(i, x, area)
  do i = 1, intervals_count
    x = 0.5D0 * width + width * dfloat(i - 1)
    area = area + width * (4.0D0 / (1.0D0 + x * x))

    !$ total_area = area
  end do
  !$OMP END PARALLEL DO

  !$ area = total_area

  !$ if (threads_count == 1) then
    call cpu_time(time_end)
  !$ else
  !$ time_end = OMP_GET_WTIME()
  !$ end if

  write(6, "(A, E15.8)") "pi = ", area
  write(6, "(A, E15.8)") "error = ", dabs(area - pi25)
  write(6, "(A, E15.8)") "Exec time = ", time_end - time_start, "sec"
end program pi_from_sectional_measurement