OpenMP and ISO_C_BINDING

OpenMP and ISO_C_BINDING

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

 

 

 

 

3 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.
jimdempseyatthecove's picture

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

www.quickthreadprogramming.com

Thank you. That did the trick. :-)

Login to leave a comment.