CODEWITHSUNDEEP
multithreadingparallel-processingfortranopenmpgfortran
user
user

Reputation: 319

Openmp: use of parallel do with omp_get_thread_num()

I'm splitting a do loop using parallel do and private clause. In this loop I add a variable to itself. Why do I get errors if I don't need a critical block or atomic statement in this case?
How can I fix it?

program trap
use omp_lib 
implicit none
double precision::suma=0.d0 ! sum is a scalar
double precision:: h,x,lima,limb
integer::n,i, istart, iend, thread_num=4, total_threads, ppt
integer(kind=8):: tic, toc, rate
double precision:: time
double precision, dimension(4):: pi= 0.d0

call system_clock(count_rate = rate)
call system_clock(tic)

lima=0.0d0; limb=1.0d0; suma=0.0d0; n=100000000
h=(limb-lima)/n

suma=h*(f(lima)+f(limb))*0.5d0 !first and last points

ppt= n/total_threads
!$ call omp_set_num_threads(total_threads)

!$omp parallel do private (istart, iend, thread_num, i)
  thread_num = omp_get_thread_num()
  !$ istart = thread_num*ppt +1
  !$ iend = min(thread_num*ppt + ppt, n)
do i=istart,iend ! this will control the loop in different threads
  x=lima+i*h
  suma=suma+f(x) 
  pi(thread_num+1)=suma
enddo
!$omp end parallel do 

suma=sum(pi) 
suma=suma*h

print *,"The value of pi is= ",suma ! print once from the first image

call system_clock(toc)
time = real(toc-tic)/real(rate)
print*, 'Time ', time, 's'

contains

double precision function f(y)
double precision:: y
f=4.0d0/(1.0d0+y*y)
end function f

end program trap

I get the following errors:

test.f90:23:35:

23 | thread_num = omp_get_thread_num()

Error: Unexpected assignment statement at (1)

test.f90:24:31:

24 | !$ istart = thread_num*ppt +1

Error: Unexpected assignment statement at (1)

test.f90:25:40:

25 | !$ iend = min(thread_num*ppt + ppt, n)

Error: Unexpected assignment statement at (1)

Compiled with:

gfortran -fopenmp -Wall -Wextra -O2 -Wall -o prog.exe test.f90 

./prog.exe

Upvotes: 1

Views: 381

Answers (2)

Ian Bush
Ian Bush

Reputation: 7434

I don't understand why you are manually splitting up the loop when the worksharing constructs in openmp, such as !$omp do, can do this automatically for you. Below is how I would do it

ian@eris:~/work/stack$ cat thread.f90
program trap
  Use, Intrinsic :: iso_fortran_env, Only :  wp => real64, li => int64
  use omp_lib 
  implicit none
  Real( wp ) ::suma=0.0_wp ! sum is a scalar
  Real( wp ) :: h,x,lima,limb
  integer(li):: tic, toc, rate
  Real( wp ) :: time
  Real( wp ) :: pi
  Integer :: i, n

  call system_clock(count_rate = rate)
  call system_clock(tic)

  lima=0.0_wp; limb=1.0_wp; suma=0.0_wp; n=100000000
  h=(limb-lima)/n

  suma=h*(f(lima)+f(limb))*0.5_wp !first and last points

  pi = 0.0_wp
  !$omp parallel default( None ) private( i, x, lima ) &
  !$omp                          shared( pi, n, h )
  !$omp do reduction( +:pi )
  do i= 1, n
     x  = lima + i * h
     pi = pi + f( x ) 
  enddo
  !$omp end do
  !$omp end parallel

  print *,"The value of pi is= ", pi / n

  call system_clock(toc)
  time = real(toc-tic)/real(rate)
  print*, 'Time ', time, 's on ', omp_get_max_threads(), ' threads'

contains

  function f(y)
    Real( wp ) :: f
    Real( wp ) :: y
    f=4.0_wp/(1.0_wp+y*y)
  end function f

end program trap
ian@eris:~/work/stack$ gfortran --version
GNU Fortran (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

ian@eris:~/work/stack$ export OMP_NUM_THREADS=1
ian@eris:~/work/stack$ ./a.out
 The value of pi is=    3.1415926435902248     
 Time    1.8548842668533325      s on            1  threads
ian@eris:~/work/stack$ export OMP_NUM_THREADS=2
ian@eris:~/work/stack$ ./a.out
 The value of pi is=    3.1415926435902120     
 Time   0.86763000488281250      s on            2  threads
ian@eris:~/work/stack$ export OMP_NUM_THREADS=4
ian@eris:~/work/stack$ ./a.out
 The value of pi is=    3.1415926435898771     
 Time   0.54704123735427856      s on            4  threads
ian@eris:~/work/stack$

Upvotes: 3

dreamcrash
dreamcrash

Reputation: 51393

Instead of

!$omp parallel do private (istart, iend, thread_num, i)
  thread_num = omp_get_thread_num()
  !$ istart = thread_num*ppt +1
  !$ iend = min(thread_num*ppt + ppt, n)

try the following:

!$omp parallel private (istart, iend, thread_num, i)
  thread_num = omp_get_thread_num()
  !$ istart = thread_num*ppt +1
  !$ iend = min(thread_num*ppt + ppt, n)
....
!$omp end parallel

Upvotes: 1

Related Questions