OpenMP 参数扫描并行

Question

我是 OpenMP 新手。我想使用并行 do 循环解决一系列参数值的刚性 ODE 系统。我在下面给出的 Fortran 中使用以下代码。但是，我不知道是否允许在并行 do 循环内调用刚性求解器（作为子程序）？另外，我想在返回主程序之前在子程序中将时间序列数据写入文件名如“r_value_s__value.txt”的文件中。任何人都可以帮忙。下面是代码和错误。我使用gfortran 和标志-fopenmp 来编译。

   PROGRAM OPENMP_PARALLEL_STIFF

           USE omp_lib
           IMPLICIT NONE

           INTEGER :: I, J
           INTEGER, PARAMETER :: RTOT=10, STOT=15
           INTEGER :: TID
           INTEGER, PARAMETER :: NUM_THREADS=8
           DOUBLE PRECISION :: T_INITIAL, T_FINAL
           CALL OMP_SET_NUM_THREADS(NUM_THREADS)
           CALL CPU_TIME(T_INITIAL)
           PRINT*, "TIME INITIAL ",T_INITIAL
!$OMP PARALLEL DO PRIVATE(I,J,TID)

           DO I=1,RTOT
              DO J=1,STOT
                TID=OMP_GET_THREAD_NUM()
                CALL STIFF_DRIVER(TID,I,J,RTOT,STOT)
              END DO
           END DO
!$OMP END PARALLEL DO

          CALL CPU_TIME(T_FINAL)
          PRINT*, "TIME FINAL ",T_FINAL
          PRINT*, "TIME ELAPSED ",(T_FINAL-T_INITIAL)/NUM_THREADS
       END PROGRAM OPENMP_PARALLEL_STIFF

       SUBROUTINE STIFF_DRIVER(TID,II,JJ,RTOT,STOT)

          USE USEFUL_PARAMETERS_N_FUNC

          USE DVODE_F90_M

!     Type declarations:

          IMPLICIT NONE


!     Number of odes for the problem:

          INTEGER :: SERIAL_NUMBER, TID
          INTEGER :: II, JJ, RTOT, STOT, IND
          INTEGER :: J, NTOUT
          INTEGER :: ITASK, ISTATE, ISTATS, I

!     parameters : declaration

          DOUBLE PRECISION, PARAMETER :: s0=0.450D0, dr=1.0D-4, ds=1.0D-2

          DOUBLE PRECISION, DIMENSION(NEQ) :: Y, YOUT
          DOUBLE PRECISION :: ATOL, RTOL, RSTATS, T, TOUT, EPS, TFINAL, DELTAT
          DIMENSION :: RSTATS(22), ISTATS(31)
          DOUBLE PRECISION :: bb, cc, ba, ba1, eta
          CHARACTER(len=45) :: filename

          TYPE (VODE_OPTS) :: OPTIONS

          SERIAL_NUMBER=3011+II+(JJ-1)*RTOT
          IND=TID+3011+II+(JJ-1)*RTOT
          WRITE (*,12)SERIAL_NUMBER,TID
    12    FORMAT ("SL. NO. ",I5," THREAD NO.",I3)

          r=(II-1)*dr
          s=s0+JJ*ds

          EPS = 1.0D-9

!         Open the output file:

          WRITE (filename,93)r,s
    93    FORMAT ("r_",f6.4,"_s_",f4.2,".txt")
          OPEN (UNIT=IND,FILE=filename,STATUS='UNKNOWN',ACTION='WRITE')

!     Parameters for the stiff ODE system

           q0 = 0.60D0;    v = 3.0D0
           Va = 20.0D-4;  Vs = 1.0D-1
           e1 = 1.0D-1;   e2 = 1.10D-5; e3 = 2.3D-3; e4=3.0D-4
          del = 1.7D-4;   mu = 5.9D-4
           al = 1.70D-4;  be = 8.9D-4;  ga = 2.5D-1

!         S and r dependent parameters

           e1s = e1/s;    e2s = e2/(s**2);   e3s = e3/s;    e4s = e4/s
          dels = del*s;    rs = r*s
           e1v = e1/v;     e2v = e2/(v**2);   e3v = e3/v;    e4v = e4/v
          delv = del*v;     rv = r*v

!         SET INITIAL PARAMETERS for INTEGRATION ROUTINES  

          T = 0.0D0
          TFINAL = 200.0D0
          DELTAT = 0.10D0
          NTOUT = INT(TFINAL/DELTAT)
          RTOL = EPS
          ATOL = EPS
          ITASK = 1
          ISTATE = 1

!         Set the initial conditions: USING MODULE USEFUL_PARAMETERS_N_FUNC

          CALL Y_INITIAL(NEQ,Y)

!        Set the VODE_F90 options:

          OPTIONS = SET_OPTS(DENSE_J=.TRUE.,USER_SUPPLIED_JACOBIAN=.FALSE., &
          RELERR=RTOL,ABSERR=ATOL,MXSTEP=100000)

!         Integration:

          DO I=1,NTOUT

          TOUT = (I-1)*DELTAT

          CALL DVODE_F90(F_FUNC,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS)

!         Stop the integration in case of an error

          IF (ISTATE<0) THEN
            WRITE (*,*)"ISTATE ", ISTATE
            STOP
          END IF

!         WRITE DATA TO FILE

          WRITE (IND,*) TOUT,T, Y(NEQ-2)

         END DO


         CLOSE(UNIT=IND)

         RETURN
    END SUBROUTINE STIFF_DRIVER

在文件 openmp_parallel_stiff.f90 的第 ** 行（单位 = 3013） Fortran 运行时错误：文件已在另一个单元中打开

Answer 1

问题在于您选择的格式：f6.4 for r will overflow for r>=10。然后，对于所有线程上的所有r>=10 值，输出将是六个星号******（取决于编译器）。 s 也是如此。

我建议要么限制/检查这些值的范围，要么扩展格式以支持更多数字。

---

正如@francescalus 提到的，另一种可能性是II 和JJ 的组合，其中r 和s 是相同的。

只是为了好玩-让我们算一下：

r=(II-1)*dr
s=s0+JJ*ds

来自r=s 关注

(II-1)*dr = s0+JJ*ds

或

II = 1 + s0/dr + JJ*ds/dr

使用常量s0=0.450D0, dr=1.0D-4, ds=1.0D-2 产生

II = 4501 + JJ*10

因此，只要这种组合同时适用于两个（或更多）线程，您就会遇到观察到的问题。

这种情况的简单解决方案：将线程号添加到文件名中。