0

我有一个使用 OpenMP 并行化的 fortran 程序。主函数调用一个 Fortran 函数，该函数包含具有常微分方程的系统定义。当跨不同线程拆分时，这运行正常。我还想要链接到具有类似 ODE 定义的 C 文件的选项，当我不使用 OpenMP 时，我可以让它运行。使用 OpenMP 时，ODE 计算不正确。

!************* Main Program OpenMP part **************
!SETUPBV contains a function that contains the ODE
!definition

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(I,IT,NT)

      IT = 0
!$    IT = OMP_GET_THREAD_NUM()
      NT = 1
!$    NT = OMP_GET_NUM_THREADS()
      IF(NLLV>=0.OR.IFST==1) &
        CALL SETUBV(NDIM,NA,NCOL,NINT,NFPR,NRC,NROW,NCLM,                &
         FUNI,ICNI,AP,PAR,NPAR,ICP,A,B,C,DD(1,1,BASE+1),DUPS,            &
         FCFC(1,BASE+1),UPS,UOLDPS,RLOLD,UDOTPS,UPOLDP,DTM,THU,          &
         IFST,IAM,IT,NT,IRF,ICF,IID,NLLV)

      I = (IT*NA+NT-1)/NT+1
      CALL BRBD(A(1,1,I),B(1,1,I),C(1,1,I),D,DD,DUPS(1,(I-1)*NCOL),FAA,FC,&
        FCFC,P0,P1,IFST,IID,NLLV,DET,NDIM,NTST,NA,NBC,NROW,NCLM,         &
        NFPR,NFC,A1,A2,BB,CC,C2,CDBC,                                    &
        SOL,S1,S2,IPR,IPC,IRF(1,I),ICF(1,I),IAM,KWT,IT,NT)

!$OMP END PARALLEL

!************* Extract from Module containing INTERFACE definition **************

SUBROUTINE FUNCINTFC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP) BIND(C,NAME='func')
   USE, INTRINSIC :: ISO_C_BINDING
   INTEGER(C_INT), INTENT(IN) :: NDIM, IJAC, ICP(*)
   REAL(C_DOUBLE), INTENT(IN), VALUE :: U(NDIM), PAR(*)
   REAL(C_DOUBLE), INTENT(OUT), VALUE :: F(NDIM)
   REAL(C_DOUBLE), INTENT(INOUT) :: DFDU(NDIM,NDIM),DFDP(NDIM,*)
END SUBROUTINE FUNCINTFC

!************* Function calling c-function with ODE **************

SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP) 
!--------- ---- 

  USE INTERFACES, ONLY: FUNCINTFC

  IMPLICIT NONE
  INTEGER, INTENT(IN) :: NDIM, IJAC, ICP(*)
  DOUBLE PRECISION, INTENT(IN) :: U(NDIM), PAR(*)
  DOUBLE PRECISION, INTENT(OUT) :: F(NDIM)
  DOUBLE PRECISION, INTENT(INOUT) :: DFDU(NDIM,*), DFDP(NDIM,*)

  CALL FUNCINTFC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP) 

END SUBROUTINE FUNC


!************* C-function containing ODE definition **************

int func(integer *ndim, const doublereal u, const integer *icp,
               const doublereal par, const integer *ijac,
               doublereal f, doublereal *dfdu, doublereal *dfdp) 

{
    /* Builtin functions */
    double exp(doublereal);

    /* Local variables */
    static doublereal e, u1, u2;

/*     ---------- ---- */

/* Evaluates the algebraic equations or ODE right hand side */

/* Input arguments : */
/*      NDIM   :   Dimension of the ODE system */
/*      U      :   State variables */
/*      ICP    :   Array indicating the free parameter(s) */
/*      PAR    :   Equation parameters */

/* Values to be returned : */
/*      F      :   ODE right hand side values */

/* Normally unused Jacobian arguments : IJAC, DFDU, DFDP (see manual) */

    /* Function Body */
    u1 = u[0];
    u2 = u[1];

    e = exp(u2);

    f[0] = -u1 + par[0] * (1 - u1) * e;
    f[1] = -u2 + par[0] * par[1] * (1 - u1) * e - par[2] * u2;

    return 0;
} /* func_ */
4

1 回答 1

0

在 C 中，static局部变量不是线程安全的。static从位置变量的声明中删除。您不需要跨调用保持值。如果你这样做了，你就必须非常小心线程安全。

于 2013-11-08T22:21:53.590 回答