Неправильный результат при параллельном выполнении кода

Question

Неправильный результат при параллельном выполнении кода

gfortran Компилятор дает неправильный ответ, когда я запускаю параллельную программу с использованием OpenMP. В то же время, ifort дает точный результат.

Это весь скомпилированный код.

!_______________________________________________________________ !
!____________MODULE SECTION_____________________________________ !

  MODULE MATRIC
    IMPLICIT NONE
    INTEGER , PARAMETER :: NG = 40  
    DOUBLE PRECISION,SAVE :: Z , PA , PB ,CMU 
    DOUBLE PRECISION , PARAMETER :: PI=2.0D0*ACOS(0.0D0) , &
             FPI=4.0D0*PI , SQFPI = SQRT(FPI), DLAM=1.0D0
    DOUBLE PRECISION , DIMENSION(450) :: DEL1,  DEL2, X,  R ,  SNLO 
    DOUBLE PRECISION :: XG(60) , WG(60) 
  END MODULE MATRIC
!_________________________________________________________________________!
!                  MODULE SECTION 
!__________________________________________________________________________!

  MODULE POTDATA
    IMPLICIT NONE
    INTEGER                            :: IA , IB , ID       
    DOUBLE PRECISION                   :: RA , RB , R1s(450)     
  END MODULE POTDATA
!__________________________________________________________________________!



!______________________________________________________________________!

  program check
    use matric
    use potdata
    implicit double precision(a-h,o-z)

    pa   = 0.72D0  ;  pb   =  0.19D0  
    mesh = 441     ;  noint=  40      ;  z   =  2.0d0    
    CALL GAULEG(-1.d0,1.d0)

    NB = MESH/NOINT
    I = 1
    X(I) = 0.0D+00
    DELTAX = 0.0025D+00*40.0D+00/DBLE(NOINT)
    DO  J=1,NB
      IMK = (J-1)*NOINT + 1
      DO K=1,NOINT
        AK=K
        I=I+1
        X(I)=X(IMK)+AK*DELTAX
      END DO
      DELTAX=2.0D+00*DELTAX
    END DO

    CMU=(9.0D00*PI*PI/(128.0D00*Z))**THIRD

    R(1)=0.0D+00 ;  SNLO(1) = 0.D00
    DO  I=2,MESH
      R(I)=CMU*X(I)
      SNLO(I) = R(I)*dexp(-Z*R(I))
      R1S(I) = SNLO(I)/(SQFPI*R(I))
    END DO

    call EFFPOT(MESH,NOINT)

  end program check


  subroutine EFFPOT(MESH,NOINT)
    USE OMP_LIB
    USE MATRIC  
    USE POTDATA 
    implicit none 
    integer, intent(in) :: MESH, NOINT 
    double precision::anorm(450)
    double precision, external :: funct
    double precision :: asum, fac, cnorm

!$omp parallel do default(none) private(del1,ia,asum,ib,ra,rb,fac) &
!$omp shared(id,mesh,r,anorm,NOINT,del2,R1s)
    do  ia = 2,mesh
      ra = r(ia)
      if(R1s(ia).lt.1.D-7.and.R1s(ia).ge.1.D-8)id = ia
      do ib = 2,mesh
         rb = r(ib)
         call QGAUSS(funct,-1.d0,1.d0,fac)
         del1(ib) = rb**2*fac*R1s(ib)**2
      end do
      CALL NCDF(del1,ASUM,r(2),mesh,NOINT)
      anorm(ia) = 2.0d0*pi*asum
      del2(ia)  = 2.0d0*pi*asum*(ra*R1s(ia))**2
    end do
!$omp end parallel do

    CALL NCDF(del2,ASUM,r(2),mesh,NOINT)
    cnorm = 1.0/dsqrt(4.*pi*ASUM)
    write(6,*)'cnorm =',cnorm

    return 
  end


  double precision function funct(x)

    USE POTDATA , ONLY : RA , RB 
    USE MATRIC  , ONLY : PA , PB  , DLAM

    implicit none
    double precision, intent(in) :: x
    double precision             :: f1, f2, ramrb

    ramrb = dsqrt(ra**2+rb**2-2.d0*ra*rb*x)
    f1 = dcosh(pa*ra)+dcosh(pa*rb)
    f2  = 1.d0+0.5*dlam*ramrb*dexp(-pb*ramrb)
    funct = (f1*f2)**2
    return
  end


  SUBROUTINE QGAUSS(func,aa,bb,ss)
    USE OMP_LIB
    USE MATRIC , ONLY : XG ,WG , NG 
    IMPLICIT DOUBLE PRECISION(A-H,O-Z)
    external func
    xm = 0.5d0*(bb+aa)
    xl = 0.5d0*(bb-aa)
    ss = 0.d0
    do  j=1,ng
      dx = xl*xg(j)
      ss = ss + wg(j)*(func(xm+dx)+func(xm-dx))
    end do
    ss = xl*ss/2.0
    return
  END


  SUBROUTINE GAULEG(x1,x2)

    USE MATRIC , ONLY : XG ,WG ,NG , PI

    IMPLICIT DOUBLE PRECISION(A-H,O-Z)
    eps = 1.d-14
    m = (ng+1)/2
    xm = 0.5D0*(x1+x2)
    xl = 0.5D0*(x2-x1)

    do i=1,m
      z = dcos(pi*(dfloat(i)-0.25d0)/(dfloat(ng)+0.5d0))
1     continue
      p1 = 1.d0
      p2 = 0.d0

      do j=1,ng
        p3 = p2
        p2 = p1
        p1 = ((2.d0*dfloat(j)-1.d0)*z*p2  &
          - (dfloat(j)-1.d0)*p3)/dfloat(j)
      end do

      pp = dfloat(ng)*(z*p1-p2)/(z*z-1.d0)
      z1 = z
      z = z1 - p1/pp
      if (dabs(z-z1).gt.eps) go to 1
      xg(i) = xm - xl*z
      xg(ng+1-i) = xm + xl*z
      wg(i) = 2.d0*xl/((1.d0-z*z)*pp*pp)
      wg(ng+1-i) = wg(i)                          
    end do
    return
  end


  SUBROUTINE NCDF(F,ASUM,H,KKK,NOINT)
    IMPLICIT DOUBLE PRECISION (A-H,O-Z)
    DIMENSION F(450)
    NBLOCK = (KKK-2)/NOINT + 1
    C2HO45 = 2.0D+00*H/45.0D+00      
    ASUM = 0.0D+00

    DO  J=1,NBLOCK
      ISTAR = NOINT*(J-1)+5
      IEND = NOINT*J + 1
      IEND = MIN0(KKK,IEND)
      DO  I=ISTAR,IEND,4
          ASUM = ASUM + C2HO45*(7.0D+00*(F(I-4)+F(I))  &
                +32.0D+00*(F(I-3)+F(I-1)) + 12.0D+00*F(I-2))
      END DO
      IF(IEND.EQ.KKK) GO TO 4
      C2HO45 = 2.0D+00*C2HO45
4   END DO

    RETURN
  END

Спасибо всем особенно @Vladimir, кто заинтересовался моей проблемой. Наконец, я получил правильный ответ, удалив ra и rb из данных модуля и определив функцию как funct(x, ra, rb), а затем удалив ra и rb из цикла. Поскольку я писал ra, rb, то читал их значения в приведенном выше коде, поэтому цикл имел зависимость от потока. Теперь я получаю точный результат от обоих компиляторов (который 8.7933767516) параллельно, последовательно оба. Точный путь это

subroutine EFFPOT(MESH,NOINT)
    USE OMP_LIB
    USE MATRIC  
    USE POTDATA 
  implicit none 
  integer, intent(in) :: MESH, NOINT 
  double precision::anorm(450)
  double precision, external :: funct
  double precision :: asum, fac, cnorm
 !$omp parallel do default(none) private(del1,ia,asum,ib,fac) &
 !$omp shared(id,mesh,r,anorm,NOINT,del2,R1s)

  do  ia = 2,mesh
      if(R1s(ia).lt.1.D-7.and.R1s(ia).ge.1.D-8)id = ia
      do ib = 2,mesh
         call QGAUSS(funct,-1.d0,1.d0,fac,r(ia),r(ib))
         del1(ib) = r(ib)**2*fac*R1s(ib)**2
      end do
      CALL NCDF(del1,ASUM,r(2),mesh,NOINT)
      anorm(ia) = 2.0d0*pi*asum
      del2(ia)  = 2.0d0*pi*asum*(r(ia)*R1s(ia))**2
  end do

 !$omp end parallel do
  CALL NCDF(del2,ASUM,r(2),mesh,NOINT)
  cnorm = 1.0/dsqrt(4.*pi*ASUM)
  write(6,*)'cnorm =',cnorm

  return 
  end


      double precision function funct(x,ra,rb)
      USE MATRIC  , ONLY : PA , PB  , DLAM

      implicit none
      double precision, intent(in) :: x, ra, rb
      double precision             :: f1, f2, ramrb

      ramrb = dsqrt(ra**2+rb**2-2.d0*ra*rb*x)
      f1 = dcosh(pa*ra)+dcosh(pa*rb)
      f2  = 1.d0+0.5*dlam*ramrb*dexp(-pb*ramrb)
      funct = (f1*f2)**2
  return
  end
  SUBROUTINE QGAUSS(func,aa,bb,ss,ra,rb)
     USE OMP_LIB
     USE MATRIC , ONLY : XG ,WG , NG 
     IMPLICIT DOUBLE PRECISION(A-H,O-Z)
     external func
     xm = 0.5d0*(bb+aa)
     xl = 0.5d0*(bb-aa)
     ss = 0.d0
   do  j=1,ng
     dx = xl*xg(j)
     ss = ss + wg(j)*(func(xm+dx,ra,rb)+func(xm-dx,ra,rb))
   end do
   ss = xl*ss/2.0
   return
  END

2

fortran openmp gfortran hpc

Источник

user4943203 05 июн '15 в 14:32

2 ответа

Другие вопросы по тегам fortran openmp gfortran hpc

user1374437 08 июн '15 в 11:41 2015-06-08 11:41 · Answer 1 · 2015-06-08 11:41

Причиной вашей проблемы является то, что стандарт OpenMP не определяет, что произойдет, если private Доступ к элементу списка осуществляется в регионе, но за пределами конструкции. См пример private.2f (см. стр. 135 стандартного дополнения OpenMP) для краткой версии той же проблемы.

В частности, переменные модуля ra а также rb объявлены private в параллельной области OpenMP внутри EFFPOT а также доступ из funct, funct находится в динамической области действия параллельной области, но (лексически) вне ее, и поэтому не указано, ra а также rb ссылается на funct являются исходными переменными модуля или их частными копиями (большинство компиляторов используют исходные переменные).

Вы уже нашли одно из решений. Другой будет объявить ra а также rbthreadprivate так как они используются только для передачи данных из EFFPOT в funct:

MODULE POTDATA
  IMPLICIT NONE
  INTEGER                            :: IA , IB , ID       
  DOUBLE PRECISION                   :: RA , RB , R1s(450)
  !$OMP THREADPRIVATE(RA,RB)
END MODULE POTDATA

Затем вы должны также удалить ra а также rb из списка private положение о параллельной области внутри EFFPOT,

На некоторых платформах, например OS X, используя threadprivate с GCC (т.е. gfortran) может быть медленнее, чем фактическая передача двух переменных в качестве аргументов из-за эмулируемого TLS.

Обратите внимание, что эту семантическую ошибку действительно трудно обнаружить, и многие инструменты OpenMP на самом деле ее не обнаружат.

user4271479 05 июн '15 в 14:52 2015-06-05 14:52 · Answer 2 · 2015-06-05 14:52

Прежде всего, очень сложно сказать что-то конкретное, не видя реального кода. Тем не менее, у меня есть некоторые комментарии по поводу вашей ситуации и выводов, которые вы делаете.

Тот факт, что ваша программа работает нормально как при параллельном, так и последовательном выполнении при компиляции с помощью "ifort", не означает, что ваша программа правильная. Так как ошибки компилятора, приводящие к программам, дающим неправильные ответы, очень редки, но, с другой стороны, параллельное программирование вручную очень подвержено ошибкам, мы должны предположить проблему с тем, как вы распараллеливали свой код. Мы, наверное, говорим о состоянии гонки.

И нет, проблема, с которой вы столкнулись, не является чем-то необычным. Когда у вас есть условие гонки, часто случается, что последовательное выполнение работает везде, и ваше параллельное выполнение работает в одних средах и дает сбой в других. Даже часто бывает, что ваш код дает разные ответы каждый раз, когда вы его вызываете (не только в зависимости от компилятора, но и от многих других факторов, которые могут со временем меняться).

Я предлагаю вам установить параллельный отладчик, такой как, например, TotalView, который поможет вам отслеживать различные потоки и их состояния. Попробуйте найти простую среду тестирования (как можно меньше потоков), которая надежно завершится неудачей.