MODULE kinds INTEGER, PARAMETER :: dp = KIND(0.0D0) END MODULE MODULE T_C_G0 USE kinds, ONLY: dp IMPLICIT NONE PUBLIC :: T_C_G0_n REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE, SAVE :: C0 INTEGER, PARAMETER :: degree=13 REAL(KIND=dp), PARAMETER :: target_error= 0.100000E-08 INTEGER, PARAMETER :: nderiv_max=21 INTEGER :: nderiv_init=-1 INTEGER, SAVE :: patches=-1 CONTAINS SUBROUTINE T_C_G0_n(RES,use_gamma,R,T,NDERIV) IMPLICIT NONE REAL(KIND=dp), INTENT(OUT) :: RES(*) LOGICAL, INTENT(OUT) :: use_gamma REAL(KIND=dp),INTENT(IN) :: R,T INTEGER, INTENT(IN) :: NDERIV REAL(KIND=dp) :: upper,lower,X1,X2,TG1,TG2 use_gamma=.FALSE. upper=R**2 + 11.0_dp*R + 50.0_dp lower=R**2 - 11.0_dp*R + 0.0_dp IF (T>upper) THEN RES(1:NDERIV+1)=0.0_dp RETURN ENDIF IF (R>11.0_dp) THEN IF (Tnderiv_max) STOP "T_C_G0 init failed" nderiv_init=Nder IF(ALLOCATED(C0)) DEALLOCATE(C0) ! round up to a multiple of 32 to give some generous alignment for each C0 ALLOCATE(C0(32*((31+(Nder+1)*(degree+1)*(degree+2)/2)/32),patches)) ALLOCATE(chunk((nderiv_max+1)*(degree+1)*(degree+2)/2)) C0=0 DEALLOCATE(chunk) END SUBROUTINE INIT SUBROUTINE FREE() IF(ALLOCATED(C0)) DEALLOCATE(C0) nderiv_init=-1 END SUBROUTINE FREE SUBROUTINE PD2VAL(RES,NDERIV,TG1,TG2,C0) IMPLICIT NONE REAL(KIND=dp), INTENT(OUT) :: res(*) INTEGER, INTENT(IN) :: NDERIV REAL(KIND=dp),INTENT(IN) :: TG1,TG2 REAL(KIND=dp) :: T1(0:13) REAL(KIND=dp) :: T2(0:13) REAL(KIND=dp), INTENT(IN) :: C0(105,*) INTEGER :: I,J,K REAL(KIND=dp), PARAMETER :: SQRT2=1.4142135623730950488016887242096980785696718753_dp T1(0)=1.0_dp T2(0)=1.0_dp T1(1)=SQRT2 * TG1 T2(1)=SQRT2 * TG2 T1(2)= 2 * TG1 * T1(1) - SQRT2 T2(2)= 2 * TG2 * T2(1) - SQRT2 T1(3) = 2 * TG1 * T1(2) - T1(1) T2(3) = 2 * TG2 * T2(2) - T2(1) T1(4) = 2 * TG1 * T1(3) - T1(2) T2(4) = 2 * TG2 * T2(3) - T2(2) T1(5) = 2 * TG1 * T1(4) - T1(3) T2(5) = 2 * TG2 * T2(4) - T2(3) T1(6) = 2 * TG1 * T1(5) - T1(4) T2(6) = 2 * TG2 * T2(5) - T2(4) T1(7) = 2 * TG1 * T1(6) - T1(5) T2(7) = 2 * TG2 * T2(6) - T2(5) T1(8) = 2 * TG1 * T1(7) - T1(6) T2(8) = 2 * TG2 * T2(7) - T2(6) T1(9) = 2 * TG1 * T1(8) - T1(7) T2(9) = 2 * TG2 * T2(8) - T2(7) T1(10) = 2 * TG1 * T1(9) - T1(8) T2(10) = 2 * TG2 * T2(9) - T2(8) T1(11) = 2 * TG1 * T1(10) - T1(9) T2(11) = 2 * TG2 * T2(10) - T2(9) T1(12) = 2 * TG1 * T1(11) - T1(10) T2(12) = 2 * TG2 * T2(11) - T2(10) T1(13) = 2 * TG1 * T1(12) - T1(11) T2(13) = 2 * TG2 * T2(12) - T2(11) DO K=1,NDERIV+1 RES(K) = 0.0_dp RES(K)=RES(K)+DOT_PRODUCT(T1(0:13),C0(1:14,K))*T2(0) RES(K)=RES(K)+DOT_PRODUCT(T1(0:12),C0(15:27,K))*T2(1) RES(K)=RES(K)+DOT_PRODUCT(T1(0:11),C0(28:39,K))*T2(2) RES(K)=RES(K)+DOT_PRODUCT(T1(0:10),C0(40:50,K))*T2(3) RES(K)=RES(K)+DOT_PRODUCT(T1(0:9),C0(51:60,K))*T2(4) RES(K)=RES(K)+DOT_PRODUCT(T1(0:8),C0(61:69,K))*T2(5) RES(K)=RES(K)+DOT_PRODUCT(T1(0:7),C0(70:77,K))*T2(6) RES(K)=RES(K)+DOT_PRODUCT(T1(0:6),C0(78:84,K))*T2(7) RES(K)=RES(K)+DOT_PRODUCT(T1(0:5),C0(85:90,K))*T2(8) RES(K)=RES(K)+DOT_PRODUCT(T1(0:4),C0(91:95,K))*T2(9) RES(K)=RES(K)+DOT_PRODUCT(T1(0:3),C0(96:99,K))*T2(10) RES(K)=RES(K)+DOT_PRODUCT(T1(0:2),C0(100:102,K))*T2(11) RES(K)=RES(K)+DOT_PRODUCT(T1(0:1),C0(103:104,K))*T2(12) RES(K)=RES(K)+DOT_PRODUCT(T1(0:0),C0(105:105,K))*T2(13) ENDDO END SUBROUTINE PD2VAL END MODULE T_C_G0 USE T_C_G0 IMPLICIT NONE INTEGER, PARAMETER :: Nder=9 REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: ref REAL(KIND=dp) :: t1,t2,R,T REAL(KIND=dp) :: RES(0:Nder),MAXERR(0:Nder) INTEGER :: I,J,Nline LOGICAL :: use_gamma CALL INIT(Nder,12) Nline=5000000 CALL CPU_TIME(T1) DO I=1,Nline R=7 T=7 CALL T_C_G0_n(RES,use_gamma,R,T,Nder) ENDDO CALL CPU_TIME(T2) write(6,'(A,F8.3)') "Time for evaluation [s]: ",t2-t1 END