Problem with Fortran + OpenMP

General OpenMP discussion

Problem with Fortran + OpenMP

Postby ftinetti » Tue Nov 15, 2011 1:33 pm

Hi,

I have the following small example reproducing a problem I can't
solve, maybe you figure out what the problem is:

$ cat synth2.f
Code: Select all
     PROGRAM COMPUTE

       use omp_lib

       Implicit none

       REAL vmax, rslt, fpar
       INTEGER nstep

       REAL(KIND = 8) :: start, finish

       vmax = 100.0
       nstep = 100000000


       PRINT *, "Sequential"
       CALL omp_set_num_threads(1)
       start  = OMP_get_wtime()
       rslt = fpar(0.0, vmax, nstep)
       finish = OMP_get_wtime()
       PRINT *, 'rslt = ', rslt
       PRINT *, 'Time: ', finish-start

       PRINT *, "//"

       PRINT *, "Two threads"
       CALL omp_set_num_threads(2)
       start  = OMP_get_wtime()
       rslt = fpar(0.0, vmax, nstep)
       finish = OMP_get_wtime()
       PRINT *, 'rslt = ', rslt
       PRINT *, 'Time: ', finish-start

     END PROGRAM COMPUTE


!-----------------------------------------------------------------------
     Function fpar(xlo, xhi, nstep)
!-----------------------------------------------------------------------
       Implicit none

       REAL fpar, xlo, xhi, step, x, lfpar, lv1, lv2
       INTEGER nstep, i


       step = (xhi-xlo)/float(nstep)

       lv1 = 0.032
       lv2 = 273.0

       lfpar = 0.0

!$OMP PARALLEL DO PRIVATE(i, x)
!$OMP&SHARED(nstep, xlo, step, lv1, lv2)
!$OMP&REDUCTION(+: lfpar)
       do i=1, nstep-1
         x     = xlo + float(i)*step
         lfpar = lfpar + (4.*(lv1/2.0/lv2)**1.5 * exp(-lv1*x*x/2.0/
    &            lv2)*x*x)*step
       enddo
!$OMP END PARALLEL DO

       fpar = lfpar
     End


$ gfortran -Wall -O3 -fopenmp -o synth2 synth2.f
$ ./synth2
Sequential
rslt = 0.25000000
Time: 39.176843863911927
//
Two threads
rslt = 0.32098320
Time: 19.592349765822291


i.e. the result with 2 threads is not the same as that of 1 thread/
sequential. AFAIK I'm using the right OpenMP directive and clauses,
but please let me know if I'm wrong.

I'm using gfortran 4.4.2, which is rather old, but it's hard for me
installing a new one right now, if anyone has a newer gfortran version
(or another Fortran compiler - implementing OpenMP, of course) please
tell me if that fails too.

If everything is fine in the code, I have some possible causes for the
error, or wild guesses, anyway...
1) float() and/or exp() are/is not "thread-safe"... is there some way
to find in the gfortran documentation the "thread-safeness" of Fortran
intrinsic funtions? (I assume the Fortran Standard does not define any
thing about this).
2) Some intermediate value is stored in memory on which there is some
data race.

(cross-posting: comp. lang. Fortran)

Thanks in advance,

Fernando.
ftinetti
 
Posts: 582
Joined: Wed Feb 10, 2010 2:44 pm

Re: Problem with Fortran + OpenMP

Postby ftinetti » Tue Nov 15, 2011 9:01 pm

Hi again,

Never mind, there is a numerical error, nothing related specifically to OpenMP or the compiler used.
ftinetti
 
Posts: 582
Joined: Wed Feb 10, 2010 2:44 pm


Return to Using OpenMP

Who is online

Users browsing this forum: No registered users and 6 guests

cron