Commits

Randy LeVeque  committed 8bd1c31

OpenMP and MPI examples of pisum for lecture 18

  • Participants
  • Parent commits 8ae36c1

Comments (0)

Files changed (2)

File codes/mpi/pisum1.f90

+! $UWHPSC/codes/mpi/pisum1.f90
+
+! Computes pi using MPI.  
+! Compare to $UWHPSC/codes/openmp/pisum2.f90 
+
+program pisum1
+    use mpi
+    implicit none
+    integer :: ierr, numprocs, proc_num, points_per_proc, n, &
+               i, istart, iend
+    real (kind=8) :: x, dx, pisum, pisum_proc, pi
+
+    call mpi_init(ierr)
+    call mpi_comm_size(MPI_COMM_WORLD, numprocs, ierr)
+    call mpi_comm_rank(MPI_COMM_WORLD, proc_num, ierr)
+
+    ! Ask the user for the number of points
+    if (proc_num == 0) then
+        print *, "Using ",numprocs," processors"
+        print *, "Input n ... "
+        read *, n
+    end if
+
+    ! Broadcast to all procs; everybody gets the value of n from proc 0
+    call mpi_bcast(n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+
+    dx = 1.d0/n
+
+    ! Determine how many points to handle with each proc
+    points_per_proc = (n + numprocs - 1)/numprocs
+    if (proc_num == 0) then   ! Only one proc should print to avoid clutter
+        print *, "points_per_proc = ", points_per_proc
+    end if
+
+    ! Determine start and end index for this proc's points
+    istart = proc_num * points_per_proc + 1
+    iend = min((proc_num + 1)*points_per_proc, n)
+
+    ! Diagnostic: tell the user which points will be handled by which proc
+    print 201, proc_num, istart, iend
+201 format("Process ",i2," will take i = ",i6," through i = ",i6)
+
+    pisum_proc = 0.d0
+    do i=istart,iend
+        x = (i-0.5d0)*dx
+        pisum_proc = pisum_proc + 1.d0 / (1.d0 + x**2)
+        enddo
+
+    call MPI_REDUCE(pisum_proc,pisum,1,MPI_DOUBLE_PRECISION,MPI_SUM,0, &
+                        MPI_COMM_WORLD,ierr)
+
+    if (proc_num == 0) then
+        pi = 4.d0 * dx * pisum 
+        print *, "The approximation to pi is ",pi
+        endif
+
+    call mpi_finalize(ierr)
+
+end program pisum1

File codes/openmp/pisum2.f90

+! $UWHPSC/codes/openmp/pisum2.f90
+! Computes pi using OpenMP.  
+! Compare to $UWHPSC/codes/mpi/pisum1.f90
+
+program pisum2
+    
+    use omp_lib
+    implicit none
+    real(kind=8) :: pi,pisum,pisum_thread,x,dx
+    integer :: n, nthreads, points_per_thread,thread_num
+    integer :: i,istart,iend
+
+    ! Specify number of threads to use:
+    nthreads = 1       ! need this value in serial mode
+    !$ nthreads = 8    
+    !$ call omp_set_num_threads(nthreads)
+    !$ print "('Using OpenMP with ',i3,' threads')", nthreads
+
+    print *, "Input n ... "
+    read *, n
+
+
+    ! First do with parallel do loop ...
+    ! ----------------------------------
+
+    print *, "WITH OMP PARALLEL DO LOOP:"
+    dx = 1.d0 / n
+    pisum = 0.d0
+    !$omp parallel do reduction(+: pisum) &
+    !$omp             private(x)
+    do i=1,n
+        x = (i-0.5d0) * dx
+        pisum = pisum + 1.d0 / (1.d0 + x**2)
+        enddo
+    pi = 4.d0 * dx * pisum
+    print *, 'The approximation to pi = ', pi
+
+
+    ! Now do with parallel blocks of code...
+    ! --------------------------------------
+
+    print *, "WITH OMP PARALLEL:"
+
+    ! Determine how many points to handle with each thread.
+    ! Note that dividing two integers and assigning to an integer will
+    ! round down if the result is not an integer.  
+    ! This, together with the min(...) in the definition of iend below,
+    ! insures that all points will get distributed to some thread.
+    points_per_thread = (n + nthreads - 1) / nthreads
+    print *, "points_per_thread = ",points_per_thread
+
+    dx = 1.d0 / real(n, kind=8)
+    pisum = 0.d0
+
+    !$omp parallel private(i,pisum_thread,x, &
+    !$omp                  istart,iend,thread_num) 
+
+    thread_num = 0     ! needed in serial mode
+    !$ thread_num = omp_get_thread_num()    ! unique for each thread
+
+    ! Determine start and end index for the set of points to be 
+    ! handled by this thread:
+    istart = thread_num * points_per_thread + 1
+    iend = min((thread_num+1) * points_per_thread, n)
+
+    !$omp critical
+    print 201, thread_num, istart, iend
+    !$omp end critical
+201 format("Thread ",i2," will take i = ",i6," through i = ",i6)
+
+    pisum_thread = 0.d0
+    do i=istart,iend
+        x = (i-0.5d0)*dx
+        pisum_thread = pisum_thread + 1.d0 / (1.d0 + x**2)
+        enddo
+
+    ! update global pisum with value from each thread:
+    !$omp critical
+      pisum = pisum + pisum_thread
+    !$omp end critical
+
+    !$omp end parallel 
+
+    pi = 4.d0 * dx * pisum
+
+    print *, 'The approximation to pi = ', pi
+
+end program pisum2
+