Skip to content
Snippets Groups Projects
Commit f47a7efb authored by RODIER Quentin's avatar RODIER Quentin
Browse files

Quentin 26/09/2024: use gamma version in PURE FUNCTIONS with array-syntax for GPU

parent 227eb179
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/MNH/gamma.f90
+ 61
− 50
View file @ f47a7efb
!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier !MNH_LIC Copyright 1995-2022 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1. !MNH_LIC for details. version 1.
!########################
!
!-------------------------------------------------------------------------- !--------------------------------------------------------------------------
! !
! !
!* 1. FUNCTION GAMMA FOR SCALAR VARIABLE !* 1. FUNCTION GAMMA FOR SCALAR VARIABLE
! !
! !
! ###################################### ! ###########################################
FUNCTION GAMMA_X0D(PX) RESULT(PGAMMA) PURE FUNCTION GAMMA_X0D(PX) RESULT(PGAMMA)
USE PARKIND1, ONLY : JPRB ! ###########################################
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
! ######################################
! !
! !
!!**** *GAMMA * - Gamma function !!**** *GAMMA * - Gamma function
...@@ -49,11 +45,19 @@ ...@@ -49,11 +45,19 @@
!! ------------- !! -------------
!! Original 7/11/95 !! Original 7/11/95
!! C. Barthe 9/11/09 add a function for 1D arguments !! C. Barthe 9/11/09 add a function for 1D arguments
! P. Wautelet 22/06/2022: GAMMA_X0D is now declared PURE
! !
!* 0. DECLARATIONS !* 0. DECLARATIONS
! ------------ ! ------------
! !
#if defined(MNH_BITREP) || defined(MNH_BITREP_OMP)
USE MODI_BITREP
#endif
!
IMPLICIT NONE IMPLICIT NONE
!$acc routine seq
! !
!* 0.1 declarations of arguments and result !* 0.1 declarations of arguments and result
! !
...@@ -66,8 +70,6 @@ INTEGER :: JJ ! Loop index ...@@ -66,8 +70,6 @@ INTEGER :: JJ ! Loop index
REAL :: ZSER,ZSTP,ZTMP,ZX,ZY,ZCOEF(6) REAL :: ZSER,ZSTP,ZTMP,ZX,ZY,ZCOEF(6)
REAL :: ZPI REAL :: ZPI
! !
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('GAMMA_X0D',0,ZHOOK_HANDLE)
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
!* 1. SOME CONSTANTS !* 1. SOME CONSTANTS
...@@ -90,12 +92,16 @@ ZPI = 3.141592654 ...@@ -90,12 +92,16 @@ ZPI = 3.141592654
! !
IF (PX .LT. 0.) THEN IF (PX .LT. 0.) THEN
ZX = 1. - PX ZX = 1. - PX
ELSE ELSE
ZX = PX ZX = PX
END IF END IF
ZY = ZX ZY = ZX
ZTMP = ZX + 5.5 ZTMP = ZX + 5.5
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
ZTMP = (ZX + 0.5) * ALOG(ZTMP) - ZTMP ZTMP = (ZX + 0.5) * ALOG(ZTMP) - ZTMP
#else
ZTMP = (ZX + 0.5) * BR_LOG(ZTMP) - ZTMP
#endif
ZSER = 1.000000000190015 ZSER = 1.000000000190015
! !
DO JJ = 1, 6 DO JJ = 1, 6
...@@ -104,11 +110,18 @@ DO JJ = 1, 6 ...@@ -104,11 +110,18 @@ DO JJ = 1, 6
END DO END DO
! !
IF (PX .LT. 0.) THEN IF (PX .LT. 0.) THEN
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
PGAMMA = ZPI / SIN(ZPI*PX) / EXP(ZTMP + ALOG(ZSTP*ZSER/ZX)) PGAMMA = ZPI / SIN(ZPI*PX) / EXP(ZTMP + ALOG(ZSTP*ZSER/ZX))
#else
PGAMMA = ZPI / SIN(ZPI*PX) / BR_EXP(ZTMP + BR_LOG(ZSTP*ZSER/ZX))
#endif
ELSE ELSE
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
PGAMMA = EXP(ZTMP + ALOG(ZSTP*ZSER/ZX)) PGAMMA = EXP(ZTMP + ALOG(ZSTP*ZSER/ZX))
#else
PGAMMA = BR_EXP(ZTMP + BR_LOG(ZSTP*ZSER/ZX))
#endif
END IF END IF
IF (LHOOK) CALL DR_HOOK('GAMMA_X0D',1,ZHOOK_HANDLE)
RETURN RETURN
! !
END FUNCTION GAMMA_X0D END FUNCTION GAMMA_X0D
...@@ -119,11 +132,9 @@ END FUNCTION GAMMA_X0D ...@@ -119,11 +132,9 @@ END FUNCTION GAMMA_X0D
!* 1. FUNCTION GAMMA FOR 1D ARRAY !* 1. FUNCTION GAMMA FOR 1D ARRAY
! !
! !
! ###################################### ! ###########################################
FUNCTION GAMMA_X1D(PX) RESULT(PGAMMA) PURE FUNCTION GAMMA_X1D(PX) RESULT(PGAMMA)
USE PARKIND1, ONLY : JPRB ! ###########################################
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
! ######################################
! !
! !
!!**** *GAMMA * - Gamma function !!**** *GAMMA * - Gamma function
...@@ -157,12 +168,17 @@ END FUNCTION GAMMA_X0D ...@@ -157,12 +168,17 @@ END FUNCTION GAMMA_X0D
!! MODIFICATIONS !! MODIFICATIONS
!! ------------- !! -------------
!! Original 7/11/95 !! Original 7/11/95
! P. Wautelet 22/06/2022: GAMMA_X1D is now declared PURE
!! !!
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
!* 0. DECLARATIONS !* 0. DECLARATIONS
! ------------ ! ------------
! !
#if defined(MNH_BITREP) || defined(MNH_BITREP_OMP)
USE MODI_BITREP
#endif
!
IMPLICIT NONE IMPLICIT NONE
! !
!* 0.1 declarations of arguments and result !* 0.1 declarations of arguments and result
...@@ -173,12 +189,10 @@ REAL, DIMENSION(SIZE(PX)) :: PGAMMA ...@@ -173,12 +189,10 @@ REAL, DIMENSION(SIZE(PX)) :: PGAMMA
!* 0.2 declarations of local variables !* 0.2 declarations of local variables
! !
INTEGER :: JJ ! Loop index INTEGER :: JJ ! Loop index
INTEGER :: JI ! Loop index REAL, DIMENSION(SIZE(PX)) :: ZSER,ZSTP,ZTMP,ZX,ZY
REAL :: ZSER, ZSTP, ZTMP, ZX, ZY, ZCOEF(6) REAL :: ZCOEF(6)
REAL :: ZPI REAL :: ZPI
! !
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('GAMMA_X1D',0,ZHOOK_HANDLE)
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
!* 1. SOME CONSTANTS !* 1. SOME CONSTANTS
...@@ -193,35 +207,32 @@ ZCOEF(6) = -0.5395239384953E-5 ...@@ -193,35 +207,32 @@ ZCOEF(6) = -0.5395239384953E-5
ZSTP = 2.5066282746310005 ZSTP = 2.5066282746310005
! !
ZPI = 3.141592654 ZPI = 3.141592654
! ZX(:) = PX(:)
!------------------------------------------------------------------------------- WHERE ( PX(:)<0.0 )
! ZX(:) = 1.- PX(:)
!* 2. COMPUTE GAMMA END WHERE
! ------------- ZY(:) = ZX(:)
! ZTMP(:) = ZX(:) + 5.5
DO JI = 1, SIZE(PX) #if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
IF (PX(JI) .LT. 0.) THEN ZTMP(:) = (ZX(:) + 0.5)*ALOG(ZTMP(:)) - ZTMP(:)
ZX = 1. - PX(JI) #else
ELSE ZTMP(:) = (ZX(:) + 0.5)*BR_LOG(ZTMP(:)) - ZTMP(:)
ZX = PX(JI) #endif
END IF ZSER(:) = 1.000000000190015
ZY = ZX !
ZTMP = ZX + 5.5 DO JJ = 1 , 6
ZTMP = (ZX + 0.5) * ALOG(ZTMP) - ZTMP ZY(:) = ZY(:) + 1.0
ZSER = 1.000000000190015 ZSER(:) = ZSER(:) + ZCOEF(JJ)/ZY(:)
!
DO JJ = 1, 6
ZY = ZY + 1.0
ZSER = ZSER + ZCOEF(JJ) / ZY
END DO
!
IF (PX(JI) .LT. 0.) THEN
PGAMMA = ZPI / SIN(ZPI*PX(JI)) / EXP(ZTMP + ALOG(ZSTP*ZSER/ZX))
ELSE
PGAMMA = EXP(ZTMP + ALOG(ZSTP*ZSER/ZX))
END IF
END DO END DO
IF (LHOOK) CALL DR_HOOK('GAMMA_X1D',1,ZHOOK_HANDLE) !
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
PGAMMA(:) = EXP( ZTMP(:) + ALOG( ZSTP*ZSER(:)/ZX(:) ) )
#else
PGAMMA(:) = BR_EXP( ZTMP(:) + BR_LOG( ZSTP*ZSER(:)/ZX(:) ) )
#endif
WHERE ( PX(:)<0.0 )
PGAMMA(:) = ZPI/SIN(ZPI*PX(:))/PGAMMA(:)
END WHERE
RETURN RETURN
! !
END FUNCTION GAMMA_X1D END FUNCTION GAMMA_X1D
src/PHYEX/aux/modi_gamma.f90
+ 2
− 2
View file @ f47a7efb
...@@ -8,12 +8,12 @@ ...@@ -8,12 +8,12 @@
! !
INTERFACE GAMMA INTERFACE GAMMA
! !
FUNCTION GAMMA_X0D(PX) RESULT(PGAMMA) PURE FUNCTION GAMMA_X0D(PX) RESULT(PGAMMA)
REAL, INTENT(IN) :: PX REAL, INTENT(IN) :: PX
REAL :: PGAMMA REAL :: PGAMMA
END FUNCTION GAMMA_X0D END FUNCTION GAMMA_X0D
! !
FUNCTION GAMMA_X1D(PX) RESULT(PGAMMA) PURE FUNCTION GAMMA_X1D(PX) RESULT(PGAMMA)
REAL, DIMENSION(:), INTENT(IN) :: PX REAL, DIMENSION(:), INTENT(IN) :: PX
REAL, DIMENSION(SIZE(PX)) :: PGAMMA REAL, DIMENSION(SIZE(PX)) :: PGAMMA
END FUNCTION GAMMA_X1D END FUNCTION GAMMA_X1D
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Documentation | Mentions légales | CGU | Accessibilité : non conforme