1!==================================================================================================================== 2! * Copyright (c) 2008, 2009 Joost VandeVondele and Manuel Guidon 3! * All rights reserved. 4! * 5! * Redistribution and use in source and binary forms, with or without 6! * modification, are permitted provided that the following conditions are met: 7! * * Redistributions of source code must retain the above copyright 8! * notice, this list of conditions and the following disclaimer. 9! * * Redistributions in binary form must reproduce the above copyright 10! * notice, this list of conditions and the following disclaimer in the 11! * documentation and/or other materials provided with the distribution. 12! * 13! * THIS SOFTWARE IS PROVIDED BY Joost VandeVondele and Manuel Guidon ''AS IS'' AND ANY 14! * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 15! * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 16! * DISCLAIMED. IN NO EVENT SHALL Joost VandeVondele or Manuel Guidon BE LIABLE FOR ANY 17! * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 18! * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 19! * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND 20! * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 21! * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 22! * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 23! 24! Joost VandeVondele and Manuel Guidon, Nov 2008. 25! 26! recursively bisect the rectangular domain of a bivariate function 27! till the lagrange interpolation of specified degree through 28! the Padua points for it converges to a preset threshold 29! generate a Fortran module for the efficient evaluation 30! This program requires the auxiliary program drvrec_mpfr.x 31! 32!==================================================================================================================== 33 MODULE RECURSE 34 USE function_types, ONLY: fun_trunc_coulomb_farfield,fun_trunc_coulomb_nearfield,fun_yukawa 35 36 IMPLICIT NONE 37 INTEGER, SAVE :: find_count=0 38 INTEGER, PARAMETER :: dp=KIND(0.0D0) 39 REAL(KIND=dp) ERR 40 integer :: functionid,digits 41 CONTAINS 42 43!==================================================================================================================== 44! 45! The GO_XXXX are function specific wrapper routines to the bisection routines. 46! 47! they must open a unit 13 (the coefficients of the interpolation) 48! they must open a unit 17 (some info on the bisected domains ) 49! 50! they must match the function definitions / domains in functions.f90 51! 52!==================================================================================================================== 53!==================================================================================================================== 54! 55! GO_truncated 56! 57! This version computes the basic integrals for the 58! truncated coulomb operator 59! 60! G_0(R,T)= ((2*erf(sqrt(t))+erf(R-sqrt(t))-erf(R+sqrt(t)))/sqrt(t)) 61! 62! and up to 21 derivatives with respect to T 63! 64! (-1)**n d^n/dT^n G_0(R,T) 65! 66! 67! This version basically splits the full interval (R>=0 T>=0) in 4 domains 68! The function is only computed for values of R,T which fulfil 69! 70! R**2 - 11.0_dp*R + 0.0_dp < T < R**2 + 11.0_dp*R + 50.0_dp 71! 72! 1) for T larger than the upper bound, 0 is returned, which is accurate at least up to 1.0E-16 73! 2) for T smaller than the lower bound, the caller is instructed to use the gamma function instead 74! 3) for R>11 (farfield), one rectangular domain is bisected 75! 4) for R<11 (nearfield), another rectangle domain is bisected 76! 77!==================================================================================================================== 78 SUBROUTINE GO_truncated() 79 IMPLICIT NONE 80 REAL(KIND=dp) A1,B1,A2,B2 81 INTEGER :: DEG,LEVEL 82 INTEGER :: nderiv 83 84 85 ! these are the selected settings for CP2K 86 ! some initial guess for the needed number of digits 87 digits=128 88 ! number of simultaneous function values to be evaluated for a given 2D point 89 nderiv=21 90 ! degree of the interpolating polynomial 91 DEG=13 92 ! target error 93 ERR=1.0E-9 94 95 OPEN(UNIT=13,FILE="T_C_G.dat") 96 LEVEL=0 97 98 CALL write_copy_right() 99 100 WRITE(6,'(A)') "!" 101 WRITE(6,'(A)') "! This module computes the basic integrals for the" 102 WRITE(6,'(A)') "! truncated coulomb operator" 103 WRITE(6,'(A)') "!" 104 WRITE(6,'(A)') "! res(1) =G_0(R,T)= ((2*erf(sqrt(t))+erf(R-sqrt(t))-erf(R+sqrt(t)))/sqrt(t))" 105 WRITE(6,'(A)') "!" 106 WRITE(6,'(A)') "! and up to 21 derivatives with respect to T" 107 WRITE(6,'(A)') "!" 108 WRITE(6,'(A)') "! res(n+1)=(-1)**n d^n/dT^n G_0(R,T)" 109 WRITE(6,'(A)') "!" 110 WRITE(6,'(A)') "!" 111 WRITE(6,'(A)') "! The function is only computed for values of R,T which fulfil" 112 WRITE(6,'(A)') "!" 113 WRITE(6,'(A)') "! R**2 - 11.0_dp*R + 0.0_dp < T < R**2 + 11.0_dp*R + 50.0_dp" 114 WRITE(6,'(A)') "! where R>=0 T>=0" 115 WRITE(6,'(A)') "!" 116 WRITE(6,'(A)') "! for T larger than the upper bound, 0 is returned" 117 WRITE(6,'(A)') "! (which is accurate at least up to 1.0E-16)" 118 WRITE(6,'(A)') "! while for T smaller than the lower bound, the caller is instructed" 119 WRITE(6,'(A)') "! to use the conventional gamma function instead" 120 WRITE(6,'(A)') "! (i.e. the limit of above expression for R to Infinity)" 121 WRITE(6,'(A)') "! instead of a module one could use: INTEGER, PARAMETER :: dp=KIND(0.0D0)" 122 123 CALL write_module_head("T_C_G0",DEG,nderiv,err) 124 125 WRITE(6,'(A,I0,A)') "SUBROUTINE T_C_G0_n(RES,use_gamma,R,T,NDERIV)" 126 WRITE(6,'(A)') " IMPLICIT NONE" 127 WRITE(6,'(A)') " REAL(KIND=dp), INTENT(OUT) :: RES(*)" 128 WRITE(6,'(A)') " LOGICAL, INTENT(OUT) :: use_gamma" 129 WRITE(6,'(A)') " REAL(KIND=dp),INTENT(IN) :: R,T" 130 WRITE(6,'(A)') " INTEGER, INTENT(IN) :: NDERIV" 131 WRITE(6,'(A)') " REAL(KIND=dp) :: upper,lower,X1,X2,TG1,TG2" 132 133 WRITE(6,'(A)') " use_gamma=.FALSE." 134 WRITE(6,'(A)') " upper=R**2 + 11.0_dp*R + 50.0_dp" 135 WRITE(6,'(A)') " lower=R**2 - 11.0_dp*R + 0.0_dp" 136 WRITE(6,'(A)') " IF (T>upper) THEN" 137 WRITE(6,'(A)') " RES(1:NDERIV+1)=0.0_dp" 138 WRITE(6,'(A)') " RETURN" 139 WRITE(6,'(A)') " ENDIF" 140 141 WRITE(6,'(A)') " IF (R<=11.0_dp) THEN" 142 WRITE(6,'(A)') " X2=R/11.0_dp" 143 WRITE(6,'(A)') " upper=R**2 + 11.0_dp*R + 50.0_dp" 144 WRITE(6,'(A)') " lower=0.0_dp" 145 WRITE(6,'(A)') " X1=(T-lower)/(upper-lower)" 146 147 functionid=fun_trunc_coulomb_nearfield 148 open(UNIT=17,FILE="nearfield_domains.dat") 149 A1=0.0D0 150 B1=1.0D0 151 A2=0.0D0 152 B2=1.0D0 153 CALL FIND_DOMAINS(A1,B1,A2,B2,DEG,NDERIV,ERR,LEVEL) 154 close(17) 155 156 WRITE(6,'(A)') " ELSE" 157 158 WRITE(6,'(A)') " IF (T<lower) THEN" 159 WRITE(6,'(A)') " use_gamma=.TRUE." 160 WRITE(6,'(A)') " RETURN" 161 WRITE(6,'(A)') " ENDIF" 162 WRITE(6,'(A)') " X2=11.0_dp/R" 163 WRITE(6,'(A)') " X1=(T-lower)/(upper-lower)" 164 165 functionid=fun_trunc_coulomb_farfield 166 open(UNIT=17,FILE="farfield_domains.dat") 167 A1=0.0D0 168 B1=1.0D0 169 A2=0.0D0 170 B2=1.0D0 171 CALL FIND_DOMAINS(A1,B1,A2,B2,DEG,NDERIV,ERR,LEVEL) 172 close(17) 173 174 WRITE(6,'(A)') " ENDIF" 175 176 WRITE(6,'(A)') "END SUBROUTINE T_C_G0_n" 177 178 CLOSE(UNIT=13) 179 180 CALL WRITE_MODULE_TAIL("T_C_G0",deg) 181 182 END SUBROUTINE GO_truncated 183!==================================================================================================================== 184! 185! GO_yukawa 186! 187! This version computes the basic integrals for the 188! yukawa coulomb operator (G_0 from S. Ten-no, JCP 126, 014108 (2007), Eq. 49) 189! 190! G_0(R,T)= exp(-T)/4 * sqrt(Pi/T) * [ exp(kappa**2)*erfc(kappa)-exp(lambda**2)*erfc(lambda)] 191! 192! kappa=-sqrt(T)+sqrt(R) 193! lambda=sqrt(T)+sqrt(R) 194! 195! and up to 21 derivatives with respect to T 196! 197! (-1)**n d^n/dT^n G_0(R,T) 198! 199! This version uses a single domain R,T>=0 200! 201!==================================================================================================================== 202 SUBROUTINE GO_yukawa() 203 IMPLICIT NONE 204 REAL(KIND=dp) A1,B1,A2,B2,ERR 205 INTEGER :: DEG,LEVEL 206 INTEGER :: nderiv 207 208 ! these are the selected settings for CP2K 209 digits=128 210 ! number of simultaneous function values to be evaluated for a given 2D point 211 nderiv=21 212 ! degree of the interpolating polynomial 213 DEG=13 214 ! target error 215 ERR=1.0E-9 216 217 IF (DEG<3) STOP "Not supported" 218 OPEN(UNIT=13,FILE="yukawa.dat") 219 LEVEL=0 220 221 CALL WRITE_COPY_RIGHT() 222 223 CALL WRITE_MODULE_HEAD("yukawa",DEG,nderiv,err) 224 225 WRITE(6,'(A,I0,A)') "SUBROUTINE yukawa_n(RES,R,T,NDERIV)" 226 WRITE(6,'(A)') " IMPLICIT NONE" 227 WRITE(6,'(A)') " REAL(KIND=dp), INTENT(OUT) :: RES(*)" 228 WRITE(6,'(A)') " REAL(KIND=dp),INTENT(IN) :: R,T" 229 WRITE(6,'(A)') " INTEGER, INTENT(IN) :: NDERIV" 230 WRITE(6,'(A)') " REAL(KIND=dp) :: upper,lower,X1,X2,TG1,TG2" 231 232 WRITE(6,'(A)') " X2=1/(1+sqrt(R))" 233 WRITE(6,'(A)') " X1=1/(1+sqrt(T))" 234 235 functionid=fun_yukawa 236 open(UNIT=17,FILE="farfield_domains.dat") 237 A1=0.00D0 238 B1=1.00D0 239 A2=0.00D0 240 B2=1.00D0 241 242 CALL FIND_DOMAINS(A1,B1,A2,B2,DEG,NDERIV,ERR,LEVEL) 243 244 WRITE(6,'(A,I0,A)') "END SUBROUTINE yukawa_n" 245 246 CLOSE(UNIT=13) 247 248 CALL WRITE_MODULE_TAIL("yukawa",deg) 249 250 END SUBROUTINE GO_yukawa 251!==================================================================================================================== 252! 253! write the top of the module 254! 255!==================================================================================================================== 256 SUBROUTINE WRITE_MODULE_HEAD(mod_name,DEG,nderiv,err) 257 CHARACTER(LEN=*) :: mod_name 258 INTEGER :: nderiv,deg 259 REAL(KIND=dp) :: err 260 261 WRITE(6,'(A,A)') "MODULE ",mod_name 262 WRITE(6,'(A)') " USE kinds, ONLY: dp" 263 WRITE(6,'(A)') " IMPLICIT NONE" 264 WRITE(6,'(A)') " REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE, SAVE :: C0" 265 WRITE(6,'(A,I0)') " INTEGER, PARAMETER :: degree=",DEG 266 WRITE(6,'(A,E18.6)')" REAL(KIND=dp), PARAMETER :: target_error=",err 267 WRITE(6,'(A,I0)') " INTEGER, PARAMETER :: nderiv_max=",nderiv 268 WRITE(6,'(A)') " INTEGER :: nderiv_init=-1" 269 WRITE(6,'(A)') " INTEGER, SAVE :: patches=-1" 270 WRITE(6,'(A)') "CONTAINS" 271 272 END SUBROUTINE WRITE_MODULE_HEAD 273!==================================================================================================================== 274! 275! write the tail of the module 276! 277!==================================================================================================================== 278 SUBROUTINE WRITE_MODULE_TAIL(mod_name,DEG) 279 CHARACTER(LEN=*) :: mod_name 280 INTEGER :: deg 281 282 CALL WRITE_INIT() 283 CALL WRITE_PD2VAL(DEG) 284 WRITE(6,'(A,A)') "END MODULE ", mod_name 285 286 END SUBROUTINE WRITE_MODULE_TAIL 287!==================================================================================================================== 288! 289! writes a BSD style copy right header 290! 291!==================================================================================================================== 292 SUBROUTINE WRITE_COPY_RIGHT() 293 294 WRITE(6,'(A)') "!" 295 WRITE(6,'(A)') "! * Copyright (c) 2008, Joost VandeVondele and Manuel Guidon" 296 WRITE(6,'(A)') "! * All rights reserved." 297 WRITE(6,'(A)') "! *" 298 WRITE(6,'(A)') "! * Redistribution and use in source and binary forms, with or without" 299 WRITE(6,'(A)') "! * modification, are permitted provided that the following conditions are met:" 300 WRITE(6,'(A)') "! * * Redistributions of source code must retain the above copyright" 301 WRITE(6,'(A)') "! * notice, this list of conditions and the following disclaimer." 302 WRITE(6,'(A)') "! * * Redistributions in binary form must reproduce the above copyright" 303 WRITE(6,'(A)') "! * notice, this list of conditions and the following disclaimer in the" 304 WRITE(6,'(A)') "! * documentation and/or other materials provided with the distribution." 305 WRITE(6,'(A)') "! *" 306 WRITE(6,'(A)') "! * THIS SOFTWARE IS PROVIDED BY Joost VandeVondele and Manuel Guidon AS IS AND ANY" 307 WRITE(6,'(A)') "! * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED" 308 WRITE(6,'(A)') "! * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE" 309 WRITE(6,'(A)') "! * DISCLAIMED. IN NO EVENT SHALL Joost VandeVondele or Manuel Guidon BE LIABLE FOR ANY" 310 WRITE(6,'(A)') "! * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES" 311 WRITE(6,'(A)') "! * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;" 312 WRITE(6,'(A)') "! * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND" 313 WRITE(6,'(A)') "! * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT" 314 WRITE(6,'(A)') "! * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS" 315 WRITE(6,'(A)') "! * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE." 316 WRITE(6,'(A)') "!" 317 WRITE(6,'(A)') "!" 318 WRITE(6,'(A)') "! Joost VandeVondele and Manuel Guidon, Nov 2008." 319 320 END SUBROUTINE WRITE_COPY_RIGHT 321!==================================================================================================================== 322! 323! write the table initialization routines 324! 325!==================================================================================================================== 326 SUBROUTINE WRITE_INIT() 327 328 WRITE(6,'(A)') "! iunit contains the data file to initialize the table" 329 WRITE(6,'(A)') "! Nder is the number of derivatives that will actually be used" 330 WRITE(6,'(A)') "SUBROUTINE INIT(Nder,iunit)" 331 WRITE(6,'(A)') " IMPLICIT NONE" 332 WRITE(6,'(A)') " INTEGER, INTENT(IN) :: Nder,iunit" 333 WRITE(6,'(A)') " INTEGER :: I" 334 WRITE(6,'(A)') " REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: chunk" 335 WRITE(6,'(A,I0)') " patches=",find_count 336 WRITE(6,'(A)') ' IF (Nder>nderiv_max) STOP "Reading data for initialization of C0 failed"' 337 WRITE(6,'(A)') " nderiv_init=Nder" 338 WRITE(6,'(A)') " IF(ALLOCATED(C0)) DEALLOCATE(C0)" 339 WRITE(6,'(A)') " ! round up to a multiple of 32 to give some generous alignment for each C0" 340 WRITE(6,'(A)') " ALLOCATE(C0(32*((31+(Nder+1)*(degree+1)*(degree+2)/2)/32),patches))" 341 WRITE(6,'(A)') " ALLOCATE(chunk((nderiv_max+1)*(degree+1)*(degree+2)/2))" 342 WRITE(6,'(A)') " DO I=1,patches" 343 WRITE(6,'(A)') " READ(iunit,*) chunk" 344 WRITE(6,'(A)') " C0(1:(Nder+1)*(degree+1)*(degree+2)/2,I)=chunk(1:(Nder+1)*(degree+1)*(degree+2)/2)" 345 WRITE(6,'(A)') " ENDDO" 346 WRITE(6,'(A)') " DEALLOCATE(chunk)" 347 WRITE(6,'(A)') "END SUBROUTINE INIT" 348 WRITE(6,'(A)') "SUBROUTINE FREE()" 349 WRITE(6,'(A)') " IF(ALLOCATED(C0)) DEALLOCATE(C0)" 350 WRITE(6,'(A)') " nderiv_init=-1" 351 WRITE(6,'(A)') "END SUBROUTINE FREE" 352 END SUBROUTINE WRITE_INIT 353!==================================================================================================================== 354! 355! write the code to evaluate the interpolation 356! 357!==================================================================================================================== 358 SUBROUTINE WRITE_PD2VAL(DEG) 359 360 INTEGER :: DEG 361 362 INTEGER :: I,J 363 364 WRITE(6,'(A)') "SUBROUTINE PD2VAL(RES,NDERIV,TG1,TG2,C0) " 365 WRITE(6,'(A)') " IMPLICIT NONE" 366 WRITE(6,'(A)') " REAL(KIND=dp), INTENT(OUT) :: res(*)" 367 WRITE(6,'(A)') " INTEGER, INTENT(IN) :: NDERIV" 368 WRITE(6,'(A)') " REAL(KIND=dp),INTENT(IN) :: TG1,TG2" 369 WRITE(6,'(A,I0,A)') " REAL(KIND=dp) :: T1(0:",DEG,")" 370 WRITE(6,'(A,I0,A)') " REAL(KIND=dp) :: T2(0:",DEG,")" 371 WRITE(6,'(A,I0,A)') " REAL(KIND=dp), INTENT(IN) :: C0(",(DEG+1)*(DEG+2)/2,",*)" 372 WRITE(6,'(A)') " INTEGER :: I,J,K" 373 WRITE(6,'(A)') " REAL(KIND=dp), PARAMETER :: SQRT2=1.4142135623730950488016887242096980785696718753_dp" 374 WRITE(6,'(A)') " T1(0)=1.0_dp" 375 WRITE(6,'(A)') " T2(0)=1.0_dp" 376 WRITE(6,'(A)') " T1(1)=SQRT2 * TG1" 377 WRITE(6,'(A)') " T2(1)=SQRT2 * TG2" 378 WRITE(6,'(A)') " T1(2)= 2 * TG1 * T1(1) - SQRT2 " 379 WRITE(6,'(A)') " T2(2)= 2 * TG2 * T2(1) - SQRT2 " 380 DO I=3,DEG 381 WRITE(6,'(A,I0,A,I0,A,I0,A)') " T1(",I,") = 2 * TG1 * T1(",I-1,") - T1(",I-2,")" 382 WRITE(6,'(A,I0,A,I0,A,I0,A)') " T2(",I,") = 2 * TG2 * T2(",I-1,") - T2(",I-2,")" 383 ENDDO 384 385 WRITE(6,'(A)') " DO K=1,NDERIV+1" 386 WRITE(6,'(A)') " RES(K) = 0.0_dp" 387 J=1 388 DO I=DEG,0,-1 389 WRITE(6,'(A,I0,A,I0,A,I0,A,I0,A)') " RES(K)=RES(K)+DOT_PRODUCT(T1(0:",I,"),C0(",& 390 J,":",J+I,",K))*T2(",DEG-I,")" 391 J=J+I+1 392 ENDDO 393 WRITE(6,'(A)') " ENDDO" 394 395 WRITE(6,'(A)') "END SUBROUTINE PD2VAL" 396 397 END SUBROUTINE WRITE_PD2VAL 398!==================================================================================================================== 399! 400! The basic piece of the code, does the bisection till convergence 401! 402!==================================================================================================================== 403 RECURSIVE SUBROUTINE FIND_DOMAINS(A1,B1,A2,B2,DEG,NDERIV,ERR,LEVEL) 404 405 IMPLICIT NONE 406 REAL(KIND=dp) A1,B1,A2,B2,ERR,R1,R2 407 INTEGER :: DEG,LEVEL,I,DIRECTION,NDERIV,K,I1,I2 408 409 REAL(KIND=dp) ESTIMATED_ERR,MIDDLE,ESTIMATED_ERR1,ESTIMATED_ERR2,ESTIMATED_ERR3,ESTIMATED_ERR4 410 REAL(KIND=dp) :: C0(0:DEG,0:DEG,0:NDERIV) 411 412 IF (DEG<3) STOP "Not supported" 413 IF (A1==B1 .AND. A2==B2) STOP "BUG 1" 414 IF (LEVEL>20) STOP "Too many bisections: is the function sufficiently differentiable in the full domain?" 415 416 ! compute the interpolating polynomial and its estimated error 417 418 CALL drvrec_wrap(A1,B1,A2,B2,DEG,NDERIV,C0,ESTIMATED_ERR) 419 420 IF (ESTIMATED_ERR<ERR) THEN 421 422 ! we have converged ... turn this into a call to pd2val, and append to the C0 table 423 find_count=find_count+1 424 WRITE(6,'(A,E25.18,A,E25.18,A)') REPEAT(" ",LEVEL+3)//"TG1= ( 2 * X1 - ",(B1+A1),"_dp)*",1.0D0/(B1-A1),"_dp" 425 WRITE(6,'(A,E25.18,A,E25.18,A)') REPEAT(" ",LEVEL+3)//"TG2= ( 2 * X2 - ",(B2+A2),"_dp)*",1.0D0/(B2-A2),"_dp" 426 427 WRITE(6,'(A,I0,A)') REPEAT(" ",LEVEL+3)//"CALL PD2VAL(RES,NDERIV,TG1,TG2,C0(1,", find_count, "))" 428 WRITE(17,'(2F20.16,I10,E16.8)') A1,A2,find_count,ESTIMATED_ERR 429 WRITE(17,'(2F20.16,I10,E16.8)') A1,B2 ,find_count,ESTIMATED_ERR 430 WRITE(17,*) 431 WRITE(17,'(2F20.16,I10,E16.8)') A1,A2 ,find_count,ESTIMATED_ERR 432 WRITE(17,'(2F20.16,I10,E16.8)') B1,A2 ,find_count,ESTIMATED_ERR 433 WRITE(17,*) 434 WRITE(17,'(2F20.16,I10,E16.8)') B1,B2,find_count,ESTIMATED_ERR 435 WRITE(17,'(2F20.16,I10,E16.8)') A1,B2,find_count,ESTIMATED_ERR 436 WRITE(17,*) 437 WRITE(17,'(2F20.16,I10,E16.8)') B1,B2,find_count,ESTIMATED_ERR 438 WRITE(17,'(2F20.16,I10,E16.8)') B1,A2,find_count,ESTIMATED_ERR 439 WRITE(17,*) 440 DO K=0,NDERIV 441 DO I=DEG,0,-1 442 write(13,*) C0(0:I,DEG-I,K) 443 ENDDO 444 ENDDO 445 446 ELSE 447 448 ! decide in which dimension to bisection. 449 ! the easiest and fastest way 450 ! DIRECTION=MOD(LEVEL,2)+1 451 ! The way leading to more efficient code 452 ! using a section that minimizes the maximum error is slightly more 453 ! efficient than alternating bisection. 454 MIDDLE=(A1+B1)/2 455 CALL drvrec_wrap(A1,MIDDLE,A2,B2,DEG,NDERIV,C0,ESTIMATED_ERR1) 456 CALL drvrec_wrap(MIDDLE,B1,A2,B2,DEG,NDERIV,C0,ESTIMATED_ERR2) 457 MIDDLE=(A2+B2)/2 458 CALL drvrec_wrap(A1,B1,A2,MIDDLE,DEG,NDERIV,C0,ESTIMATED_ERR3) 459 CALL drvrec_wrap(A1,B1,MIDDLE,B2,DEG,NDERIV,C0,ESTIMATED_ERR4) 460 IF (MAX(ESTIMATED_ERR1,ESTIMATED_ERR2)<MAX(ESTIMATED_ERR3,ESTIMATED_ERR4)) THEN 461 DIRECTION=1 462 ELSE 463 DIRECTION=2 464 ENDIF 465 466 IF (DIRECTION==1) THEN 467 MIDDLE=(A1+B1)/2 468 WRITE(6,'(A,E25.18,A)') REPEAT(" ",LEVEL)//" IF (X1<=",MIDDLE,"_dp) THEN " 469 CALL FIND_DOMAINS(A1,MIDDLE,A2,B2,DEG,NDERIV,ERR,LEVEL+1) 470 WRITE(6,'(A)') REPEAT(" ",LEVEL)//" ELSE" 471 CALL FIND_DOMAINS(MIDDLE,B1,A2,B2,DEG,NDERIV,ERR,LEVEL+1) 472 WRITE(6,'(A)') REPEAT(" ",LEVEL)//" ENDIF" 473 ELSE 474 MIDDLE=(A2+B2)/2 475 WRITE(6,'(A,E25.18,A)') REPEAT(" ",LEVEL)//" IF (X2<=",MIDDLE,"_dp) THEN " 476 CALL FIND_DOMAINS(A1,B1,A2,MIDDLE,DEG,NDERIV,ERR,LEVEL+1) 477 WRITE(6,'(A)') REPEAT(" ",LEVEL)//" ELSE" 478 CALL FIND_DOMAINS(A1,B1,MIDDLE,B2,DEG,NDERIV,ERR,LEVEL+1) 479 WRITE(6,'(A)') REPEAT(" ",LEVEL)//" ENDIF" 480 ENDIF 481 ENDIF 482 483 END SUBROUTINE FIND_DOMAINS 484!==================================================================================================================== 485! 486! evaluates C0 coefs increasing the number of digits, till the resulting C0 is converged. 487! This is a simple wrapper to the program drvrec_mpfr.x 488! 489!==================================================================================================================== 490 SUBROUTINE drvrec_wrap(A1,B1,A2,B2,DEG,NDERIV,C0,ESTIMATED_ERR) 491 492 REAL(KIND=dp) :: a1,b1,a2,b2,estimated_err,estimated_err2,estimated_err3,esterr 493 INTEGER :: deg,nderiv 494 REAL(KIND=dp) :: C0(0:DEG,0:DEG,0:NDERIV), C02(0:DEG,0:DEG,0:NDERIV) 495 REAL(KIND=dp) :: R,T,vals(0:nderiv) 496 integer :: digits_local 497 498 digits_local=digits 499 C02 = 0.0_dp 500 501 DO 502 OPEN(UNIT=51,FILE="drvrec_mpfr.in") 503 WRITE(51,*) digits_local 504 WRITE(51,*) deg 505 WRITE(51,*) nderiv 506 WRITE(51,*) A1,B1 507 WRITE(51,*) A2,B2 508 WRITE(51,*) functionid 509 CLOSE(51) 510 CALL SYSTEM("./drvrec_mpfr.x") 511 OPEN(UNIT=51,FILE="drvrec_mpfr.out") 512 READ(51,*) ESTIMATED_ERR 513 READ(51,*) C0 514 CLOSE(51) 515 516 ! do not allow for differences larger than ERR times a given factor 517 IF(MAXVAL(ABS(C0-C02))>ERR*1.0E-5) THEN 518 C02=C0 519 digits_local=digits_local*2 520 ELSE 521 OPEN(UNIT=97,FILE="function_vals.dat") 522 DO 523 READ(97,*,END=999) R,T,vals 524 WRITE(98,*) R,T,vals 525 ENDDO 526999 CONTINUE 527 CLOSE(97) 528 EXIT 529 ENDIF 530 ENDDO 531 532 END SUBROUTINE 533END MODULE Recurse 534 535!==================================================================================================================== 536! 537! The main program ... small and easy 538! 539! output files.... 540! std out : a program that can evaluate the function in a given domain 541! history.dat : a list of function values, evaluated during the construction of the approximation, 542! that can be used as a reference data. 543! 544!==================================================================================================================== 545PROGRAM Fun2D 546 547 USE Recurse 548 549 OPEN(UNIT=98,FILE="history.dat") 550 551#if defined(__T_C_G0) 552 CALL GO_truncated() 553#endif 554#if defined(__YUKAWA) 555 CALL GO_yukawa() 556#endif 557 558 CLOSE(98) 559 560END PROGRAM Fun2D 561