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