comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-27 21:55:16 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Computing 2D functions: source generator for doing Lagrange interpolation through the Padua points in automatically bisected domains.

Browse source 
svn-origin-rev: 8033
This commit is contained in:
Joost VandeVondele 2008-12-09 11:11:43 +00:00
parent 2445f7e935
commit 93e8156b48
10 changed files with 5702 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

34
tools/Fun2D/Makefile Normal file
View file

@ -0,0 +1,34 @@
#
# requires a sufficiently recent gfortran and mpfr ....
#
FC = gfortran -ffree-line-length-512 -g
FFLAGS = -O3 -march=native -funroll-loops
LIBSBLAS = /ext/software/64/opt4/gfortran-suite/lib/libgoto.a -lpthread
LIBSMPFR = -lmpfr -lgmp
mpfr_cutoff_gamma.o: mpfr_cutoff_gamma.f90
$(FC) -c mpfr.f90 mpfr_cutoff_gamma.f90
evalf: mpfr_cutoff_gamma.o
$(FC) -o evalf.x mpfr.f90 mpfr_cutoff_gamma.o evalf.f90 $(LIBSMPFR)
drvrec_mpfr: mpfr_cutoff_gamma.o
$(FC) -o drvrec_mpfr.x mpfr.f90 padua2_mpfr.f90 mpfr_cutoff_gamma.o drvrec_mpfr.f90 $(LIBSMPFR)
recurse: drvrec_mpfr evalf
$(FC) -o recurse.x recurse.f90 $(LIBSBLAS)
t_c_g0: recurse
./recurse.x > t_c_g0.f90
tester: t_c_g0
$(FC) $(FFLAGS) -o tester.x kinds.f90 t_c_g0.f90 tester.f90
test: tester
./tester.x
clean:
rm -f *.o *.x *.mod *.MOD
realclean: clean
rm -f *.dat *.in *.out t_c_g0.f90 fort.*

22
tools/Fun2D/README Normal file
View file

@ -0,0 +1,22 @@
A package to generate Fortran code for the evaluation of bivariate functions.
Using recursive bisection and lagrane interpolation at padua points, the code
and necessary tables are being constructed.
The function (f) must be defined in drvrec_mpfr.f90,
and its domain must be (transformed to) rectangular.
The settings for the bisection, accuracy, ... are set in recurse.f90. Most important
are the degree of the interpolation, and the target error.
The actual error is usually smaller than the estimate of the error.
The degree of the interpolation and the target error influence the size of the tables,
and speed of evaluation.
The code is interfaced to the arbitrary precision library MPFR, with a convenient Fortran
wrapper. However, as long as derived type finalization (a Fortran2003 feature) is not available,
memory is not freed during MPFR math, and thus a significant amount of memory is needed. To reduce the
impact, drvrec_mpfr has been transformed into a little program, so that the OS frees the memory in between
calls.
The code for interpolation through the Padua points is based on an older version of Padua2D,
kindly made available by the authors (Marco Caliari and Marco Vianello and Stefano De Marchi)
see top of padua2_mpfr.f90 for further credits.

177
tools/Fun2D/drvrec_mpfr.f90 Normal file
View file

@ -0,0 +1,177 @@
!
! a simple wrapper program to
! obtain the C0 coefficients for an NDERIV+1 valued function in the
! rectangular domain A1<=X1<=B1, A2<=X2<=B2
!
! input from 'drvrec_mpfr.in'
! output to 'drvrec_mpfr.out'
!
SUBROUTINE drvrec(a1,b1,a2,b2,deg,nderiv,c0,esterr_out)
USE padua2_mpfr
IMPLICIT NONE
INTEGER deg,npd,nderiV
TYPE(mpfr_type) pd1((deg + 1) * (deg + 2) / 2)
TYPE(mpfr_type) pd2((deg + 1) * (deg + 2) / 2)
TYPE(mpfr_type) wpd((deg + 1) * (deg + 2) / 2)
TYPE(mpfr_type) fpd((deg + 1) * (deg + 2) / 2,0:nderiv)
TYPE(mpfr_type) t1(0:deg,(deg + 1) * (deg + 2) / 2)
TYPE(mpfr_type) t2(0:deg,(deg + 1) * (deg + 2) / 2)
TYPE(mpfr_type) c0(0:deg,0:deg,0:nderiv)
TYPE(mpfr_type) res(0:nderiv),esterr(0:nderiv)
TYPE(mpfr_type) a1,b1,a2,b2,esterr_out,x1,x2,zero
INTEGER i,ntg1
CALL pdpts(deg,pd1,pd2,wpd,npd)
DO i = 1,npd
x1=((b1 - a1) * pd1(i) + (b1 + a1)) / 2
x2=((b2 - a2) * pd2(i) + (b2 + a2)) / 2
CALL f(res,nderiv,x1,x2)
FPD(I,:) = RES(0:NDERIV)
ENDDO
DO i=0,nderiv
CALL padua2(deg,npd,pd1,pd2,wpd,fpd(:,i),t1,t2,c0(:,:,i),esterr(i))
ENDDO
esterr_out=0
zero=0
DO i=0,nderiv
IF (mpfr_cmp(esterr(i),zero)>0) THEN
IF (mpfr_cmp(esterr(i),esterr_out)>0) THEN
esterr_out=esterr(i)
ENDIF
ELSE
IF (mpfr_cmp(-esterr(i),esterr_out)>0) esterr_out=-esterr(i)
ENDIF
ENDDO
END SUBROUTINE drvrec
! evaluate the function, + nderiv other functions in the point X1,X2
! write additionally to a file 'history.res' for later checking the quality of the interpolation
SUBROUTINE f(res,nderiv,x1,x2)
USE mpfr
USE mpfr_ops
USE mpfr_cutoff_gamma
IMPLICIT NONE
INTEGER :: nderiv
TYPE(mpfr_type) :: x1,x2
TYPE(mpfr_type) :: res(0:nderiv)
TYPE(mpfr_type) :: T,r,upper,lower,zero
TYPE(mpfr_type) :: dummy(0:21)
INTEGER :: i
INTEGER :: functionid,should_output
COMMON/functiontype/functionid,should_output
zero=.CONVERT."0"
! functionid in this case just selects a different
! domain of the same function
! 1 = farfield (R>11)
! 2 = nearfield (R<11)
SELECT CASE(functionid)
CASE (1)
! if R=+Infinity the result is zero
IF (mpfr_cmp(x2,zero)<=0) THEN
res=.CONVERT."0"
return
ENDIF
r=11/x2
upper=r**2 + 11*r + 50
lower=r**2 - 11*r
CASE (2)
R=x2*11
upper=r**2 + 11*r + 50
lower=.CONVERT."0"
CASE DEFAULT
STOP "Function ID not implemented"
END SELECT
t=lower+x1*(upper-lower)
!t is zero, return the limiting expansion
IF (mpfr_cmp(t,zero)<=0) THEN
CALL cutoff_gamma_T0(21,R,dummy)
ELSE
CALL cutoff_gamma(21,t,r,dummy)
END IF
res(0:nderiv)=dummy(0:nderiv)
IF (should_output>0) THEN
WRITE(should_output,*) REAL(R),REAL(T),REAL(res)
ENDIF
END SUBROUTINE f
!
! this program is a simple interface to obtain the C0 coefficients
! with high accuracy (i.e. all digits & correctly rounded),
! for a function programmed in subroutine f
! this version uses double precision (~16 digits)
! input/output but does all calculations
! internally using a user specified number of digits
!
PROGRAM drvrec_padua
USE mpfr
USE mpfr_ops
USE padua2_mpfr, ONLY: pd2val_mp=>pd2val
IMPLICIT NONE
INTEGER :: deg,nderiv,i,j,digits,deg_ref,k,l
INTEGER*2 :: precision
REAL(KIND=dp) :: esterr,A1,B1,A2,B2
REAL(KIND=dp), DIMENSION(:,:,:), ALLOCATABLE :: c0
TYPE(mpfr_type) :: A1_mp,B1_mp,A2_mp,B2_mp,esterr_mp
TYPE(mpfr_type), DIMENSION(:,:,:), ALLOCATABLE :: c0_mp
integer :: functionid,should_output
COMMON/functiontype/functionid,should_output
! keeps a history of evaluated function values.
! useful to investigate the error later on
should_output=31
IF (should_output>0) THEN
OPEN(should_output,FILE="history.dat",ACCESS="APPEND")
ENDIF
! input definition
OPEN(17,FILE="drvrec_mpfr.in")
READ(17,*) digits
READ(17,*) deg
READ(17,*) nderiv
READ(17,*) A1,B1
READ(17,*) A2,B2
READ(17,*) functionid
CLOSE(17)
precision = log(10.0)/log(2.0)*digits
CALL mpfr_set_default_precision(precision)
A1_mp=A1
B1_mp=B1
A2_mp=A2
B2_mp=B2
ALLOCATE(c0_mp(0:deg,0:deg,0:nderiv))
CALL drvrec(a1_mp,b1_mp,a2_mp,b2_mp,deg,nderiv,c0_mp,esterr_mp)
! convert ot double precision and output
ALLOCATE(c0(0:deg,0:deg,0:nderiv))
c0=REAL(c0_mp)
esterr=REAL(esterr_mp)
OPEN(17,FILE="drvrec_mpfr.out")
write(17,*) esterr
write(17,*) c0
CLOSE(17)
IF (should_output>0) THEN
CLOSE(should_output)
ENDIF
END PROGRAM drvrec_padua

32
tools/Fun2D/evalf.f90 Normal file
View file

@ -0,0 +1,32 @@
!
! a simple interface to the cutoff_gamma,
! reads from 'evalf.in' writes to 'evalf.out'
!
PROGRAM main
USE mpfr_cutoff_gamma
IMPLICIT NONE
REAL(dp) :: T,R,dummy(0:21)
INTEGER*2 :: precision
TYPE(mpfr_type):: Tval_mp, Rval_mp, dummy_mp(0:21)
INTEGER :: i
OPEN(UNIT=77,FILE="evalf.in")
READ(77,*) R,T
CLOSE(77)
precision = log(10.0)/log(2.0)*4096
CALL mpfr_set_default_precision(precision)
Tval_mp = T
Rval_mp = R
IF(T==0.0_dp) THEN
CALL cutoff_gamma_T0(21,Rval_mp,dummy_mp)
ELSE
CALL cutoff_gamma(21,Tval_mp,Rval_mp,dummy_mp)
END IF
DO i=0,21
dummy(i)=mp_to_real(dummy_mp(i))
END DO
OPEN(UNIT=78,FILE="evalf.out")
write(78,*) R,T,dummy(:)
CLOSE(78)
END PROGRAM main

3
tools/Fun2D/kinds.f90 Normal file
View file

@ -0,0 +1,3 @@
MODULE kinds
INTEGER, PARAMETER :: dp = KIND(0.0D0)
END MODULE

2053
tools/Fun2D/mpfr.f90 Normal file
View file

File diff suppressed because it is too large Load diff

2644
tools/Fun2D/mpfr_cutoff_gamma.f90 Normal file
View file

File diff suppressed because it is too large Load diff

266
tools/Fun2D/padua2_mpfr.f90 Normal file
View file

@ -0,0 +1,266 @@
!
! These routines have been converted to use mpfr variables
! from an 'older' version of Padua2D.
!
! In this way, arbitrary precision can be used to construct the
! coefficients and evaluate the interpolation
!
! Joost VandeVondele, November 2008
!
! ******************************************************************
!
! From Padua2D
!
! Marco Caliari and Marco Vianello
! Department of Pure and Applied Mathematics
! University of Padua (Italy)
! {mcaliari,marcov}@math.unipd.it
!
! Stefano De Marchi
! Department of Computer Science
! University of Verona (Italy)
! stefano.demarchi@univr.it
!
! June 19th, 2007
!
MODULE padua2_mpfr
USE mpfr
USE mpfr_ops
IMPLICIT NONE
CONTAINS
! ******************************************************************
!
! PARAMS type i=input, o=output, a=auxiliary
!
! DEG integer i maximum degree
! PT mpfr i point
! TCHEB mpfr(0:DEG) o Normalized Chebyshev polynomials at
! the point PT
! S1 mpfr i scale factor for TCHEB
!
! **********************************************************************
SUBROUTINE cheb(deg,pt,tcheb,s1)
INTEGER :: deg
TYPE(mpfr_type) :: pt,s1,tcheb(0:deg)
TYPE(mpfr_type) :: sqrt2
INTEGER :: j
sqrt2="2"
sqrt2=sqrt(sqrt2)
IF (deg .GE. 0) THEN
tcheb(0) = s1
END IF
IF (deg .GE. 1) THEN
tcheb(1) = sqrt2 * pt * s1
END IF
IF (deg .GE. 2) THEN
tcheb(2) = 2 * pt * tcheb(1) - sqrt2 * tcheb(0)
END IF
DO J = 3,DEG
tcheb(j) = 2 * pt * tcheb(j - 1) - tcheb(j - 2)
ENDDO
END SUBROUTINE cheb
! ******************************************************************
!
! PARAMS type i=input, o=output, a=auxiliary
!
! DEG integer i degree of approximation
! PD1 mpfr(*) o array of the first coordinates of
! Padua points
! PD2 mpfr(*) o array of the second coordinates of
! Padua points
! WPD mpfr(*) o array of the weights
! NPD integer o number of Padua points =
! (DEG + 1) * (DEG + 2) / 2
!
! ******************************************************************
SUBROUTINE pdpts(deg,pd1,pd2,wpd,npd)
INTEGER :: deg,npd
TYPE(mpfr_type) :: pd1(*),pd2(*),wpd(*)
INTEGER j,k,deg1,itemp0,kold,retval
TYPE(mpfr_type) :: pi,rtemp0
pi=1.0D0
pi=atan(pi)*4
IF (deg .EQ. 0) THEN
pd1(1) = -1
pd2(1) = -1
WPD(1) = 2
npd = 1
RETURN
END IF
npd = 0
deg1 = deg + 1
itemp0 = deg * deg1
rtemp0 = pi / itemp0
k = 0
npd = npd + 1
!
! The points are generated starting from bottom-left
! and following the generating curve
!
pd1(npd) = -1
pd2(npd) = -1
wpd(npd) = .CONVERT."1" / (2*itemp0)
DO k = 1,deg - 1
npd = npd + 1
pd1(npd) = -COS(deg1 * k * rtemp0)
pd2(npd) = -COS(deg * k * rtemp0)
wpd(npd) = .CONVERT."2" / itemp0
ENDDO
k = deg
npd = npd + 1
pd1(npd) = -COS(deg1 * k * rtemp0)
pd2(npd) = -COS(deg * k * rtemp0)
wpd(npd) = .CONVERT."1" / itemp0
kold = deg1
DO j = 0,deg - 2
k = kold
npd = npd + 1
pd1(npd) = -COS(deg1 * k * rtemp0)
pd2(npd) = -COS(deg * k * rtemp0)
wpd(npd) = .CONVERT."1" / itemp0
DO k = kold + 1,kold + deg - j - 2
npd = npd + 1
pd1(npd) = -COS(deg1 * k * rtemp0)
pd2(npd) = -COS(deg * k * rtemp0)
wpd(npd) = .CONVERT."2" / itemp0
ENDDO
k = kold + deg - j - 1
npd = npd + 1
pd1(npd) = -COS(deg1 * k * rtemp0)
pd2(npd) = -COS(deg * k * rtemp0)
wpd(npd) = .CONVERT."1" / itemp0
kold = kold + deg1
ENDDO
k = itemp0
npd = npd + 1
pd1(npd) = -(2 * MOD(DEG,2) - 1)
pd2(npd) = -pd1(npd)
wpd(npd) = .CONVERT."1" / (2 * itemp0)
END SUBROUTINE pdpts
! ******************************************************************
!
! PARAMS type i=input, o=output, a=auxiliary
!
! DEG integer i degree of approximation
! DEGMAX integer i maximum degree allowed. DEGMAX + 1 is
! the leading dimension of T1, T2 and
! C0
! NPD integer i number of Padua points
! PD1 mpfr(NPD) i array of the first coordinates of
! Padua points
! PD2 mpfr(NPD) i array of the second coordinates of
! Padua points
! WPD mpfr(NPD) i array of the weights
! FPD mpfr(NPD) i array of the values of the function
! at Padua points
! T1 mpfr(0:DEG,NPD) a Normalized Chebyshev polynomials at
! the first coordinate of Padua points
! T2 mpfr(0:DEG,NPD) a Normalized Chebyshev polynomials at
! the second coordinate of Padua points
! C0 mpfr(0:DEG,0:DEG) o matrix C0 of the coefficients
! ESTERR mpfr o estimated error
!
! ******************************************************************
SUBROUTINE padua2(deg,npd,pd1,pd2,wpd,fpd,t1,t2,c0,esterr)
INTEGER :: deg,npd
TYPE(mpfr_type) :: esterr
TYPE(mpfr_type) :: pd1(npd),pd2(npd),wpd(npd),fpd(npd), &
t1(0:deg,npd),t2(0:deg,npd),c0(0:deg,0:deg),rtmp
INTEGER :: i,j,k,retval
rtmp=1
DO i = 1,npd
CALL cheb(deg,pd1(i),t1(0,i),wpd(i) * fpd(i))
CALL cheb(deg,pd2(i),t2(0,i),rtmp)
ENDDO
! B = T1 * T2'
rtmp=0
DO i=0,deg
DO j=0,deg
c0(i,j)=0
DO k=1,npd
! reduce the impact of the memory leak problem, code explicitly
! in mpfr operations the O(n^4) flops in this section
! c0(i,j)=c0(i,j)+t1(i,k)*t2(j,k)
retval = mpfr_mul(rtmp,t1(i,k),t2(j,k),GMP_RNDN)
retval = mpfr_add(c0(i,j),c0(i,j),rtmp,GMP_RNDN)
ENDDO
ENDDO
ENDDO
! Halve the correct term in C0
c0(deg,0) = c0(deg,0) / 2
esterr = 0
rtmp = 0
DO j = 0,2
DO i = 0,deg - j
IF (mpfr_cmp(c0(i,deg - i - j),rtmp)>0) THEN
esterr = esterr + c0(i,deg - i - j)
ELSE
esterr = esterr - c0(i,deg - i - j)
ENDIF
ENDDO
ENDDO
esterr = 2 * esterr
END SUBROUTINE padua2
! ******************************************************************
!
! PARAMS type i=input, o=output, a=auxiliary
!
! DEG integer i degree of approximation
! DEGMAX integer i maximum degree allowed. DEGMAX + 1 is
! the leading dimension of C0
! C0 double(0:DEG,0:DEG) i matrix C0 of the coefficients
! TG1 double i first coordinate of the
! target point
! TG2 double i second coordinate of the
! target point
! T1 double(0:DEG) a Normalized Chebyshev polynomials at
! the first coordinate of target point
! T2 double(0:DEG) a Normalized Chebyshev polynomials at
! the second coordinate of target point
!
! **********************************************************************
FUNCTION pd2val(deg,c0,tg1,tg2,t1,t2) RESULT(res)
INTEGER :: deg
TYPE(mpfr_type) :: res,tg1,tg2,c0(0:deg,0:deg),t1(0:deg),t2(0:deg)
INTEGER :: i,k
TYPE(mpfr_type) :: dot,one
one = 1
CALL cheb(deg,tg1,t1,one)
CALL cheb(deg,tg2,t2,one)
res = 0
DO i = deg,0,-1
dot = 0
DO k=0,i
dot = dot + t1(k)*c0(k,deg-i)
ENDDO
res = res + dot * t2(deg-i)
ENDDO
END FUNCTION pd2val
END MODULE padua2_mpfr

403
tools/Fun2D/recurse.f90 Normal file
View file

@ -0,0 +1,403 @@
! * Copyright (c) 2008, Joost VandeVondele and Manuel Guidon
! * All rights reserved.
! *
! * Redistribution and use in source and binary forms, with or without
! * modification, are permitted provided that the following conditions are met:
! * * Redistributions of source code must retain the above copyright
! * notice, this list of conditions and the following disclaimer.
! * * Redistributions in binary form must reproduce the above copyright
! * notice, this list of conditions and the following disclaimer in the
! * documentation and/or other materials provided with the distribution.
! *
! * THIS SOFTWARE IS PROVIDED BY Joost VandeVondele and Manuel Guidon ''AS IS'' AND ANY
! * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! * DISCLAIMED. IN NO EVENT SHALL Joost VandeVondele or Manuel Guidon BE LIABLE FOR ANY
! * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
! * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! recursively bisect the rectangular domain of a bivariate function
! till the lagrange interpolation of specified degree through
! the Padua points for it converges to a preset threshold
! generate a Fortran module for the efficient evaluation
!
! This version computes the basic integrals for the
! truncated coulomb operator
!
! G_0(R,T)= ((2*erf(sqrt(t))+erf(R-sqrt(t))-erf(R+sqrt(t)))/sqrt(t))
!
! and up to 21 derivatives with respect to T
!
! (-1)**n d^n/dT^n G_0(R,T)
!
!
! The function is only computed for values of R,T which fulfil
!
! R**2 - 11.0_dp*R + 0.0_dp < T < R**2 + 11.0_dp*R + 50.0_dp
! where R>=0 T>=0
!
! for T larger than the upper bound, 0 is returned
! (which is accurate at least up to 1.0E-16)
! while for T smaller than the lower bound, the caller is instructed
! to use the gamma function instead
!
! Joost VandeVondele and Manuel Guidon, Nov 2008.
!
! This program requires the auxiliary program drvrec_mpfr.x
!
MODULE RECURSE
INTEGER, SAVE :: find_count=0
INTEGER, PARAMETER :: dp=KIND(0.0D0)
CONTAINS
SUBROUTINE GO()
IMPLICIT NONE
REAL(KIND=dp) A1,B1,A2,B2,ERR
INTEGER :: DEG,LEVEL
INTEGER :: nderiv
integer :: functionid,digits
COMMON/functiontype/functionid,digits
! these are the selected settings for CP2K
digits=1024
! number of simultaneous function values to be evaluated for a given 2D point
nderiv=21
! degree of the interpolating polynomial
DEG=13
! target error
ERR=1.0E-9
IF (DEG<3) STOP "Not supported"
OPEN(UNIT=13,FILE="T_C_G.dat")
LEVEL=0
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "! * Copyright (c) 2008, Joost VandeVondele and Manuel Guidon"
WRITE(6,'(A)') "! * All rights reserved."
WRITE(6,'(A)') "! *"
WRITE(6,'(A)') "! * Redistribution and use in source and binary forms, with or without"
WRITE(6,'(A)') "! * modification, are permitted provided that the following conditions are met:"
WRITE(6,'(A)') "! * * Redistributions of source code must retain the above copyright"
WRITE(6,'(A)') "! * notice, this list of conditions and the following disclaimer."
WRITE(6,'(A)') "! * * Redistributions in binary form must reproduce the above copyright"
WRITE(6,'(A)') "! * notice, this list of conditions and the following disclaimer in the"
WRITE(6,'(A)') "! * documentation and/or other materials provided with the distribution."
WRITE(6,'(A)') "! *"
WRITE(6,'(A)') "! * THIS SOFTWARE IS PROVIDED BY Joost VandeVondele and Manuel Guidon AS IS AND ANY"
WRITE(6,'(A)') "! * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED"
WRITE(6,'(A)') "! * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE"
WRITE(6,'(A)') "! * DISCLAIMED. IN NO EVENT SHALL Joost VandeVondele or Manuel Guidon BE LIABLE FOR ANY"
WRITE(6,'(A)') "! * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES"
WRITE(6,'(A)') "! * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;"
WRITE(6,'(A)') "! * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND"
WRITE(6,'(A)') "! * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT"
WRITE(6,'(A)') "! * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS"
WRITE(6,'(A)') "! * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE."
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "! This module computes the basic integrals for the"
WRITE(6,'(A)') "! truncated coulomb operator"
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "! res(1) =G_0(R,T)= ((2*erf(sqrt(t))+erf(R-sqrt(t))-erf(R+sqrt(t)))/sqrt(t))"
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "! and up to 21 derivatives with respect to T"
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "! res(n+1)=(-1)**n d^n/dT^n G_0(R,T)"
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "! The function is only computed for values of R,T which fulfil"
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "! R**2 - 11.0_dp*R + 0.0_dp < T < R**2 + 11.0_dp*R + 50.0_dp"
WRITE(6,'(A)') "! where R>=0 T>=0"
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "! for T larger than the upper bound, 0 is returned"
WRITE(6,'(A)') "! (which is accurate at least up to 1.0E-16)"
WRITE(6,'(A)') "! while for T smaller than the lower bound, the caller is instructed"
WRITE(6,'(A)') "! to use the conventional gamma function instead"
WRITE(6,'(A)') "! (i.e. the limit of above expression for R to Infinity)"
WRITE(6,'(A)') "!"
WRITE(6,'(A)') "! Joost VandeVondele and Manuel Guidon, Nov 2008."
WRITE(6,'(A)') "! instead of a module one could use: INTEGER, PARAMETER :: dp=KIND(0.0D0)"
WRITE(6,'(A,I0,A)') "MODULE T_C_G0"
WRITE(6,'(A)') " USE kinds, ONLY: dp"
WRITE(6,'(A)') " IMPLICIT NONE"
WRITE(6,'(A)') " PUBLIC :: T_C_G0_n"
WRITE(6,'(A)') " REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE, SAVE :: C0"
WRITE(6,'(A,I0)') " INTEGER, PARAMETER :: degree=",DEG
WRITE(6,'(A,E18.6)') " REAL(KIND=dp), PARAMETER :: target_error=",err
WRITE(6,'(A,I0)') " INTEGER, PARAMETER :: nderiv_max=",nderiv
WRITE(6,'(A,I0)') " INTEGER :: nderiv_init=",-1
WRITE(6,'(A)') " INTEGER, SAVE :: patches=-1"
WRITE(6,'(A,I0,A)') "CONTAINS"
WRITE(6,'(A,I0,A)') "SUBROUTINE T_C_G0_n(RES,use_gamma,R,T,NDERIV)"
WRITE(6,'(A)') " IMPLICIT NONE"
WRITE(6,'(A)') " REAL(KIND=dp), INTENT(OUT) :: RES(*)"
WRITE(6,'(A)') " LOGICAL, INTENT(OUT) :: use_gamma"
WRITE(6,'(A)') " REAL(KIND=dp),INTENT(IN) :: R,T"
WRITE(6,'(A)') " INTEGER, INTENT(IN) :: NDERIV"
WRITE(6,'(A)') " REAL(KIND=dp) :: upper,lower,X1,X2,TG1,TG2"
WRITE(6,'(A)') " use_gamma=.FALSE."
WRITE(6,'(A)') " upper=R**2 + 11.0_dp*R + 50.0_dp"
WRITE(6,'(A)') " lower=R**2 - 11.0_dp*R + 0.0_dp"
WRITE(6,'(A)') " IF (T>upper) THEN"
WRITE(6,'(A)') " RES(1:NDERIV+1)=0.0_dp"
WRITE(6,'(A)') " RETURN"
WRITE(6,'(A)') " ENDIF"
WRITE(6,'(A)') " IF (R<=11.0_dp) THEN"
WRITE(6,'(A)') " X2=R/11.0_dp"
WRITE(6,'(A)') " upper=R**2 + 11.0_dp*R + 50.0_dp"
WRITE(6,'(A)') " lower=0.0_dp"
WRITE(6,'(A)') " X1=(T-lower)/(upper-lower)"
functionid=2
open(UNIT=17,FILE="nearfield_domains.dat")
A1=0.0D0
B1=1.0D0
A2=0.0D0
B2=1.0D0
CALL FIND_DOMAINS(A1,B1,A2,B2,DEG,NDERIV,ERR,LEVEL)
close(17)
WRITE(6,'(A)') " ELSE"
WRITE(6,'(A)') " IF (T<lower) THEN"
WRITE(6,'(A)') " use_gamma=.TRUE."
WRITE(6,'(A)') " RETURN"
WRITE(6,'(A)') " ENDIF"
WRITE(6,'(A)') " X2=11.0_dp/R"
WRITE(6,'(A)') " X1=(T-lower)/(upper-lower)"
functionid=1
open(UNIT=17,FILE="farfield_domains.dat")
A1=0.0D0
B1=1.0D0
A2=0.0D0
B2=1.0D0
CALL FIND_DOMAINS(A1,B1,A2,B2,DEG,NDERIV,ERR,LEVEL)
close(17)
WRITE(6,'(A)') " ENDIF"
WRITE(6,'(A,I0,A)') "END SUBROUTINE T_C_G0_n"
CLOSE(UNIT=13)
CALL WRITE_INIT()
CALL WRITE_PD2VAL(DEG)
WRITE(6,'(A)') "END MODULE T_C_G0"
END SUBROUTINE GO
SUBROUTINE WRITE_INIT()
WRITE(6,'(A)') "! iunit contains the data file to initialize the table"
WRITE(6,'(A)') "! Nder is the number of derivatives that will actually be used"
WRITE(6,'(A)') "SUBROUTINE INIT(Nder,iunit)"
WRITE(6,'(A)') " IMPLICIT NONE"
WRITE(6,'(A)') " INTEGER, INTENT(IN) :: Nder,iunit"
WRITE(6,'(A)') " INTEGER :: I"
WRITE(6,'(A)') " REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: chunk"
WRITE(6,'(A,I0)') " patches=",find_count
WRITE(6,'(A)') ' IF (Nder>nderiv_max) STOP "T_C_G0 init failed"'
WRITE(6,'(A)') " nderiv_init=Nder"
WRITE(6,'(A)') " IF(ALLOCATED(C0)) DEALLOCATE(C0)"
WRITE(6,'(A)') " ! round up to a multiple of 32 to give some generous alignment for each C0"
WRITE(6,'(A)') " ALLOCATE(C0(32*((31+(Nder+1)*(degree+1)*(degree+2)/2)/32),patches))"
WRITE(6,'(A)') " ALLOCATE(chunk((nderiv_max+1)*(degree+1)*(degree+2)/2))"
WRITE(6,'(A)') " DO I=1,patches"
WRITE(6,'(A)') " READ(iunit,*) chunk"
WRITE(6,'(A)') " C0(1:(Nder+1)*(degree+1)*(degree+2)/2,I)=chunk(1:(Nder+1)*(degree+1)*(degree+2)/2)"
WRITE(6,'(A)') " ENDDO"
WRITE(6,'(A)') " DEALLOCATE(chunk)"
WRITE(6,'(A)') "END SUBROUTINE INIT"
WRITE(6,'(A)') "SUBROUTINE FREE()"
WRITE(6,'(A)') " IF(ALLOCATED(C0)) DEALLOCATE(C0)"
WRITE(6,'(A)') " nderiv_init=-1"
WRITE(6,'(A)') "END SUBROUTINE FREE"
END SUBROUTINE WRITE_INIT
SUBROUTINE WRITE_PD2VAL(DEG)
INTEGER :: DEG
INTEGER :: I,J
WRITE(6,'(A)') "SUBROUTINE PD2VAL(RES,NDERIV,TG1,TG2,C0) "
WRITE(6,'(A)') " IMPLICIT NONE"
WRITE(6,'(A)') " REAL(KIND=dp), INTENT(OUT) :: res(*)"
WRITE(6,'(A)') " INTEGER, INTENT(IN) :: NDERIV"
WRITE(6,'(A)') " REAL(KIND=dp),INTENT(IN) :: TG1,TG2"
WRITE(6,'(A,I0,A)') " REAL(KIND=dp) :: T1(0:",DEG,")"
WRITE(6,'(A,I0,A)') " REAL(KIND=dp) :: T2(0:",DEG,")"
WRITE(6,'(A,I0,A)') " REAL(KIND=dp), INTENT(IN) :: C0(",(DEG+1)*(DEG+2)/2,",*)"
WRITE(6,'(A,I0,A)') " INTEGER :: I,J,K"
WRITE(6,'(A,I0,A)') " REAL(KIND=dp), PARAMETER :: SQRT2=1.4142135623730950488016887242096980785696718753_dp"
WRITE(6,'(A,I0,A)') " T1(0)=1.0_dp"
WRITE(6,'(A,I0,A)') " T2(0)=1.0_dp"
WRITE(6,'(A,I0,A)') " T1(1)=SQRT2 * TG1"
WRITE(6,'(A,I0,A)') " T2(1)=SQRT2 * TG2"
WRITE(6,'(A,I0,A)') " T1(2)= 2 * TG1 * T1(1) - SQRT2 "
WRITE(6,'(A,I0,A)') " T2(2)= 2 * TG2 * T2(1) - SQRT2 "
DO I=3,DEG
WRITE(6,'(A,I0,A,I0,A,I0,A)') " T1(",I,") = 2 * TG1 * T1(",I-1,") - T1(",I-2,")"
WRITE(6,'(A,I0,A,I0,A,I0,A)') " T2(",I,") = 2 * TG2 * T2(",I-1,") - T2(",I-2,")"
ENDDO
WRITE(6,'(A)') " DO K=1,NDERIV+1"
WRITE(6,'(A)') " RES(K) = 0.0_dp"
J=1
DO I=DEG,0,-1
WRITE(6,'(A,I0,A,I0,A,I0,A,I0,A)') " RES(K)=RES(K)+DOT_PRODUCT(T1(0:",I,"),C0(",J,":",J+I,",K))*T2(",DEG-I,")"
J=J+I+1
ENDDO
WRITE(6,'(A)') " ENDDO"
WRITE(6,'(A)') "END SUBROUTINE PD2VAL"
END SUBROUTINE
!
! the basic piece of the code, does the bisection till convergence
!
RECURSIVE SUBROUTINE FIND_DOMAINS(A1,B1,A2,B2,DEG,NDERIV,ERR,LEVEL)
IMPLICIT NONE
REAL(KIND=dp) A1,B1,A2,B2,ERR,R1,R2
INTEGER :: DEG,LEVEL,I,DIRECTION,NDERIV,K,I1,I2
REAL(KIND=dp) ESTIMATED_ERR,MIDDLE,ESTIMATED_ERR1,ESTIMATED_ERR2,ESTIMATED_ERR3,ESTIMATED_ERR4
REAL(KIND=dp) :: C0(0:DEG,0:DEG,0:NDERIV)
IF (A1==B1 .AND. A2==B2) STOP "BUG 1"
IF (LEVEL>20) STOP "Too many bisections "
! compute the interpolating polynomial and its estimated error
CALL drvrec_wrap(A1,B1,A2,B2,DEG,NDERIV,C0,ESTIMATED_ERR)
IF (ESTIMATED_ERR<ERR) THEN
! we have converged ... turn this into a call to pd2val, and append to the C0 table
find_count=find_count+1
WRITE(6,'(A,E25.18,A,E25.18,A)') REPEAT(" ",LEVEL+3)//"TG1= ( 2 * X1 - ",(B1+A1),"_dp)*",1.0D0/(B1-A1),"_dp"
WRITE(6,'(A,E25.18,A,E25.18,A)') REPEAT(" ",LEVEL+3)//"TG2= ( 2 * X2 - ",(B2+A2),"_dp)*",1.0D0/(B2-A2),"_dp"
WRITE(6,'(A,I0,A)') REPEAT(" ",LEVEL+3)//"CALL PD2VAL(RES,NDERIV,TG1,TG2,C0(1,", find_count, "))"
WRITE(17,'(2F20.16,I10,E16.8)') A1,A2,find_count,ESTIMATED_ERR
WRITE(17,'(2F20.16,I10,E16.8)') A1,B2 ,find_count,ESTIMATED_ERR
WRITE(17,*)
WRITE(17,'(2F20.16,I10,E16.8)') A1,A2 ,find_count,ESTIMATED_ERR
WRITE(17,'(2F20.16,I10,E16.8)') B1,A2 ,find_count,ESTIMATED_ERR
WRITE(17,*)
WRITE(17,'(2F20.16,I10,E16.8)') B1,B2,find_count,ESTIMATED_ERR
WRITE(17,'(2F20.16,I10,E16.8)') A1,B2,find_count,ESTIMATED_ERR
WRITE(17,*)
WRITE(17,'(2F20.16,I10,E16.8)') B1,B2,find_count,ESTIMATED_ERR
WRITE(17,'(2F20.16,I10,E16.8)') B1,A2,find_count,ESTIMATED_ERR
WRITE(17,*)
DO K=0,NDERIV
DO I=DEG,0,-1
write(13,*) C0(0:I,DEG-I,K)
ENDDO
ENDDO
ELSE
! decide in which dimension to bisection.
! using a section that minimizes the maximum error is slightly more
! efficient than alternating bisection.
MIDDLE=(A1+B1)/2
CALL drvrec_wrap(A1,MIDDLE,A2,B2,DEG,NDERIV,C0,ESTIMATED_ERR1)
CALL drvrec_wrap(MIDDLE,B1,A2,B2,DEG,NDERIV,C0,ESTIMATED_ERR2)
MIDDLE=(A2+B2)/2
CALL drvrec_wrap(A1,B1,A2,MIDDLE,DEG,NDERIV,C0,ESTIMATED_ERR3)
CALL drvrec_wrap(A1,B1,MIDDLE,B2,DEG,NDERIV,C0,ESTIMATED_ERR4)
IF (MAX(ESTIMATED_ERR1,ESTIMATED_ERR2)<MAX(ESTIMATED_ERR3,ESTIMATED_ERR4)) THEN
DIRECTION=1
ELSE
DIRECTION=2
ENDIF
IF (DIRECTION==1) THEN
MIDDLE=(A1+B1)/2
WRITE(6,'(A,E25.18,A)') REPEAT(" ",LEVEL)//" IF (X1<=",MIDDLE,"_dp) THEN "
CALL FIND_DOMAINS(A1,MIDDLE,A2,B2,DEG,NDERIV,ERR,LEVEL+1)
WRITE(6,'(A)') REPEAT(" ",LEVEL)//" ELSE"
CALL FIND_DOMAINS(MIDDLE,B1,A2,B2,DEG,NDERIV,ERR,LEVEL+1)
WRITE(6,'(A)') REPEAT(" ",LEVEL)//" ENDIF"
ELSE
MIDDLE=(A2+B2)/2
WRITE(6,'(A,E25.18,A)') REPEAT(" ",LEVEL)//" IF (X2<=",MIDDLE,"_dp) THEN "
CALL FIND_DOMAINS(A1,B1,A2,MIDDLE,DEG,NDERIV,ERR,LEVEL+1)
WRITE(6,'(A)') REPEAT(" ",LEVEL)//" ELSE"
CALL FIND_DOMAINS(A1,B1,MIDDLE,B2,DEG,NDERIV,ERR,LEVEL+1)
WRITE(6,'(A)') REPEAT(" ",LEVEL)//" ENDIF"
ENDIF
ENDIF
END SUBROUTINE FIND_DOMAINS
SUBROUTINE drvrec_wrap(A1,B1,A2,B2,DEG,NDERIV,C0,ESTIMATED_ERR)
REAL(KIND=dp) :: a1,b1,a2,b2,estimated_err,estimated_err2,estimated_err3,esterr
INTEGER :: deg,nderiv
REAL(KIND=dp) :: C0(0:DEG,0:DEG,0:NDERIV), C02(0:DEG,0:DEG,0:NDERIV)
integer :: functionid,digits
COMMON/functiontype/functionid,digits
OPEN(UNIT=51,FILE="drvrec_mpfr.in")
WRITE(51,*) digits
WRITE(51,*) deg
WRITE(51,*) nderiv
WRITE(51,*) A1,B1
WRITE(51,*) A2,B2
WRITE(51,*) functionid
CLOSE(51)
CALL SYSTEM("./drvrec_mpfr.x")
OPEN(UNIT=51,FILE="drvrec_mpfr.out")
READ(51,*) ESTIMATED_ERR
READ(51,*) C0
CLOSE(51)
IF (.FALSE.) THEN
CALL drvrec(A1,B1,A2,B2,DEG,NDERIV,C02,ESTIMATED_ERR2)
estimated_err3=0
DO K=0,NDERIV
ESTERR = 0.0D0
DO J = 0,2
DO I = 0,DEG - J
ESTERR = ESTERR + ABS(C0(I,DEG - I - J,K))
ENDDO
ENDDO
ESTERR = ESTERR * 2
estimated_err3=MAX(ESTERR,estimated_err3)
ENDDO
write(100,*) A1,B1,A2,B2,estimated_err,estimated_err2,estimated_err3,MAXVAL(ABS(C0-C02))
ENDIF
END SUBROUTINE
END MODULE
!
! The main program ... small and easy
!
USE RECURSE
CALL GO()
END

68
tools/Fun2D/tester.f90 Normal file
View file

@ -0,0 +1,68 @@
USE T_C_G0
IMPLICIT NONE
INTEGER :: Nder
REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: ref
REAL(KIND=dp) :: t1,t2,R,T
REAL(KIND=dp), ALLOCATABLE :: RES(:),MAXERR(:)
INTEGER :: I,J,Nline
LOGICAL :: use_gamma
Nder=MIN(9,nderiv_max)
OPEN(12,FILE="T_C_G.dat")
CALL INIT(Nder,12)
CLOSE(12)
ALLOCATE(RES(0:Nder),MAXERR(0:Nder))
Nline=0
OPEN(13,FILE="history.dat")
DO
READ(13,*,END=999) R,T
Nline=Nline+1
ENDDO
999 CONTINUE
REWIND(13)
ALLOCATE(ref(-2:Nder,Nline))
DO I=1,Nline
READ(13,*) ref(:,I)
ENDDO
CLOSE(13)
maxerr=0
write(6,'(A)') "List of points with largest errors so far"
write(6,'(A)') "-----------------------------------------"
CALL CPU_TIME(T1)
! reference points can be outside the expected domain due to
! the fact that they are rounded from arbitrary precision to
! double precision
DO I=1,Nline
R=ref(-2,I)
T=ref(-1,I)
CALL T_C_G0_n(RES,use_gamma,R,T,Nder)
IF (.NOT.use_gamma) THEN
IF (ANY(RES.NE.0)) THEN
IF (ANY((abs(RES-REF(0:Nder,I)))>(maxerr))) THEN
maxerr=MAX(maxerr,abs(RES-REF(0:Nder,I)))
write(6,*) R,T,RES(0:Nder),REF(0:Nder,I)
ENDIF
ENDIF
ELSE
! write(6,*) "oops",R,T,R**2-11*R-T
ENDIF
ENDDO
CALL CPU_TIME(T2)
write(6,'(A)') "Summary of generation"
write(6,'(A)') "---------------------"
write(6,'(A,I8)') "Degree of the interpolation: ",degree
write(6,'(A,E8.3)') "Target absolute error: ",target_error
write(6,'(A,I8)') "Number of domains after bisection: ",patches
write(6,'(A,I8)') "Table size for polynomial coefficients [Kb]: ",(8*patches*(Nder+1)*(degree+1)*(degree+2)/2)/(1024)
write(6,'(A,I8)') "Number of reference points: ",Nline
write(6,'(A,I8)') "Tested up to derivative: ",Nder
write(6,'(A,F8.3)') "Time for evaluation [s]: ",t2-t1
write(6,'(A,E8.3)') "Maximum error over all reference points: ",MAXVAL(maxerr)
write(6,'(A)') "Maximum error per derivative:"
write(6,'(5E18.4)') maxerr
END