OpenMP を使用して並列化された Fortran プログラムがあります。main 関数は、常微分方程式を含む系の定義を含む 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_ */