OpenMP and ISO_C_BINDING

[Ed_M_]

I have a fortran program that is parallelised using OpenMP. The main function calls a Fortran function that contains the definition of a system with Ordinary Differential Equations. This runs okay when split across different threads. I also want the option to link in a C file with similar definition for the ODE's, and I can get this to run when I don't use OpenMP. The ODE is calculated incorrectly though when using OpenMP.

!************* 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_ */

[jimdempseyatthecove]

Your C function has:

static doublereal e, u1, u2;

Remove the "static". The consequences of "static" inside the scope of a function (or other scope) is somewhat equivilent to Fortran's SAVE attribute. IOW one copy is stored in a static area that is not on stack, meaning the last value on exit of scope is valid on entry of scope. Think of this as shared(e,u1,u2).

Jim Dempsey

[Ed_M_]

Thank you. That did the trick. :-)