OpenMP Tutorial

Serial Pi:

function  funct(dummy) result(dumdum)

    real(kind=8) :: dummy, dumdum

    dumdum = 4.d0 / (1.d0 + dummy * dummy)

end function funct

!

    program compute_omp_pi

!

! compile line:  ifort -integer-size 64 -o pi pi.f90

!

 integer ::  i

 integer(kind = 8)  :: interval

 real(kind=8) ::  width, partial, sum, pi, dummy

real(kind=8) ::  funct

!

! Sum of intervals equals one

!

interval = 800000000000

width = 1.d0 / interval

sum = 0.d0

do i = 1, interval

    partial = width * (i - 0.5d0)

    sum = sum + funct(partial)

enddo

pi = width * sum

print *, "computed pi =", pi

print *, "reference pi = 3.1415926535897932385"

stop

end program compute_omp_pi

PBS Script:

#PBS -l nodes=1:ppn=1,walltime=3:30:00

#PBS -m ae

#PBS -N test_pbs

#PBS -o /N/dc/scratch/hpstrnXX/pi/pi_test.out

#PBS -e /N/dc/scratch/hpstrnXX/pi/pi_test.err

#

cd /N/dc/scratch/hpstrnXX/pi

time ./pi

Note, you must change “XX” to your number.

 

OpenMP  Pi:

function  funct(dummy) result(dumdum)

    real(kind=8) :: dummy, dumdum

    dumdum = 4.d0 / (1.d0 + dummy * dummy)

end function funct

    program compute_omp_pi

!

    USE OMP_LIB

!

! compile line:  ifort -integer-size 64 -openmp -par-num-threads=6 -o omp_pi omp_pi.f90

!

 integer, parameter :: nthreads = 6

 integer ::  i

 integer(kind = 8)  :: interval

 real(kind=8) ::  width, partial, sum, pi, dummy

real(kind=8) ::  funct

!

 call OMP_SET_NUM_THREADS(nthreads)

! Sum of intervals equals one

interval = 800000000000

width = 1.d0 / interval

sum = 0.d0

!$OMP PARALLEL DO PRIVATE(partial), SHARED(width), REDUCTION(+: sum)

do i = 1, interval

    partial = width * (i - 0.5d0)

    sum = sum + funct(partial)

enddo

!! !$OMP END PARALLEL DO

pi = width * sum

print *, "computed pi =", pi

print *, "reference pi = 3.1415926535897932385"

stop

end program compute_omp_pi

Notes:  Notice “reduction.”

 

PBS Script:

#PBS -l nodes=1:ppn=6,walltime=1:00:00

#PBS -m ae

#PBS -N test_omp_pi

#PBS -o /N/dc/scratch/hpstrnXX/threads/omp_test.out

#PBS -e /N/dc/scratch/hpstrnXX/threads/omp_test.err

#

cd /N/dc/scratch/hpstrnXX/threads

time ./omp_pi

Notes:  Same editing of “XX” required here too.

 

PThread example (from Cray User Group):

module pi_module

   use fpthread

   implicit none

   type (fpthread_mutex_t)                      :: reduction_mutex

   type (fpthread_t), dimension(:), allocatable :: tid

   integer :: intervals, num_threads

   real    :: pi, n_1

end module pi_module

 

program compute_pi

   use fpthread

   use pi_module

   implicit none

   integer :: i, ierr

   external PIworker

   external fpthread_mutex_init, fpthread_create, fpthread_join

##   print*, 'How many intervals and threads?'

 ##   read*, intervals, num_threads

   allocate( tid(num_threads) )

   n_1 = 1.0 / real( intervals )

   pi = 0.0

   call fpthread_mutex_init(reduction_mutex, NULL, ierr)

   do i = 1, num_threads

      call fpthread_create(tid(i), NULL, PI_worker, NULL, ierr)

   enddo

   do i = 1, num_threads

      call fpthread_join(tid(i), NULL, ierr)

   enddo

   deallocate( tid )

   print*, 'Computed pi = ', pi

end program compute_pi

 

subroutine PI_worker

   use fpthread

   use pi_module

   implicit none

   type (fpthread_t) :: my_num

   integer           :: i, my_id, ierr

   real              :: sum, my_pi, x, f

   external f

   external fpthread_self, fpthread_mutex_lock, fpthread_mutex_unlock

 

   call fpthread_self( my_num )

   my_id = my_num%thread - minval( tid(:)%thread ) + 1

   sum = 0.0

   do i = myid, intervals, num_threads

      x = n_1 * ( real(i) - 0.5 )

      sum = sum + f(x)

   enddo

   my_pi = n_1 * sum

   call fpthread_mutex_lock(reduction_mutex, ierr)

      pi = pi + my_pi

   call fpthread_mutex_unlock(reduction_mutex, ierr)

end subroutine PI_worker

 

function f(a)

   real :: f, a

   f = 4.0 / (1.0 + a * a)

end function f