Page MenuHomec4science

pi_hybrid.f90
No OneTemporary

File Metadata

Created
Sun, Nov 17, 11:44

pi_hybrid.f90

! ==========================================================================
! This exercise is taken from the class Parallel Programming Workshop (MPI,
! OpenMP and Advanced Topics) at HLRS given by Rolf Rabenseifner
! ==========================================================================
program pi_hybrid
use mpi
implicit none
integer, parameter :: n = 10000000
double precision :: omp_t1, omp_t2
double precision :: dx, x, sm, pif
double precision, external :: f
integer :: i
integer :: ierror, prank, psize, error_code, nthreads, nrank
integer :: nlocal, istart, iend
!$ integer, external :: omp_get_max_threads
!$ integer, external :: omp_get_thread_num
integer :: provided
call MPI_init_thread(MPI_THREAD_MULTIPLE, provided, ierror)
if (ierror .ne. MPI_SUCCESS) then
print *, 'Cannot continue with MPI not accepting multi-threaded communications'
call MPI_ABORT(MPI_COMM_WORLD, error_code, ierror)
end if
call MPI_Comm_size(MPI_COMM_WORLD, psize, ierror)
call MPI_Comm_rank(MPI_COMM_WORLD, prank, ierror)
! calculate pi = integral [0..1] 4 / (1 + x**2) dx */
omp_t1 = MPI_Wtime()
nlocal = n / psize
istart = 1 + nlocal * prank
iend = nlocal * (prank + 1);
dx = 1. / n
sm = 0.0;
!$omp parallel private(x)
!$ nthreads = omp_get_max_threads()
!$ nrank = omp_get_thread_num()
!$omp do reduction(+:sm)
do i = istart, iend
x = (1. * i - 0.5) * dx
sm = sm + f(x)
enddo
!$omp end do
!$omp single
call MPI_Allreduce(MPI_IN_PLACE, sm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD,ierror);
!$omp end single nowait
!$omp end parallel
pif = dx * sm
omp_t2 = MPI_Wtime()
if (prank.eq.0) then
print *, 'computed pi =',pif
print *, 'Running time = ', (omp_t2-omp_t1)
endif
call MPI_Finalize(ierror)
end program pi_hybrid
double precision function f(a)
implicit none
double precision, intent(in) :: a
f = 4. / (1. + (a**2))
end function f

Event Timeline