1!
2!     a simple wrapper program to
3!     obtain the C0 coefficients for an NDERIV+1 valued function in the
4!     rectangular domain A1<=X1<=B1, A2<=X2<=B2
5!
6!     input from 'drvrec_mpfr.in'
7!     output to 'drvrec_mpfr.out'
8!
9
10      SUBROUTINE drvrec(a1,b1,a2,b2,deg,nderiv,c0,esterr_out)
11        USE padua2_mpfr
12        USE functions, ONLY: evalf=>f
13        IMPLICIT NONE
14
15        INTEGER deg,npd,nderiV
16        TYPE(mpfr_type) pd1((deg + 1) * (deg + 2) / 2)
17        TYPE(mpfr_type) pd2((deg + 1) * (deg + 2) / 2)
18        TYPE(mpfr_type) wpd((deg + 1) * (deg + 2) / 2)
19        TYPE(mpfr_type) fpd((deg + 1) * (deg + 2) / 2,0:nderiv)
20        TYPE(mpfr_type) t1(0:deg,(deg + 1) * (deg + 2) / 2)
21        TYPE(mpfr_type) t2(0:deg,(deg + 1) * (deg + 2) / 2)
22        TYPE(mpfr_type) c0(0:deg,0:deg,0:nderiv)
23        TYPE(mpfr_type) res(0:nderiv),esterr(0:nderiv)
24        TYPE(mpfr_type) a1,b1,a2,b2,esterr_out,x1,x2,zero
25        INTEGER i,ntg1
26
27        CALL pdpts(deg,pd1,pd2,wpd,npd)
28
29!$OMP PARALLEL DO &
30!$OMP DEFAULT(NONE) &
31!$OMP PRIVATE(i,x1,x2,res) &
32!$OMP SHARED(npd,a1,b1,a2,b2,pd1,pd2,fpd,nderiv)
33        DO i = 1,npd
34           x1=((b1 - a1) * pd1(i) + (b1 + a1)) / 2
35           x2=((b2 - a2) * pd2(i) + (b2 + a2)) / 2
36           CALL evalf(res,nderiv,x1,x2)
37           FPD(I,:) = RES(0:NDERIV)
38        ENDDO
39
40!$OMP PARALLEL DO &
41!$OMP DEFAULT(NONE) &
42!$OMP PRIVATE(i,t1,t2) &
43!$OMP SHARED(npd,deg,pd1,pd2,wpd,fpd,c0,esterr)
44        DO i=0,nderiv
45           CALL padua2(deg,npd,pd1,pd2,wpd,fpd(:,i),t1,t2,c0(:,:,i),esterr(i))
46        ENDDO
47
48        esterr_out=0
49        zero=0
50        DO i=0,nderiv
51           IF (mpfr_cmp(esterr(i),zero)>0) THEN
52              IF (mpfr_cmp(esterr(i),esterr_out)>0) THEN
53                   esterr_out=esterr(i)
54              ENDIF
55           ELSE
56              IF (mpfr_cmp(-esterr(i),esterr_out)>0) esterr_out=-esterr(i)
57           ENDIF
58        ENDDO
59
60      END SUBROUTINE drvrec
61
62      !
63      ! this program is a simple interface to obtain the C0 coefficients
64      ! with high accuracy (i.e. all digits & correctly rounded),
65      ! for a function programmed in subroutine f
66      ! this version uses double precision (~16 digits)
67      ! input/output but does all calculations
68      ! internally using a user specified number of digits
69      !
70      PROGRAM drvrec_padua
71         USE mpfr
72         USE mpfr_ops
73         USE padua2_mpfr, ONLY: pd2val_mp=>pd2val
74         USE functions, ONLY: functionid,should_output
75         IMPLICIT NONE
76         INTEGER :: deg,nderiv,i,j,digits,deg_ref,k,l
77         INTEGER*2 :: precision
78
79         REAL(KIND=dp) :: esterr,A1,B1,A2,B2
80         REAL(KIND=dp), DIMENSION(:,:,:), ALLOCATABLE :: c0
81
82         TYPE(mpfr_type) ::  A1_mp,B1_mp,A2_mp,B2_mp,esterr_mp
83         TYPE(mpfr_type), DIMENSION(:,:,:), ALLOCATABLE :: c0_mp
84
85         ! keeps a history of evaluated function values.
86         ! useful to investigate the error later on
87         should_output=31
88         IF (should_output>0) THEN
89            OPEN(should_output,FILE="function_vals.dat")
90         ENDIF
91
92         ! input definition
93         OPEN(17,FILE="drvrec_mpfr.in")
94         READ(17,*) digits
95         READ(17,*) deg
96         READ(17,*) nderiv
97         READ(17,*) A1,B1
98         READ(17,*) A2,B2
99         READ(17,*) functionid
100         CLOSE(17)
101
102         precision = log(10.0)/log(2.0)*digits
103!$OMP PARALLEL DEFAULT(NONE) SHARED(precision)
104         CALL mpfr_set_default_precision(precision)
105!$OMP END PARALLEL
106
107         A1_mp=A1
108         B1_mp=B1
109         A2_mp=A2
110         B2_mp=B2
111
112         ALLOCATE(c0_mp(0:deg,0:deg,0:nderiv))
113         CALL drvrec(a1_mp,b1_mp,a2_mp,b2_mp,deg,nderiv,c0_mp,esterr_mp)
114
115         ! convert ot double precision and output
116         ALLOCATE(c0(0:deg,0:deg,0:nderiv))
117         c0=REAL(c0_mp)
118         esterr=REAL(esterr_mp)
119
120         OPEN(17,FILE="drvrec_mpfr.out")
121         write(17,*) esterr
122         write(17,*) c0
123         CLOSE(17)
124
125         IF (should_output>0) THEN
126           CLOSE(should_output)
127         ENDIF
128
129      END PROGRAM drvrec_padua
130