C Does not work: f2py -m mvnprdmod -h mvnprdmod.pyf mvnprodcorrprb.f only: mvnprodcorrprb
C gfortran -fPIC -c mvnprodcorrprb.f
C f2py -m mvnprdmod -c mvnprodcorrprb.o mvnprodcorrprb_interface.f --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71
C f2py -m mvnprdmod -c mvnprodcorrprb.o mvnprodcorrprb_interface.f --build-dir tmp1 --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71
* This is a MEX-file for MATLAB.
* and contains a mex-interface to, mvnprodcorrprb a subroutine
* for computing multivariate normal probabilities with product
* correlation structure.
* The file should compile without errors on (Fortran90)
* standard Fortran compilers.
*
* The mex-interface and mvnprodcorrprb was written by
* Per Andreas Brodtkorb
* Norwegian Defence Research Establishment
* P.O. Box 115m
* N-3191 Horten
* Norway
* Email: Per.Brodtkorb@ffi.no
*
*
* MVNPRODCORRPRBMEX Computes multivariate normal probability
* with product correlation structure.
*
* CALL [value,error,inform]=mvnprodcorrprbmex(rho,A,B,abseps,releps,useBreakPoints);
*
* RHO REAL, array of coefficients defining the correlation
* coefficient by:
* correlation(I,J) = RHO(I)*RHO(J) for J/=I
* where
* 1 <= RHO(I) <= 1
* A REAL, array of lower integration limits.
* B REAL, array of upper integration limits.
* NOTE: any values greater the 10, are considered as
* infinite values.
* ABSEPS REAL absolute error tolerance.
* RELEPS REAL relative error tolerance.
* USEBREAKPOINTS = 1 If extra integration points should be used
* around possible singularities
* 0 If no extra
*
* ERROR REAL estimated absolute error, with 99% confidence level.
* VALUE REAL estimated value for the integral
* INFORM INTEGER, termination status parameter:
* if INFORM = 0, normal completion with ERROR < EPS;
* if INFORM = 1, completion with ERROR > EPS and MAXPTS
* function vaules used; increase MAXPTS to
* decrease ERROR;
*
* MVNPRODCORRPRB calculates multivariate normal probability
* with product correlation structure for rectangular regions.
* The accuracy is up to almost double precision, i.e., about 1e-14.
*
* This file was successfully compiled for matlab 5.3
* using Compaq Visual Fortran 6.1, and Windows 2000.
* The example here uses Fortran77 source.
* First, you will need to modify your mexopts.bat file.
* To find it, issue the command prefdir(1) from the Matlab command line,
* the directory it answers with will contain your mexopts.bat file.
* Open it for editing. The first section will look like:
*
*rem ********************************************************************
*rem General parameters
*rem ********************************************************************
*set MATLAB=%MATLAB%
*set DF_ROOT=C:\Program Files\Microsoft Visual Studio
*set VCDir=%DF_ROOT%\VC98
*set MSDevDir=%DF_ROOT%\Common\msdev98
*set DFDir=%DF_ROOT%\DF98
*set PATH=%MSDevDir%\bin;%DFDir%\BIN;%VCDir%\BIN;%PATH%
*set INCLUDE=%DFDir%\INCLUDE;%DFDir%\IMSL\INCLUDE;%INCLUDE%
*set LIB=%DFDir%\LIB;%VCDir%\LIB
*
* then you are ready to compile this file at the matlab prompt using the
* following command:
* mex -O mvnprodcorrprbmex.f
MODULE ERFCOREMOD
IMPLICIT NONE
INTERFACE CALERF
MODULE PROCEDURE CALERF
END INTERFACE
INTERFACE DERF
MODULE PROCEDURE DERF
END INTERFACE
INTERFACE DERFC
MODULE PROCEDURE DERFC
END INTERFACE
INTERFACE DERFCX
MODULE PROCEDURE DERFCX
END INTERFACE
CONTAINS
C--------------------------------------------------------------------
C
C DERF subprogram computes approximate values for erf(x).
C (see comments heading CALERF).
C
C Author/date: W. J. Cody, January 8, 1985
C
C--------------------------------------------------------------------
FUNCTION DERF( X ) RESULT (VALUE)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: X
DOUBLE PRECISION :: VALUE
INTEGER, PARAMETER :: JINT = 0
CALL CALERF(X,VALUE,JINT)
RETURN
END FUNCTION DERF
C--------------------------------------------------------------------
C
C DERFC subprogram computes approximate values for erfc(x).
C (see comments heading CALERF).
C
C Author/date: W. J. Cody, January 8, 1985
C
C--------------------------------------------------------------------
FUNCTION DERFC( X ) RESULT (VALUE)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: X
DOUBLE PRECISION :: VALUE
INTEGER, PARAMETER :: JINT = 1
CALL CALERF(X,VALUE,JINT)
RETURN
END FUNCTION DERFC
C------------------------------------------------------------------
C
C DERFCX subprogram computes approximate values for exp(x*x) * erfc(x).
C (see comments heading CALERF).
C
C Author/date: W. J. Cody, March 30, 1987
C
C------------------------------------------------------------------
FUNCTION DERFCX( X ) RESULT (VALUE)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: X
DOUBLE PRECISION :: VALUE
INTEGER, PARAMETER :: JINT = 2
CALL CALERF(X,VALUE,JINT)
RETURN
END FUNCTION DERFCX
SUBROUTINE CALERF(ARG,RESULT,JINT)
IMPLICIT NONE
C------------------------------------------------------------------
C
C CALERF packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x)
C for a real argument x. It contains three FUNCTION type
C subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX),
C and one SUBROUTINE type subprogram, CALERF. The calling
C statements for the primary entries are:
C
C Y=ERF(X) (or Y=DERF(X)),
C
C Y=ERFC(X) (or Y=DERFC(X)),
C and
C Y=ERFCX(X) (or Y=DERFCX(X)).
C
C The routine CALERF is intended for internal packet use only,
C all computations within the packet being concentrated in this
C routine. The function subprograms invoke CALERF with the
C statement
C
C CALL CALERF(ARG,RESULT,JINT)
C
C where the parameter usage is as follows
C
C Function Parameters for CALERF
C call ARG Result JINT
C
C ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0
C ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1
C ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2
C
C The main computation evaluates near-minimax approximations
C from "Rational Chebyshev approximations for the error function"
C by W. J. Cody, Math. Comp., 1969, PP. 631-638. This
C transportable program uses rational functions that theoretically
C approximate erf(x) and erfc(x) to at least 18 significant
C decimal digits. The accuracy achieved depends on the arithmetic
C system, the compiler, the intrinsic functions, and proper
C selection of the machine-dependent constants.
C
C*******************************************************************
C*******************************************************************
C
C Explanation of machine-dependent constants
C
C XMIN = the smallest positive floating-point number.
C XINF = the largest positive finite floating-point number.
C XNEG = the largest negative argument acceptable to ERFCX;
C the negative of the solution to the equation
C 2*exp(x*x) = XINF.
C XSMALL = argument below which erf(x) may be represented by
C 2*x/sqrt(pi) and above which x*x will not underflow.
C A conservative value is the largest machine number X
C such that 1.0 + X = 1.0 to machine precision.
C XBIG = largest argument acceptable to ERFC; solution to
C the equation: W(x) * (1-0.5/x**2) = XMIN, where
C W(x) = exp(-x*x)/[x*sqrt(pi)].
C XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to
C machine precision. A conservative value is
C 1/[2*sqrt(XSMALL)]
C XMAX = largest acceptable argument to ERFCX; the minimum
C of XINF and 1/[sqrt(pi)*XMIN].
C
C Approximate values for some important machines are:
C
C XMIN XINF XNEG XSMALL
C
C C 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15
C CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15
C IEEE (IBM/XT,
C SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8
C IEEE (IBM/XT,
C SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16
C IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17
C UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18
C VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17
C VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16
C
C
C XBIG XHUGE XMAX
C
C C 7600 (S.P.) 25.922 8.39E+6 1.80X+293
C CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465
C IEEE (IBM/XT,
C SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37
C IEEE (IBM/XT,
C SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307
C IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75
C UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307
C VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38
C VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307
C
C*******************************************************************
C*******************************************************************
C
C Error returns
C
C The program returns ERFC = 0 for ARG .GE. XBIG;
C
C ERFCX = XINF for ARG .LT. XNEG;
C and
C ERFCX = 0 for ARG .GE. XMAX.
C
C
C Intrinsic functions required are:
C
C ABS, AINT, EXP
C
C
C Author: W. J. Cody
C Mathematics and Computer Science Division
C Argonne National Laboratory
C Argonne, IL 60439
C
C Latest modification: March 19, 1990
C Updated to F90 by pab 23.03.2003
C Revised pab Dec 2008
C updated parameter statements in CALERF so that it works when
C compiling with gfortran.
C
C------------------------------------------------------------------
DOUBLE PRECISION, INTENT(IN) :: ARG
INTEGER, INTENT(IN) :: JINT
DOUBLE PRECISION, INTENT(INOUT):: RESULT
! Local variables
INTEGER :: I
DOUBLE PRECISION :: DEL,X,XDEN,XNUM,Y,YSQ
C------------------------------------------------------------------
C Mathematical constants
C------------------------------------------------------------------
DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0
DOUBLE PRECISION, PARAMETER :: HALF = 0.05D0
DOUBLE PRECISION, PARAMETER :: ONE = 1.0D0
DOUBLE PRECISION, PARAMETER :: TWO = 2.0D0
DOUBLE PRECISION, PARAMETER :: FOUR = 4.0D0
DOUBLE PRECISION, PARAMETER :: SIXTEN = 16.0D0
DOUBLE PRECISION, PARAMETER :: SQRPI = 5.6418958354775628695D-1
DOUBLE PRECISION, PARAMETER :: THRESH = 0.46875D0
C------------------------------------------------------------------
C Machine-dependent constants
C------------------------------------------------------------------
DOUBLE PRECISION, PARAMETER :: XNEG = -26.628D0
DOUBLE PRECISION, PARAMETER :: XSMALL = 1.11D-16
DOUBLE PRECISION, PARAMETER :: XBIG = 26.543D0
DOUBLE PRECISION, PARAMETER :: XHUGE = 6.71D7
DOUBLE PRECISION, PARAMETER :: XMAX = 2.53D307
DOUBLE PRECISION, PARAMETER :: XINF = 1.79D308
!---------------------------------------------------------------
! Coefficents to the rational polynomials
!--------------------------------------------------------------
C DOUBLE PRECISION, DIMENSION(5) :: A, Q
C DOUBLE PRECISION, DIMENSION(4) :: B
C DOUBLE PRECISION, DIMENSION(9) :: C
C DOUBLE PRECISION, DIMENSION(8) :: D
C DOUBLE PRECISION, DIMENSION(6) :: P
C------------------------------------------------------------------
C Coefficients for approximation to erf in first interval
C------------------------------------------------------------------
DOUBLE PRECISION, PARAMETER, DIMENSION(5) ::
& A = (/ 3.16112374387056560D00,
& 1.13864154151050156D02,3.77485237685302021D02,
& 3.20937758913846947D03, 1.85777706184603153D-1/)
DOUBLE PRECISION, PARAMETER, DIMENSION(4) ::
& B = (/2.36012909523441209D01,2.44024637934444173D02,
& 1.28261652607737228D03,2.84423683343917062D03/)
C------------------------------------------------------------------
C Coefficients for approximation to erfc in second interval
C------------------------------------------------------------------
DOUBLE PRECISION, DIMENSION(9) ::
& C=(/5.64188496988670089D-1,8.88314979438837594D0,
1 6.61191906371416295D01,2.98635138197400131D02,
2 8.81952221241769090D02,1.71204761263407058D03,
3 2.05107837782607147D03,1.23033935479799725D03,
4 2.15311535474403846D-8/)
DOUBLE PRECISION, DIMENSION(8) ::
& D =(/1.57449261107098347D01,1.17693950891312499D02,
1 5.37181101862009858D02,1.62138957456669019D03,
2 3.29079923573345963D03,4.36261909014324716D03,
3 3.43936767414372164D03,1.23033935480374942D03/)
C------------------------------------------------------------------
C Coefficients for approximation to erfc in third interval
C------------------------------------------------------------------
DOUBLE PRECISION, parameter,
& DIMENSION(6) :: P =(/3.05326634961232344D-1,
& 3.60344899949804439D-1,
1 1.25781726111229246D-1,1.60837851487422766D-2,
2 6.58749161529837803D-4,1.63153871373020978D-2/)
DOUBLE PRECISION, parameter,
& DIMENSION(5) :: Q =(/2.56852019228982242D00,
& 1.87295284992346047D00,
1 5.27905102951428412D-1,6.05183413124413191D-2,
2 2.33520497626869185D-3/)
C------------------------------------------------------------------
X = ARG
Y = ABS(X)
IF (Y .LE. THRESH) THEN
C------------------------------------------------------------------
C Evaluate erf for |X| <= 0.46875
C------------------------------------------------------------------
YSQ = ZERO
IF (Y .GT. XSMALL) THEN
YSQ = Y * Y
XNUM = A(5)*YSQ
XDEN = YSQ
DO I = 1, 3
XNUM = (XNUM + A(I)) * YSQ
XDEN = (XDEN + B(I)) * YSQ
END DO
RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
ELSE
RESULT = X * A(4) / B(4)
ENDIF
IF (JINT .NE. 0) RESULT = ONE - RESULT
IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
GO TO 800
C------------------------------------------------------------------
C Evaluate erfc for 0.46875 <= |X| <= 4.0
C------------------------------------------------------------------
ELSE IF (Y .LE. FOUR) THEN
XNUM = C(9)*Y
XDEN = Y
DO I = 1, 7
XNUM = (XNUM + C(I)) * Y
XDEN = (XDEN + D(I)) * Y
END DO
RESULT = (XNUM + C(8)) / (XDEN + D(8))
IF (JINT .NE. 2) THEN
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
END IF
C------------------------------------------------------------------
C Evaluate erfc for |X| > 4.0
C------------------------------------------------------------------
ELSE
RESULT = ZERO
IF (Y .GE. XBIG) THEN
IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300
IF (Y .GE. XHUGE) THEN
RESULT = SQRPI / Y
GO TO 300
END IF
END IF
YSQ = ONE / (Y * Y)
XNUM = P(6)*YSQ
XDEN = YSQ
DO I = 1, 4
XNUM = (XNUM + P(I)) * YSQ
XDEN = (XDEN + Q(I)) * YSQ
ENDDO
RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
RESULT = (SQRPI - RESULT) / Y
IF (JINT .NE. 2) THEN
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
END IF
END IF
C------------------------------------------------------------------
C Fix up for negative argument, erf, etc.
C------------------------------------------------------------------
300 IF (JINT .EQ. 0) THEN
RESULT = (HALF - RESULT) + HALF
IF (X .LT. ZERO) RESULT = -RESULT
ELSE IF (JINT .EQ. 1) THEN
IF (X .LT. ZERO) RESULT = TWO - RESULT
ELSE
IF (X .LT. ZERO) THEN
IF (X .LT. XNEG) THEN
RESULT = XINF
ELSE
YSQ = AINT(X*SIXTEN)/SIXTEN
DEL = (X-YSQ)*(X+YSQ)
Y = EXP(YSQ*YSQ) * EXP(DEL)
RESULT = (Y+Y) - RESULT
END IF
END IF
END IF
800 RETURN
END SUBROUTINE CALERF
END MODULE ERFCOREMOD
module functionInterface
INTERFACE
FUNCTION F(Z) result (VAL)
DOUBLE PRECISION, INTENT(IN) :: Z
DOUBLE PRECISION :: VAL
END FUNCTION F
END INTERFACE
end module functionInterface
module AdaptiveGaussKronrod
implicit none
private
public :: dqagpe,dqagp
INTERFACE dqagpe
MODULE PROCEDURE dqagpe
END INTERFACE
INTERFACE dqagp
MODULE PROCEDURE dqagp
END INTERFACE
INTERFACE dqelg
MODULE PROCEDURE dqelg
END INTERFACE
INTERFACE dqpsrt
MODULE PROCEDURE dqpsrt
END INTERFACE
INTERFACE dqk21
MODULE PROCEDURE dqk21
END INTERFACE
INTERFACE dqk15
MODULE PROCEDURE dqk15
END INTERFACE
INTERFACE dqk9
MODULE PROCEDURE dqk9
END INTERFACE
INTERFACE d1mach
MODULE PROCEDURE d1mach
END INTERFACE
contains
subroutine dea3(E0,E1,E2,abserr,result)
!***PURPOSE Given a slowly convergent sequence, this routine attempts
! to extrapolate nonlinearly to a better estimate of the
! sequence's limiting value, thus improving the rate of
! convergence. Routine is based on the epsilon algorithm
! of P. Wynn. An estimate of the absolute error is also
! given.
double precision, intent(in) :: E0,E1,E2
double precision, intent(out) :: abserr, result
!locals
double precision, parameter :: ten = 10.0d0
double precision, parameter :: one = 1.0d0
double precision :: small, delta2, delta1
double precision :: tol2, tol1, err2, err1,ss
small = spacing(one)
delta2 = E2 - E1
delta1 = E1 - E0
err2 = abs(delta2)
err1 = abs(delta1)
tol2 = max(abs(E2),abs(E1)) * small
tol1 = max(abs(E1),abs(E0)) * small
if ( ( err1 <= tol1 ) .or. err2 <= tol2) then
C IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
C ACCURACY, CONVERGENCE IS ASSUMED.
result = E2
abserr = err1 + err2 + E2*small*ten
else
ss = one/delta2 - one/delta1
if (abs(ss*E1) <= 1.0d-3) then
result = E2
abserr = err1 + err2 + E2*small*ten
else
result = E1 + one/ss
abserr = err1 + err2 + abs(result-E2)
endif
endif
end subroutine dea3
subroutine dqagp(f,a,b,npts,points,epsabs,epsrel,limit,result1,
* abserr,neval,ier)
! use functionInterface
implicit none
integer, intent(in) :: npts,limit
double precision,dimension(npts), intent(in) :: points
double precision, intent(in) :: a, b, epsabs,epsrel
double precision, intent(out) :: result1,abserr
integer, intent(out) :: neval,ier
double precision :: f
!Locals
double precision,dimension(limit) :: alist, blist, rlist, elist
double precision,dimension(npts+2) :: pts
integer, dimension(limit) :: iord, level
integer, dimension(npts+2) :: ndin
integer ::last
external f
CALL dqagpe(f,a,b,npts,points,epsabs,epsrel,limit,result1,
* abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin
$ ,last)
end subroutine dqagp
subroutine dqagpe(f,a,b,npts,points,epsabs,epsrel,limit,result,
* abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin,
* last)
! use functionInterface
implicit none
c***begin prologue dqagpe
c***date written 800101 (yymmdd)
c***revision date 830518 (yymmdd)
c***category no. h2a2a1
c***keywords automatic integrator, general-purpose,
! singularities at user specified points,
! extrapolation, globally adaptive.
c***author piessens,robert ,appl. math. & progr. div. - k.u.leuven
! de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose the routine calculates an approximation result to a given
! definite integral i = integral of f over (a,b), hopefully
! satisfying following claim for accuracy abs(i-result).le.
! max(epsabs,epsrel*abs(i)). break points of the integration
! interval, where local difficulties of the integrand may
! occur(e.g. singularities,discontinuities),provided by user.
c***description
!
! computation of a definite integral
! standard fortran subroutine
! double precision version
!
! parameters
! on entry
! f - double precision
! function subprogram defining the integrand
! function f(x). the actual name for f needs to be
! declared e x t e r n a l in the driver program.
!
! a - double precision
! lower limit of integration
!
! b - double precision
! upper limit of integration
!
! npts2 - integer
! number equal to two more than the number of
! user-supplied break points within the integration
! range, npts2.ge.2.
! if npts2.lt.2, the routine will end with ier = 6.
!
! points - double precision
! vector of dimension npts2, the first (npts2-2)
! elements of which are the user provided break
! points. if these points do not constitute an
! ascending sequence there will be an automati!
! sorting.
!
! epsabs - double precision
! absolute accuracy requested
! epsrel - double precision
! relative accuracy requested
! if epsabs.le.0
! and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
! the routine will end with ier = 6.
!
! limit - integer
! gives an upper bound on the number of subintervals
! in the partition of (a,b), limit.ge.npts2
! if limit.lt.npts2, the routine will end with
! ier = 6.
!
! on return
! result - double precision
! approximation to the integral
!
! abserr - double precision
! estimate of the modulus of the absolute error,
! which should equal or exceed abs(i-result)
!
! neval - integer
! number of integrand evaluations
!
! ier - integer
! ier = 0 normal and reliable termination of the
! routine. it is assumed that the requested
! accuracy has been achieved.
! ier.gt.0 abnormal termination of the routine.
! the estimates for integral and error are
! less reliable. it is assumed that the
! requested accuracy has not been achieved.
! error messages
! ier = 1 maximum number of subdivisions allowed
! has been achieved. one can allow more
! subdivisions by increasing the value of
! limit (and taking the according dimension
! adjustments into account). however, if
! this yields no improvement it is advised
! to analyze the integrand in order to
! determine the integration difficulties. if
! the position of a local difficulty can be
! determined (i.e. singularity,
! discontinuity within the interval), it
! should be supplied to the routine as an
! element of the vector points. if necessary
! an appropriate special-purpose integrator
! must be used, which is designed for
! handling the type of difficulty involved.
! = 2 the occurrence of roundoff error is
! detected, which prevents the requested
! tolerance from being achieved.
! the error may be under-estimated.
! = 3 extremely bad integrand behaviour occurs
! at some points of the integration
! interval.
! = 4 the algorithm does not converge.
! roundoff error is detected in the
! extrapolation table. it is presumed that
! the requested tolerance cannot be
! achieved, and that the returned result is
! the best which can be obtained.
! = 5 the integral is probably divergent, or
! slowly convergent. it must be noted that
! divergence can occur with any other value
! of ier.gt.0.
! = 6 the input is invalid because
! npts2.lt.2 or
! break points are specified outside
! the integration range or
! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.npts2.
! result, abserr, neval, last, rlist(1),
! and elist(1) are set to zero. alist(1) and
! blist(1) are set to a and b respectively.
!
! alist - double precision
! vector of dimension at least limit, the first
! last elements of which are the left end points
! of the subintervals in the partition of the given
! integration range (a,b)
!
! blist - double precision
! vector of dimension at least limit, the first
! last elements of which are the right end points
! of the subintervals in the partition of the given
! integration range (a,b)
!
! rlist - double precision
! vector of dimension at least limit, the first
! last elements of which are the integral
! approximations on the subintervals
!
! elist - double precision
! vector of dimension at least limit, the first
! last elements of which are the moduli of the
! absolute error estimates on the subintervals
!
! pts - double precision
! vector of dimension at least npts2, containing the
! integration limits and the break points of the
! interval in ascending sequence.
!
! level - integer
! vector of dimension at least limit, containing the
! subdivision levels of the subinterval, i.e. if
! (aa,bb) is a subinterval of (p1,p2) where p1 as
! well as p2 is a user-provided break point or
! integration limit, then (aa,bb) has level l if
! abs(bb-aa) = abs(p2-p1)*2**(-l).
!
! ndin - integer
! vector of dimension at least npts2, after first
! integration over the intervals (pts(i)),pts(i+1),
! i = 0,1, ..., npts2-2, the error estimates over
! some of the intervals may have been increased
! artificially, in order to put their subdivision
! forward. if this happens for the subinterval
! numbered k, ndin(k) is put to 1, otherwise
! ndin(k) = 0.
!
! iord - integer
! vector of dimension at least limit, the first k
! elements of which are pointers to the
! error estimates over the subintervals,
! such that elist(iord(1)), ..., elist(iord(k))
! form a decreasing sequence, with k = last
! if last.le.(limit/2+2), and k = limit+1-last
! otherwise
!
! last - integer
! number of subintervals actually produced in the
! subdivisions process
!
c***references (none)
c***routines called d1mach,dqelg,dqk21,dqpsrt
c***end prologue dqagpe
integer, intent(in) :: npts,limit
double precision,dimension(npts), intent(in) :: points
double precision, intent(in) :: a, b, epsabs,epsrel
double precision, intent(out) :: result,abserr
integer, intent(out) :: neval,ier
double precision,dimension(limit), intent(out) :: alist, blist
double precision,dimension(limit), intent(out) :: rlist, elist
double precision,dimension(npts+2),intent(out) :: pts
integer, dimension(limit), intent(out) :: iord, level
integer, dimension(npts+2), intent(out) :: ndin
integer ::last
double precision :: f
! locals
double precision :: area,area1,area12,area2,a1,
* a2,b1,b2,correc,abseps,defabs,defab1,defab2,
* dres,epmach,erlarg,erlast,errbnd,
* errmax,error1,erro12,error2,errsum,ertest,oflow,
* resa,resabs,reseps,sign,temp,uflow, hSplit
double precision, dimension(3) :: res3la(3)
double precision, dimension(52) :: rlist2(52)
integer :: i,id,ierro,ind1,ind2,ip1,iroff1,iroff2,iroff3,j,
* jlow,jupbnd,k,ksgn,ktmin,levcur,levmax,maxerr,
* nint,nintp1,npts2,nres,nrmax,numrl2
logical :: extrap,noext
external f
!
!
!
!
! the dimension of rlist2 is determined by the value of
! limexp in subroutine epsalg (rlist2 should be of dimension
! (limexp+2) at least).
!
!
! list of major variables
! -----------------------
!
! alist - list of left end points of all subintervals
! considered up to now
! blist - list of right end points of all subintervals
! considered up to now
! rlist(i) - approximation to the integral over
! (alist(i),blist(i))
! rlist2 - array of dimension at least limexp+2
! containing the part of the epsilon table which
! is still needed for further computations
! elist(i) - error estimate applying to rlist(i)
! maxerr - pointer to the interval with largest error
! estimate
! errmax - elist(maxerr)
! erlast - error on the interval currently subdivided
! (before that subdivision has taken place)
! area - sum of the integrals over the subintervals
! errsum - sum of the errors over the subintervals
! errbnd - requested accuracy max(epsabs,epsrel*
! abs(result))
! *****1 - variable for the left subinterval
! *****2 - variable for the right subinterval
! last - index for subdivision
! nres - number of calls to the extrapolation routine
! numrl2 - number of elements in rlist2. if an appropriate
! approximation to the compounded integral has
! been obtained, it is put in rlist2(numrl2) after
! numrl2 has been increased by one.
! erlarg - sum of the errors over the intervals larger
! than the smallest interval considered up to now
! extrap - logical variable denoting that the routine
! is attempting to perform extrapolation. i.e.
! before subdividing the smallest interval we
! try to decrease the value of erlarg.
! noext - logical variable denoting that extrapolation is
! no longer allowed (true-value)
!
! machine dependent constants
! ---------------------------
!
! epmach is the largest relative spacing.
! uflow is the smallest positive magnitude.
! oflow is the largest positive magnitude.
!
c***first executable statement dqagpe
epmach = d1mach(4)
uflow = d1mach(1)
oflow = d1mach(2)
!
! test on validity of parameters
! -----------------------------
!
hSplit = 0.2D0
ier = 0
neval = 0
last = 0
result = 0.0d+00
abserr = 0.0d+00
alist(1) = a
blist(1) = b
rlist(1) = 0.0d+00
elist(1) = 0.0d+00
iord(1) = 0
level(1) = 0
npts2 = npts+2
if((npts2.lt.2).or.(limit.le.npts).or.
& ((epsabs.le.0.0d+00).and.
& (epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)))) then
ier = 6
go to 999
endif
sign = 1.0d+00
if(a.gt.b) then
go to 999
endif
if (npts>0) then
if(any(points(1:npts)<=a).or.any(b<=points(1:npts))) then
ier = 6
go to 999
endif
endif
!
! if any break points are provided, sort them into an
! ascending sequence.
!
pts(1) = a
pts(npts+2) = b
do i = 1,npts
pts(i+1) = minval(points(i:npts))
enddo
!
! compute first integral and error approximations.
! ------------------------------------------------
!
nint = npts+1;
a1 = pts(1);
resabs = 0.0d+00
do i = 1,nint
b1 = pts(i+1)
if (b1-a1 > hSplit) then
call dqk21(f,a1,b1,area1,error1,defabs,resa)
!call dqk15(f,a1,b1,area1,error1,defabs,resa)
else
call dqkl9(f,a1,b1,area1,error1,defabs,resa)
endif
abserr = abserr + error1
result = result + area1
ndin(i) = 0
if(error1.eq.resa.and.error1.ne.0.0d+00) ndin(i) = 1
resabs = resabs + defabs
level(i) = 0
elist(i) = error1
alist(i) = a1
blist(i) = b1
rlist(i) = area1
iord(i) = i
a1 = b1
enddo !50 continue
errsum = 0.0d+00
do i = 1,nint
if(ndin(i).eq.1) elist(i) = abserr
errsum = errsum+elist(i)
enddo !55 continue
!
! test on accuracy.
!
last = nint
neval = 21*nint
dres = dabs(result)
errbnd = dmax1(epsabs,epsrel*dres)
if(abserr.le.0.1d+03*epmach*resabs.and.abserr.gt.errbnd) ier = 2
if(nint.eq.1) go to 80
do 70 i = 1,npts
jlow = i+1
ind1 = iord(i)
do 60 j = jlow,nint
ind2 = iord(j)
if(elist(ind1).gt.elist(ind2)) go to 60
ind1 = ind2
k = j
60 continue
if(ind1.eq.iord(i)) go to 70
iord(k) = iord(i)
iord(i) = ind1
70 continue
if(limit.lt.npts2) ier = 1
80 if(ier.ne.0.or.abserr.le.errbnd) go to 210
!
! initialization
! --------------
!
rlist2(1) = result
maxerr = iord(1)
errmax = elist(maxerr)
area = result
nrmax = 1
nres = 0
numrl2 = 1
ktmin = 0
extrap = .false.
noext = .false.
erlarg = errsum
ertest = errbnd
levmax = 1
iroff1 = 0
iroff2 = 0
iroff3 = 0
ierro = 0
abserr = oflow
ksgn = -1
if(dres.ge.(0.1d+01-0.5d+02*epmach)*resabs) ksgn = 1
!
! main do-loop
! ------------
!
do 160 last = npts2,limit
!
! bisect the subinterval with the nrmax-th largest error
! estimate.
!
levcur = level(maxerr)+1
a1 = alist(maxerr)
b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
a2 = b1
b2 = blist(maxerr)
erlast = errmax
if (b1-a1 > hSplit) then
call dqk21(f,a1,b1,area1,error1,resa,defab1)
call dqk21(f,a2,b2,area2,error2,resa,defab2)
!call dqk15(f,a1,b1,area1,error1,resa,defab1)
!call dqk15(f,a2,b2,area2,error2,resa,defab2)
else
call dqkl9(f,a1,b1,area1,error1,resa,defab1)
call dqkl9(f,a2,b2,area2,error2,resa,defab2)
endif
!
! improve previous approximations to integral
! and error and test for accuracy.
!
neval = neval+42
area12 = area1+area2
erro12 = error1+error2
errsum = errsum+erro12-errmax
area = area+area12-rlist(maxerr)
if(defab1.eq.error1.or.defab2.eq.error2) go to 95
if(dabs(rlist(maxerr)-area12).gt.0.1d-04*dabs(area12)
* .or.erro12.lt.0.99d+00*errmax) go to 90
if(extrap) iroff2 = iroff2+1
if(.not.extrap) iroff1 = iroff1+1
90 if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1
95 level(maxerr) = levcur
level(last) = levcur
rlist(maxerr) = area1
rlist(last) = area2
errbnd = dmax1(epsabs,epsrel*dabs(area))
!
! test for roundoff error and eventually set error flag.
!
if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2
if(iroff2.ge.5) ierro = 3
!
! set error flag in the case that the number of
! subintervals equals limit.
!
if(last.eq.limit) ier = 1
!
! set error flag in the case of bad integrand behaviour
! at a point of the integration range
!
if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)*
* (dabs(a2)+0.1d+04*uflow)) ier = 4
!
! append the newly-created intervals to the list.
!
if(error2.gt.error1) go to 100
alist(last) = a2
blist(maxerr) = b1
blist(last) = b2
elist(maxerr) = error1
elist(last) = error2
go to 110
100 alist(maxerr) = a2
alist(last) = a1
blist(last) = b1
rlist(maxerr) = area2
rlist(last) = area1
elist(maxerr) = error2
elist(last) = error1
!
! call subroutine dqpsrt to maintain the descending ordering
! in the list of error estimates and select the subinterval
! with nrmax-th largest error estimate (to be bisected next).
!
110 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
! ***jump out of do-loop
if(errsum.le.errbnd) go to 190
! ***jump out of do-loop
if(ier.ne.0) go to 170
if(noext) go to 160
erlarg = erlarg-erlast
if(levcur+1.le.levmax) erlarg = erlarg+erro12
if(extrap) go to 120
!
! test whether the interval to be bisected next is the
! smallest interval.
!
if(level(maxerr)+1.le.levmax) go to 160
extrap = .true.
nrmax = 2
120 if(ierro.eq.3.or.erlarg.le.ertest) go to 140
!
! the smallest interval has the largest error.
! before bisecting decrease the sum of the errors over
! the larger intervals (erlarg) and perform extrapolation.
!
id = nrmax
jupbnd = last
if(last.gt.(2+limit/2)) jupbnd = limit+3-last
do 130 k = id,jupbnd
maxerr = iord(nrmax)
errmax = elist(maxerr)
! ***jump out of do-loop
if(level(maxerr)+1.le.levmax) go to 160
nrmax = nrmax+1
130 continue
!
! perform extrapolation.
!
140 numrl2 = numrl2+1
rlist2(numrl2) = area
if(numrl2.le.2) go to 155
call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres)
ktmin = ktmin+1
if(ktmin.gt.5.and.abserr.lt.0.1d-02*errsum) ier = 5
if(abseps.ge.abserr) go to 150
ktmin = 0
abserr = abseps
result = reseps
correc = erlarg
ertest = dmax1(epsabs,epsrel*dabs(reseps))
! ***jump out of do-loop
if(abserr.lt.ertest) go to 170
!
! prepare bisection of the smallest interval.
!
150 if(numrl2.eq.1) noext = .true.
if(ier.ge.5) go to 170
155 maxerr = iord(1)
errmax = elist(maxerr)
nrmax = 1
extrap = .false.
levmax = levmax + 1
erlarg = errsum
160 continue
!
! set the final result.
! ---------------------
!
!
170 if(abserr.eq.oflow) go to 190
if((ier+ierro).eq.0) go to 180
if(ierro.eq.3) abserr = abserr+correc
if(ier.eq.0) ier = 3
if(result.ne.0.0d+00.and.area.ne.0.0d+00)go to 175
if(abserr.gt.errsum)go to 190
if(area.eq.0.0d+00) go to 210
go to 180
175 if(abserr/dabs(result).gt.errsum/dabs(area))go to 190
!
! test on divergence.
!
180 if(ksgn.eq.(-1).and.dmax1(dabs(result),dabs(area)).le.
* resabs*0.1d-01) go to 210
if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03.or.
* errsum.gt.dabs(area)) ier = 6
go to 210
!
! compute global integral sum.
!
190 result = 0.0d+00
do 200 k = 1,last
result = result+rlist(k)
200 continue
abserr = errsum
210 if(ier.gt.2) ier = ier-1
result = result*sign
999 return
end subroutine dqagpe
subroutine dqk21(f,a,b,result,abserr,resabs,resasc)
! use functionInterface
implicit none
c***begin prologue dqk21
c***date written 800101 (yymmdd)
c***revision date 830518 (yymmdd)
c***category no. h2a1a2
c***keywords 21-point gauss-kronrod rules
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose to compute i = integral of f over (a,b), with error
c estimate
c j = integral of abs(f) over (a,b)
c***description
c
c integration rules
c standard fortran subroutine
c double precision version
c
c parameters
c on entry
c f - double precision
c function subprogram defining the integrand
c function f(x). the actual name for f needs to be
c declared e x t e r n a l in the driver program.
c
c a - double precision
c lower limit of integration
c
c b - double precision
c upper limit of integration
c
c on return
c result - double precision
c approximation to the integral i
c result is computed by applying the 21-point
c kronrod rule (resk) obtained by optimal addition
c of abscissae to the 10-point gauss rule (resg).
c
c abserr - double precision
c estimate of the modulus of the absolute error,
c which should not exceed abs(i-result)
c
c resabs - double precision
c approximation to the integral j
c
c resasc - double precision
c approximation to the integral of abs(f-i/(b-a))
c over (a,b)
c
c***references (none)
c***routines called d1mach
c***end prologue dqk21
c
double precision :: f, a,absc,abserr,b,centr,dhlgth,
* epmach,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
* resg,resk,reskh,result,uflow,wg,wgk,xgk
integer j,jtw,jtwm1
external f
c
dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
c
c the abscissae and weights are given for the interval (-1,1).
c because of symmetry only the positive abscissae and their
c corresponding weights are given.
c
c xgk - abscissae of the 21-point kronrod rule
c xgk(2), xgk(4), ... abscissae of the 10-point
c gauss rule
c xgk(1), xgk(3), ... abscissae which are optimally
c added to the 10-point gauss rule
c
c wgk - weights of the 21-point kronrod rule
c
c wg - weights of the 10-point gauss rule
c
c
c gauss quadrature weights and kronron quadrature abscissae and weights
c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
c bell labs, nov. 1981.
c
data wg ( 1) / 0.0666713443 0868813759 3568809893 332 d0 /
data wg ( 2) / 0.1494513491 5058059314 5776339657 697 d0 /
data wg ( 3) / 0.2190863625 1598204399 5534934228 163 d0 /
data wg ( 4) / 0.2692667193 0999635509 1226921569 469 d0 /
data wg ( 5) / 0.2955242247 1475287017 3892994651 338 d0 /
c
data xgk ( 1) / 0.9956571630 2580808073 5527280689 003 d0 /
data xgk ( 2) / 0.9739065285 1717172007 7964012084 452 d0 /
data xgk ( 3) / 0.9301574913 5570822600 1207180059 508 d0 /
data xgk ( 4) / 0.8650633666 8898451073 2096688423 493 d0 /
data xgk ( 5) / 0.7808177265 8641689706 3717578345 042 d0 /
data xgk ( 6) / 0.6794095682 9902440623 4327365114 874 d0 /
data xgk ( 7) / 0.5627571346 6860468333 9000099272 694 d0 /
data xgk ( 8) / 0.4333953941 2924719079 9265943165 784 d0 /
data xgk ( 9) / 0.2943928627 0146019813 1126603103 866 d0 /
data xgk ( 10) / 0.1488743389 8163121088 4826001129 720 d0 /
data xgk ( 11) / 0.0000000000 0000000000 0000000000 000 d0 /
c
data wgk ( 1) / 0.0116946388 6737187427 8064396062 192 d0 /
data wgk ( 2) / 0.0325581623 0796472747 8818972459 390 d0 /
data wgk ( 3) / 0.0547558965 7435199603 1381300244 580 d0 /
data wgk ( 4) / 0.0750396748 1091995276 7043140916 190 d0 /
data wgk ( 5) / 0.0931254545 8369760553 5065465083 366 d0 /
data wgk ( 6) / 0.1093871588 0229764189 9210590325 805 d0 /
data wgk ( 7) / 0.1234919762 6206585107 7958109831 074 d0 /
data wgk ( 8) / 0.1347092173 1147332592 8054001771 707 d0 /
data wgk ( 9) / 0.1427759385 7706008079 7094273138 717 d0 /
data wgk ( 10) / 0.1477391049 0133849137 4841515972 068 d0 /
data wgk ( 11) / 0.1494455540 0291690566 4936468389 821 d0 /
c
c
c list of major variables
c -----------------------
c
c centr - mid point of the interval
c hlgth - half-length of the interval
c absc - abscissa
c fval* - function value
c resg - result of the 10-point gauss formula
c resk - result of the 21-point kronrod formula
c reskh - approximation to the mean value of f over (a,b),
c i.e. to i/(b-a)
c
c
c machine dependent constants
c ---------------------------
c
c epmach is the largest relative spacing.
c uflow is the smallest positive magnitude.
c
c***first executable statement dqk21
epmach = d1mach(4)
uflow = d1mach(1)
c
centr = 0.5d+00*(a+b)
hlgth = 0.5d+00*(b-a)
dhlgth = dabs(hlgth)
c
c compute the 21-point kronrod approximation to
c the integral, and estimate the absolute error.
c
resg = 0.0d+00
fc = f(centr)
resk = wgk(11)*fc
resabs = dabs(resk)
do 10 j=1,5
jtw = 2*j
absc = hlgth*xgk(jtw)
fval1 = f(centr-absc)
fval2 = f(centr+absc)
fv1(jtw) = fval1
fv2(jtw) = fval2
fsum = fval1+fval2
resg = resg+wg(j)*fsum
resk = resk+wgk(jtw)*fsum
resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
10 continue
do 15 j = 1,5
jtwm1 = 2*j-1
absc = hlgth*xgk(jtwm1)
fval1 = f(centr-absc)
fval2 = f(centr+absc)
fv1(jtwm1) = fval1
fv2(jtwm1) = fval2
fsum = fval1+fval2
resk = resk+wgk(jtwm1)*fsum
resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
15 continue
reskh = resk*0.5d+00
resasc = wgk(11)*dabs(fc-reskh)
do 20 j=1,10
resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
20 continue
result = resk*hlgth
resabs = resabs*dhlgth
resasc = resasc*dhlgth
abserr = dabs((resk-resg)*hlgth)*10.0d0
if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00) then
abserr = resasc*dmin1(0.1d+01,
& (0.2d+03*abserr/resasc)**1.5d+00)
endif
if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
* ((epmach*0.5d+02)*resabs,abserr)
return
end subroutine dqk21
subroutine dqk15(f,a,b,result,abserr,resabs,resasc)
! use functionInterface
implicit none
c***begin prologue dqk15
c***date written 800101 (yymmdd)
c***revision date 830518 (yymmdd)
c***category no. h2a1a2
c***keywords 15-point gauss-kronrod rules
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose to compute i = integral of f over (a,b), with error
c estimate
c j = integral of abs(f) over (a,b)
c***description
c
c integration rules
c standard fortran subroutine
c double precision version
c
c parameters
c on entry
c f - double precision
c function subprogram defining the integrand
c function f(x). the actual name for f needs to be
c declared e x t e r n a l in the driver program.
c
c a - double precision
c lower limit of integration
c
c b - double precision
c upper limit of integration
c
c on return
c result - double precision
c approximation to the integral i
c result is computed by applying the 21-point
c kronrod rule (resk) obtained by optimal addition
c of abscissae to the 10-point gauss rule (resg).
c
c abserr - double precision
c estimate of the modulus of the absolute error,
c which should not exceed abs(i-result)
c
c resabs - double precision
c approximation to the integral j
c
c resasc - double precision
c approximation to the integral of abs(f-i/(b-a))
c over (a,b)
c
c***references (none)
c***routines called d1mach
c***end prologue dqk15
c
double precision :: f, a,absc,abserr,b,centr,dhlgth,
* epmach,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
* resg,resk,reskh,result,uflow,wg,wgk,xgk
integer j,jtw,jtwm1
external f
c
dimension fv1(7),fv2(7),wg(4),wgk(8),xgk(8)
c
c the abscissae and weights are given for the interval (-1,1).
c because of symmetry only the positive abscissae and their
c corresponding weights are given.
c
c xgk - abscissae of the 15-point kronrod rule
c xgk(2), xgk(4), ... abscissae of the 7-point
c gauss rule
c xgk(1), xgk(3), ... abscissae which are optimally
c added to the 7-point gauss rule
c
c wgk - weights of the 15-point kronrod rule
c
c wg - weights of the 7-point gauss rule
c
c
c gauss quadrature weights and kronron quadrature abscissae and weights
c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
c bell labs, nov. 1981.
c
data wg ( 1) / 0.129484966168869693270611432679082d0 /
data wg ( 2) / 0.279705391489276667901467771423780d0 /
data wg ( 3) / 0.381830050505118944950369775488975d0 /
data wg ( 4) / 0.417959183673469387755102040816327d0 /
data xgk ( 1) / 0.991455371120812639206854697526329d0 /
data xgk ( 2) / 0.949107912342758524526189684047851d0 /
data xgk ( 3) / 0.864864423359769072789712788640926d0 /
data xgk ( 4) / 0.741531185599394439863864773280788d0 /
data xgk ( 5) / 0.586087235467691130294144838258730d0 /
data xgk ( 6) / 0.405845151377397166906606412076961d0 /
data xgk ( 7) / 0.207784955007898467600689403773245d0 /
data xgk ( 8) / 0.000000000000000000000000000000000d0 /
data wgk ( 1) / 0.022935322010529224963732008058970d0/
data wgk ( 2) / 0.063092092629978553290700663189204d0 /
data wgk ( 3) / 0.104790010322250183839876322541518d0 /
data wgk ( 4) / 0.140653259715525918745189590510238d0 /
data wgk ( 5) / 0.169004726639267902826583426598550d0 /
data wgk ( 6) / 0.190350578064785409913256402421014d0 /
data wgk ( 7) / 0.204432940075298892414161999234649d0 /
data wgk ( 8) / 0.209482141084727828012999174891714d0 /
c
c
c list of major variables
c -----------------------
c
c centr - mid point of the interval
c hlgth - half-length of the interval
c absc - abscissa
c fval* - function value
c resg - result of the 10-point gauss formula
c resk - result of the 21-point kronrod formula
c reskh - approximation to the mean value of f over (a,b),
c i.e. to i/(b-a)
c
c
c machine dependent constants
c ---------------------------
c
c epmach is the largest relative spacing.
c uflow is the smallest positive magnitude.
c
c***first executable statement dqk15
epmach = d1mach(4)
uflow = d1mach(1)
c
centr = 0.5d+00*(a+b)
hlgth = 0.5d+00*(b-a)
dhlgth = dabs(hlgth)
c
c compute the 15-point kronrod approximation to
c the integral, and estimate the absolute error.
c
fc = f(centr)
resk = wgk(8)*fc
resg = wg(4)*fc
resabs = dabs(resk)
do 10 j=1,3
jtw = 2*j
absc = hlgth*xgk(jtw)
fval1 = f(centr-absc)
fval2 = f(centr+absc)
fv1(jtw) = fval1
fv2(jtw) = fval2
fsum = fval1+fval2
resg = resg+wg(j)*fsum
resk = resk+wgk(jtw)*fsum
resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
10 continue
do 15 j = 1,4
jtwm1 = 2*j-1
absc = hlgth*xgk(jtwm1)
fval1 = f(centr-absc)
fval2 = f(centr+absc)
fv1(jtwm1) = fval1
fv2(jtwm1) = fval2
fsum = fval1+fval2
resk = resk+wgk(jtwm1)*fsum
resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
15 continue
reskh = resk*0.5d+00
resasc = wgk(8)*dabs(fc-reskh)
do 20 j=1,7
resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
20 continue
result = resk*hlgth
resabs = resabs*dhlgth
resasc = resasc*dhlgth
abserr = dabs((resk-resg)*hlgth)*10.0D0
if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00) then
abserr = resasc*dmin1(0.1d+01,
& (0.2d+03*abserr/resasc)**1.5d+00)
endif
if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
* ((epmach*0.5d+02)*resabs,abserr)
return
end subroutine dqk15
subroutine dqk9(f,a,b,result,abserr,resabs,resasc)
! use functionInterface
implicit none
c***begin prologue dqk15
c***date written 800101 (yymmdd)
c***revision date 830518 (yymmdd)
c***category no. h2a1a2
c***keywords 15-point gauss-kronrod rules extended from a 3 point gaus rule
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose to compute i = integral of f over (a,b), with error
c estimate
c j = integral of abs(f) over (a,b)
c***description
c
c integration rules
c standard fortran subroutine
c double precision version
c
c parameters
c on entry
c f - double precision
c function subprogram defining the integrand
c function f(x). the actual name for f needs to be
c declared e x t e r n a l in the driver program.
c
c a - double precision
c lower limit of integration
c
c b - double precision
c upper limit of integration
c
c on return
c result - double precision
c approximation to the integral i
c result is computed by applying the 21-point
c kronrod rule (resk) obtained by optimal addition
c of abscissae to the 10-point gauss rule (resg).
c
c abserr - double precision
c estimate of the modulus of the absolute error,
c which should not exceed abs(i-result)
c
c resabs - double precision
c approximation to the integral j
c
c resasc - double precision
c approximation to the integral of abs(f-i/(b-a))
c over (a,b)
c
c***references (none)
c***routines called d1mach
c***end prologue dqk15
c
double precision :: f, a,absc,abserr,b,centr,dhlgth,
* epmach,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
* resg,resk0,resk,reskh,result,uflow,wg,wgk0,wgk,xgk
integer j,jtw,jtwm1
external f
c
dimension fv1(7),fv2(7),wg(2),wgk0(4),wgk(8),xgk(8)
c
c the abscissae and weights are given for the interval (-1,1).
c because of symmetry only the positive abscissae and their
c corresponding weights are given.
c
c xgk - abscissae of the 15-point kronrod rule
! xgk(4), xgk(8) abscissae of the 3-point gauss rule
c xgk(2), xgk(4),xgk(6), xgk(8) ... abscissae of the 7-point
c kronrod rule
c xgk(1), xgk(3), ... abscissae which are optimally
c added to the 7-point kronrod rule
c
c wgk - weights of the 15-point kronrod rule
!
! wgk0 - weights of the 7-point kronrod rule
c
c wg - weights of the 3-point gauss rule
c
c
c gauss quadrature weights and kronrod quadrature abscissae and weights
c as evaluated in quadruple precision by Patterson
c
data wg ( 1) / 0.5555555555555555D+00/
data wg ( 2) / 0.8888888888888889D+00/
data wgk0 ( 1) / 0.1046562260264673D+00/
data wgk0 ( 2) / 0.2684880898683335D+00/
data wgk0 ( 3) / 0.4013974147759622D+00/
data wgk0 ( 4) / 0.4509165386584741D+00/
data xgk ( 1) / 0.9938319632127550D+00/
data xgk ( 2) / 0.9604912687080203D+00/
data xgk ( 3) / 0.8884592328722570D+00 /
data xgk ( 4) / 0.7745966692414834D+00/
data xgk ( 5) / 0.6211029467372264D+00/
data xgk ( 6) / 0.4342437493468026D+00/
data xgk ( 7) / 0.2233866864289669D+00 /
data xgk ( 8) / 0.000000000000000000000000000000000d0 /
data wgk ( 1) / 0.1700171962994028D-01/
data wgk ( 2) / 0.5160328299707982D-01/
data wgk ( 3) / 0.9292719531512452D-01/
data wgk ( 4) / 0.1344152552437843D+00/
data wgk ( 5) / 0.1715119091363914D+00/
data wgk ( 6) / 0.2006285293769890D+00/
data wgk ( 7) / 0.2191568584015875D+00/
data wgk ( 8) / 0.2255104997982067D+00/
c
c
c list of major variables
c -----------------------
c
c centr - mid point of the interval
c hlgth - half-length of the interval
c absc - abscissa
c fval* - function value
c resg - result of the 10-point gauss formula
c resk - result of the 21-point kronrod formula
c reskh - approximation to the mean value of f over (a,b),
c i.e. to i/(b-a)
c
c
c machine dependent constants
c ---------------------------
c
c epmach is the largest relative spacing.
c uflow is the smallest positive magnitude.
c
c***first executable statement dqk15
epmach = d1mach(4)
uflow = d1mach(1)
c
centr = 0.5d+00*(a+b)
hlgth = 0.5d+00*(b-a)
dhlgth = dabs(hlgth)
c
c compute the 15-point kronrod approximation to
c the integral, and estimate the absolute error.
c
fc = f(centr)
resk = wgk(8)*fc
resk0 = wgk0(4)*fc
resabs = dabs(resk)
do 10 j=1,3
jtw = 2*j
absc = hlgth * xgk(jtw)
fval1 = f(centr-absc)
fval2 = f(centr+absc)
fv1(jtw) = fval1
fv2(jtw) = fval2
fsum = fval1 + fval2
resk0 = resk0 + wgk0(j) * fsum
resk = resk+wgk(jtw)*fsum
resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
10 continue
resg = wg(2)*fc + wg(1)*(fv1(4) + fv2(4))
do 15 j = 1,4
jtwm1 = 2*j-1
absc = hlgth * xgk(jtwm1)
fval1 = f( centr - absc )
fval2 = f( centr + absc )
fv1(jtwm1) = fval1
fv2(jtwm1) = fval2
fsum = fval1 + fval2
resk = resk + wgk(jtwm1) * fsum
resabs = resabs + wgk(jtwm1) * (dabs(fval1) + dabs(fval2))
15 continue
reskh = resk*0.5d+00
resasc = wgk(8)*dabs(fc-reskh)
do 20 j=1,7
resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
20 continue
resg = resg * hlgth
resk0 = resk0 * hlgth
resk = resk * hlgth
resabs = resabs * dhlgth
resasc = resasc * dhlgth
result = resk
call dea3(resg,resk0,resk,abserr,result)
abserr = max((dabs(resk-resk0) + dabs(resg-resk0))
& * 10.0D0, abserr)
if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00) then
abserr = resasc * dmin1(0.1d+01,
& (0.2d+03*abserr/resasc)**1.5d+00)
endif
if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
* ((epmach*0.5d+02)*resabs,abserr)
return
end subroutine dqk9
subroutine dqkl9(f,a,b,result,abserr,resabs,resasc)
! use functionInterface
implicit none
c***begin prologue dqk15
c***date written 800101 (yymmdd)
c***revision date 830518 (yymmdd)
c***category no. h2a1a2
c***keywords 15-point gauss-kronrod rules extended from a 3 point gaus rule
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose to compute i = integral of f over (a,b), with error
c estimate
c j = integral of abs(f) over (a,b)
c***description
c
c integration rules
c standard fortran subroutine
c double precision version
c
c parameters
c on entry
c f - double precision
c function subprogram defining the integrand
c function f(x). the actual name for f needs to be
c declared e x t e r n a l in the driver program.
c
c a - double precision
c lower limit of integration
c
c b - double precision
c upper limit of integration
c
c on return
c result - double precision
c approximation to the integral i
c result is computed by applying the 21-point
c kronrod rule (resk) obtained by optimal addition
c of abscissae to the 10-point gauss rule (resg).
c
c abserr - double precision
c estimate of the modulus of the absolute error,
c which should not exceed abs(i-result)
c
c resabs - double precision
c approximation to the integral j
c
c resasc - double precision
c approximation to the integral of abs(f-i/(b-a))
c over (a,b)
c
c***references (none)
c***routines called d1mach
c***end prologue dqk15
c
double precision :: f, a,absc,abserr,b,centr,dhlgth,
* epmach,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
* resg,resk0,resk,reskh,result,uflow,wg,wgk0,wgk,xgk
integer j,jtw,jtwm1
external f
c
dimension fv1(7),fv2(7),wg(2),wgk0(3),wgk(5),xgk(5)
c
c the abscissae and weights are given for the interval (-1,1).
c because of symmetry only the positive abscissae and their
c corresponding weights are given.
c
c xgk - abscissae of the 9-point Gauss-kronrod-lobatto rule
! xgk(1), xgk(5) abscissae of the 3-point gauss-lobatto rule
c xgk(1), xgk(3),xgk(5) abscissae of the 5-point
c kronrod rule
c xgk(2), xgk(4), ... abscissae which are optimally
c added to the 5-point kronrod rule
c
c wgk - weights of the 9-point kronrod rule
!
! wgk0 - weights of the 5-point kronrod rule
c
c wg - weights of the 3-point gauss rule
c
c
c gauss quadrature weights and kronrod quadrature abscissae and weights
c as evaluated in quadruple precision by Patterson
c
data wg ( 1) / 0.33333333333333333333333333333333333D+00/
data wg ( 2) / 0.13333333333333333333333333333333333D+01/
data wgk0 ( 1) / 0.1000000000000000D+00/
data wgk0 ( 2) / 0.5444444444444445D+00/
data wgk0 ( 3) / 0.7111111111111111D+00/
data xgk ( 1) / 0.1000000000000000D+01/
data xgk ( 2) / 0.8904055275126688D+00/
data xgk ( 3) / 0.6546536707079772D+00/
data xgk ( 4) / 0.3409822659109930D+00/
data xgk ( 5) / 0.000000000000000000000000000000000d0 /
data wgk ( 1) / 0.3064373897707232D-01/
data wgk ( 2) / 0.1792626995532074D+00/
data wgk ( 3) / 0.2839787780481211D+00/
data wgk ( 4) / 0.3342337398164177D+00/
data wgk ( 5) / 0.3437620872103631D+00/
c
c
c list of major variables
c -----------------------
c
c centr - mid point of the interval
c hlgth - half-length of the interval
c absc - abscissa
c fval* - function value
c resg - result of the 10-point gauss formula
c resk - result of the 21-point kronrod formula
c reskh - approximation to the mean value of f over (a,b),
c i.e. to i/(b-a)
c
c
c machine dependent constants
c ---------------------------
c
c epmach is the largest relative spacing.
c uflow is the smallest positive magnitude.
c
c***first executable statement dqk15
epmach = d1mach(4)
uflow = d1mach(1)
c
centr = 0.5d+00*(a+b)
hlgth = 0.5d+00*(b-a)
dhlgth = dabs(hlgth)
c
c compute the 15-point kronrod approximation to
c the integral, and estimate the absolute error.
c
fc = f(centr)
resk = wgk(5)*fc
resk0 = wgk0(3)*fc
resabs = dabs(resk)
do 10 j=1,2
jtw = 2*j - 1
absc = hlgth * xgk(jtw)
fval1 = f(centr-absc)
fval2 = f(centr+absc)
fv1(jtw) = fval1
fv2(jtw) = fval2
fsum = fval1 + fval2
resk0 = resk0 + wgk0(j) * fsum
resk = resk+wgk(jtw)*fsum
resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
10 continue
resg = wg(2)*fc + wg(1)*(fv1(1) + fv2(1))
do 15 j = 1,2
jtwm1 = 2*j
absc = hlgth * xgk(jtwm1)
fval1 = f( centr - absc )
fval2 = f( centr + absc )
fv1(jtwm1) = fval1
fv2(jtwm1) = fval2
fsum = fval1 + fval2
resk = resk + wgk(jtwm1) * fsum
resabs = resabs + wgk(jtwm1) * (dabs(fval1) + dabs(fval2))
15 continue
reskh = resk*0.5d+00
resasc = wgk(5)*dabs(fc-reskh)
do 20 j=1,4
resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
20 continue
resg = resg * hlgth
resk0 = resk0 * hlgth
resk = resk * hlgth
resabs = resabs * dhlgth
resasc = resasc * dhlgth
result = resk
call dea3(resg,resk0,resk,abserr,result)
abserr = max((dabs(resk-resk0) + dabs(resg-resk0))* 10.0D0,abserr)
if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00) then
abserr = resasc * dmin1(0.1d+01,
& (0.2d+03*abserr/resasc)**1.5d+00)
endif
if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
* ((epmach*0.5d+02)*resabs,abserr)
return
end subroutine dqkl9
subroutine dqpsrt(limit,last,maxerr,ermax,elist,iord,nrmax)
implicit none
c***begin prologue dqpsrt
c***refer to dqage,dqagie,dqagpe,dqawse
c***routines called (none)
c***revision date 810101 (yymmdd)
c***keywords sequential sorting
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose this routine maintains the descending ordering in the
c list of the local error estimated resulting from the
c interval subdivision process. at each call two error
c estimates are inserted using the sequential search
c method, top-down for the largest error estimate and
c bottom-up for the smallest error estimate.
c***description
c
c ordering routine
c standard fortran subroutine
c double precision version
c
c parameters (meaning at output)
c limit - integer
c maximum number of error estimates the list
c can contain
c
c last - integer
c number of error estimates currently in the list
c
c maxerr - integer
c maxerr points to the nrmax-th largest error
c estimate currently in the list
c
c ermax - double precision
c nrmax-th largest error estimate
c ermax = elist(maxerr)
c
c elist - double precision
c vector of dimension last containing
c the error estimates
c
c iord - integer
c vector of dimension last, the first k elements
c of which contain pointers to the error
c estimates, such that
c elist(iord(1)),..., elist(iord(k))
c form a decreasing sequence, with
c k = last if last.le.(limit/2+2), and
c k = limit+1-last otherwise
c
c nrmax - integer
c maxerr = iord(nrmax)
c
c***end prologue dqpsrt
c
double precision elist,ermax,errmax,errmin
integer i,ibeg,ido,iord,isucc,j,jbnd,jupbn,k,last,limit,maxerr,
* nrmax
dimension elist(last),iord(last)
c
c check whether the list contains more than
c two error estimates.
c
c***first executable statement dqpsrt
if(last.gt.2) go to 10
iord(1) = 1
iord(2) = 2
go to 90
c
c this part of the routine is only executed if, due to a
c difficult integrand, subdivision increased the error
c estimate. in the normal case the insert procedure should
c start after the nrmax-th largest error estimate.
c
10 errmax = elist(maxerr)
if(nrmax.eq.1) go to 30
ido = nrmax-1
do 20 i = 1,ido
isucc = iord(nrmax-1)
c ***jump out of do-loop
if(errmax.le.elist(isucc)) go to 30
iord(nrmax) = isucc
nrmax = nrmax-1
20 continue
c
c compute the number of elements in the list to be maintained
c in descending order. this number depends on the number of
c subdivisions still allowed.
c
30 jupbn = last
if(last.gt.(limit/2+2)) jupbn = limit+3-last
errmin = elist(last)
c
c insert errmax by traversing the list top-down,
c starting comparison from the element elist(iord(nrmax+1)).
c
jbnd = jupbn-1
ibeg = nrmax+1
if(ibeg.gt.jbnd) go to 50
do 40 i=ibeg,jbnd
isucc = iord(i)
c ***jump out of do-loop
if(errmax.ge.elist(isucc)) go to 60
iord(i-1) = isucc
40 continue
50 iord(jbnd) = maxerr
iord(jupbn) = last
go to 90
c
c insert errmin by traversing the list bottom-up.
c
60 iord(i-1) = maxerr
k = jbnd
do 70 j=i,jbnd
isucc = iord(k)
c ***jump out of do-loop
if(errmin.lt.elist(isucc)) go to 80
iord(k+1) = isucc
k = k-1
70 continue
iord(i) = last
go to 90
80 iord(k+1) = last
c
c set maxerr and ermax.
c
90 maxerr = iord(nrmax)
ermax = elist(maxerr)
return
end subroutine dqpsrt
subroutine dqelg(n,epstab,result,abserr,res3la,nres)
implicit none
c***begin prologue dqelg
c***refer to dqagie,dqagoe,dqagpe,dqagse
c***routines called d1mach
c***revision date 830518 (yymmdd)
c***keywords epsilon algorithm, convergence acceleration,
c extrapolation
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
c de doncker,elise,appl. math & progr. div. - k.u.leuven
c***purpose the routine determines the limit of a given sequence of
c approximations, by means of the epsilon algorithm of
c p.wynn. an estimate of the absolute error is also given.
c the condensed epsilon table is computed. only those
c elements needed for the computation of the next diagonal
c are preserved.
c***description
c
c epsilon algorithm
c standard fortran subroutine
c double precision version
c
c parameters
c n - integer
c epstab(n) contains the new element in the
c first column of the epsilon table.
c
c epstab - double precision
c vector of dimension 52 containing the elements
c of the two lower diagonals of the triangular
c epsilon table. the elements are numbered
c starting at the right-hand corner of the
c triangle.
c
c result - double precision
c resulting approximation to the integral
c
c abserr - double precision
c estimate of the absolute error computed from
c result and the 3 previous results
c
c res3la - double precision
c vector of dimension 3 containing the last 3
c results
c
c nres - integer
c number of calls to the routine
c (should be zero at first call)
c
c***end prologue dqelg
c
double precision abserr,dabs,delta1,delta2,delta3,dmax1,
* epmach,epsinf,epstab,error,err1,err2,err3,e0,e1,e1abs,e2,e3,
* oflow,res,result,res3la,ss,tol1,tol2,tol3
integer i,ib,ib2,ie,indx,k1,k2,k3,limexp,n,newelm,nres,num
dimension epstab(52),res3la(3)
c
c list of major variables
c -----------------------
c
c e0 - the 4 elements on which the computation of a new
c e1 element in the epsilon table is based
c e2
c e3 e0
c e3 e1 new
c e2
c newelm - number of elements to be computed in the new
c diagonal
c error - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
c result - the element in the new diagonal with least value
c of error
c
c machine dependent constants
c ---------------------------
c
c epmach is the largest relative spacing.
c oflow is the largest positive magnitude.
c limexp is the maximum number of elements the epsilon
c table can contain. if this number is reached, the upper
c diagonal of the epsilon table is deleted.
c
c***first executable statement dqelg
epmach = d1mach(4)
oflow = d1mach(2)
nres = nres+1
abserr = oflow
result = epstab(n)
if(n.lt.3) go to 100
limexp = 50
epstab(n+2) = epstab(n)
newelm = (n-1)/2
epstab(n) = oflow
num = n
k1 = n
do 40 i = 1,newelm
k2 = k1-1
k3 = k1-2
res = epstab(k1+2)
e0 = epstab(k3)
e1 = epstab(k2)
e2 = res
e1abs = dabs(e1)
delta2 = e2-e1
err2 = dabs(delta2)
tol2 = dmax1(dabs(e2),e1abs)*epmach
delta3 = e1-e0
err3 = dabs(delta3)
tol3 = dmax1(e1abs,dabs(e0))*epmach
if(err2.gt.tol2.or.err3.gt.tol3) go to 10
c
c if e0, e1 and e2 are equal to within machine
c accuracy, convergence is assumed.
c result = e2
c abserr = abs(e1-e0)+abs(e2-e1)
c
result = res
abserr = err2+err3
c ***jump out of do-loop
go to 100
10 e3 = epstab(k1)
epstab(k1) = e1
delta1 = e1-e3
err1 = dabs(delta1)
tol1 = dmax1(e1abs,dabs(e3))*epmach
c
c if two elements are very close to each other, omit
c a part of the table by adjusting the value of n
c
if(err1.le.tol1.or.err2.le.tol2.or.err3.le.tol3) go to 20
ss = 0.1d+01/delta1+0.1d+01/delta2-0.1d+01/delta3
epsinf = dabs(ss*e1)
c
c test to detect irregular behaviour in the table, and
c eventually omit a part of the table adjusting the value
c of n.
c
if(epsinf.gt.0.1d-03) go to 30
20 n = i+i-1
c ***jump out of do-loop
go to 50
c
c compute a new element and eventually adjust
c the value of result.
c
30 res = e1+0.1d+01/ss
epstab(k1) = res
k1 = k1-2
error = err2+dabs(res-e2)+err3
if(error.gt.abserr) go to 40
abserr = error
result = res
40 continue
c
c shift the table.
c
50 if(n.eq.limexp) n = 2*(limexp/2)-1
ib = 1
if((num/2)*2.eq.num) ib = 2
ie = newelm+1
do 60 i=1,ie
ib2 = ib+2
epstab(ib) = epstab(ib2)
ib = ib2
60 continue
if(num.eq.n) go to 80
indx = num-n+1
do 70 i = 1,n
epstab(i)= epstab(indx)
indx = indx+1
70 continue
80 if(nres.ge.4) go to 90
res3la(nres) = result
abserr = oflow
go to 100
c
c compute error estimate
c
90 abserr = dabs(result-res3la(3))+dabs(result-res3la(2))
* +dabs(result-res3la(1))
res3la(1) = res3la(2)
res3la(2) = res3la(3)
res3la(3) = result
100 abserr = dmax1(abserr,0.5d+01*epmach*dabs(result))
return
end subroutine dqelg
DOUBLE PRECISION FUNCTION D1MACH(I)
implicit none
C
C Double-precision machine constants.
C
C D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude.
C D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
C D1MACH( 3) = B**(-T), the smallest relative spacing.
C D1MACH( 4) = B**(1-T), the largest relative spacing.
C D1MACH( 5) = LOG10(B)
C
C Two more added much later:
C
C D1MACH( 6) = Infinity.
C D1MACH( 7) = Not-a-Number.
C
C Reference: Fox P.A., Hall A.D., Schryer N.L.,"Framework for a
C Portable Library", ACM Transactions on Mathematical
C Software, Vol. 4, no. 2, June 1978, PP. 177-188.
C
INTEGER , INTENT(IN) :: I
DOUBLE PRECISION, SAVE :: DMACH(7)
DOUBLE PRECISION :: B, EPS
DOUBLE PRECISION :: ONE = 1.0D0
DOUBLE PRECISION :: ZERO = 0.0D0
INTEGER :: EMAX,EMIN,T
DATA DMACH /7*0.0D0/
! First time through, get values from F90 INTRINSICS:
IF (DMACH(1) .EQ. 0.0D0) THEN
T = DIGITS(ONE)
B = DBLE(RADIX(ONE)) ! base number
EPS = SPACING(ONE)
EMIN = MINEXPONENT(ONE)
EMAX = MAXEXPONENT(ONE)
DMACH(1) = B**(EMIN-1) !TINY(ONE)
DMACH(2) = (B**(EMAX-1)) * (B-B*EPS) !HUGE(ONE)
DMACH(3) = EPS/B ! EPS/B
DMACH(4) = EPS
DMACH(5) = LOG10(B)
DMACH(6) = B**(EMAX+5) !infinity
DMACH(7) = ZERO/ZERO !nan
ENDIF
C
D1MACH = DMACH(I)
RETURN
END FUNCTION D1MACH
end module AdaptiveGaussKronrod
module Integration1DModule
implicit none
interface AdaptiveSimpson
module procedure AdaptiveSimpson2, AdaptiveSimpsonWithBreaks
end interface
! interface AdaptiveSimpson1
! module procedure AdaptiveSimpson1
! end interface
interface AdaptiveTrapz
module procedure AdaptiveTrapz1, AdaptiveTrapzWithBreaks
end interface
interface Romberg
module procedure Romberg1, RombergWithBreaks
end interface
INTERFACE DEA
MODULE PROCEDURE DEA
END INTERFACE
INTERFACE d1mach
MODULE PROCEDURE d1mach
END INTERFACE
contains
DOUBLE PRECISION FUNCTION D1MACH(I)
implicit none
C
C Double-precision machine constants.
C
C D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude.
C D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
C D1MACH( 3) = B**(-T), the smallest relative spacing.
C D1MACH( 4) = B**(1-T), the largest relative spacing.
C D1MACH( 5) = LOG10(B)
C
C Two more added much later:
C
C D1MACH( 6) = Infinity.
C D1MACH( 7) = Not-a-Number.
C
C Reference: Fox P.A., Hall A.D., Schryer N.L.,"Framework for a
C Portable Library", ACM Transactions on Mathematical
C Software, Vol. 4, no. 2, June 1978, PP. 177-188.
C
INTEGER , INTENT(IN) :: I
DOUBLE PRECISION, SAVE :: DMACH(7)
DOUBLE PRECISION :: B, EPS
DOUBLE PRECISION :: ONE = 1.0D0
DOUBLE PRECISION :: ZERO = 0.0D0
INTEGER :: EMAX,EMIN,T
DATA DMACH /7*0.0D0/
! First time through, get values from F90 INTRINSICS:
IF (DMACH(1) .EQ. 0.0D0) THEN
T = DIGITS(ONE)
B = DBLE(RADIX(ONE)) ! base number
EPS = SPACING(ONE)
EMIN = MINEXPONENT(ONE)
EMAX = MAXEXPONENT(ONE)
DMACH(1) = B**(EMIN-1) !TINY(ONE)
DMACH(2) = (B**(EMAX-1)) * (B-B*EPS) !HUGE(ONE)
DMACH(3) = EPS/B ! EPS/B
DMACH(4) = EPS
DMACH(5) = LOG10(B)
DMACH(6) = B**(EMAX+5) !infinity
DMACH(7) = ZERO/ZERO !nan
ENDIF
C
D1MACH = DMACH(I)
RETURN
END FUNCTION D1MACH
subroutine dea3(E0,E1,E2,abserr,result)
!***PURPOSE Given a slowly convergent sequence, this routine attempts
! to extrapolate nonlinearly to a better estimate of the
! sequence's limiting value, thus improving the rate of
! convergence. Routine is based on the epsilon algorithm
! of P. Wynn. An estimate of the absolute error is also
! given.
double precision, intent(in) :: E0,E1,E2
double precision, intent(out) :: abserr, result
!locals
double precision, parameter :: ten = 10.0d0
double precision, parameter :: one = 1.0d0
double precision :: small, delta2, delta1
double precision :: tol2, tol1, err2, err1,ss
small = spacing(one)
delta2 = E2 - E1
delta1 = E1 - E0
err2 = abs(delta2)
err1 = abs(delta1)
tol2 = max(abs(E2),abs(E1)) * small
tol1 = max(abs(E1),abs(E0)) * small
if ( ( err1 <= tol1 ) .or. err2 <= tol2) then
C IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
C ACCURACY, CONVERGENCE IS ASSUMED.
result = E2
abserr = err1 + err2 + E2*small*ten
else
ss = one/delta2 - one/delta1
if (abs(ss*E1) <= 1.0d-3) then
result = E2
abserr = err1 + err2 + E2*small*ten
else
result = E1 + one/ss
abserr = err1 + err2 + abs(result-E2)
endif
endif
end subroutine dea3
SUBROUTINE DEA(NEWFLG,SVALUE,LIMEXP,RESULT,ABSERR,EPSTAB,IERR)
C***BEGIN PROLOGUE DEA
C***DATE WRITTEN 800101 (YYMMDD)
C***REVISION DATE 871208 (YYMMDD)
C***CATEGORY NO. E5
C***KEYWORDS CONVERGENCE ACCELERATION,EPSILON ALGORITHM,EXTRAPOLATION
C***AUTHOR PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. -
C K. U. LEUVEN
C DE DONCKER-KAPENGA, ELISE,WESTERN MICHIGAN UNIVERSITY
C KAHANER, DAVID K., NATIONAL BUREAU OF STANDARDS
C STARKENBURG, C. B., NATIONAL BUREAU OF STANDARDS
C***PURPOSE Given a slowly convergent sequence, this routine attempts
C to extrapolate nonlinearly to a better estimate of the
C sequence's limiting value, thus improving the rate of
C convergence. Routine is based on the epsilon algorithm
C of P. Wynn. An estimate of the absolute error is also
C given.
C***DESCRIPTION
C
C Epsilon algorithm. Standard fortran subroutine.
C Double precision version.
C
C A R G U M E N T S I N T H E C A L L S E Q U E N C E
C
C NEWFLG - LOGICAL (INPUT and OUTPUT)
C On the first call to DEA set NEWFLG to .TRUE.
C (indicating a new sequence). DEA will set NEWFLG
C to .FALSE.
C
C SVALUE - DOUBLE PRECISION (INPUT)
C On the first call to DEA set SVALUE to the first
C term in the sequence. On subsequent calls set
C SVALUE to the subsequent sequence value.
C
C LIMEXP - INTEGER (INPUT)
C An integer equal to or greater than the total
C number of sequence terms to be evaluated. Do not
C change the value of LIMEXP until a new sequence
C is evaluated (NEWFLG=.TRUE.). LIMEXP .GE. 3
C
C RESULT - DOUBLE PRECISION (OUTPUT)
C Best approximation to the sequence's limit.
C
C ABSERR - DOUBLE PRECISION (OUTPUT)
C Estimate of the absolute error.
C
C EPSTAB - DOUBLE PRECISION (OUTPUT)
C Workvector of DIMENSION at least (LIMEXP+7).
C
C IERR - INTEGER (OUTPUT)
C IERR=0 Normal termination of the routine.
C IERR=1 The input is invalid because LIMEXP.LT.3.
C
C T Y P I C A L P R O B L E M S E T U P
C
C This sample problem uses the trapezoidal rule to evaluate the
C integral of the sin function from 0.0 to 0.5*PI (value = 1.0). The
C program implements the trapezoidal rule 8 times creating an
C increasingly accurate sequence of approximations to the integral.
C Each time the trapezoidal rule is used, it uses twice as many
C panels as the time before. DEA is called to obtain even more
C accurate estimates.
C
C PROGRAM SAMPLE
C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C DOUBLE PRECISION EPSTAB(57)
CC [57 = LIMEXP + 7]
C LOGICAL NEWFLG
C EXTERNAL F
C DATA LIMEXP/50/
C WRITE(*,*) ' NO. PANELS TRAP. APPROX'
C * ,' APPROX W/EA ABSERR'
C WRITE(*,*)
C HALFPI = DASIN(1.0D+00)
CC [UPPER INTEGRATION LIMIT = PI/2]
C NEWFLG = .TRUE.
CC [SET FLAG - 1ST DEA CALL]
C DO 10 I = 0,7
C NPARTS = 2 ** I
C WIDTH = HALFPI/NPARTS
C APPROX = 0.5D+00 * WIDTH * (F(0.0D+00) + F(HALFPI))
C DO 11 J = 1,NPARTS-1
C APPROX = APPROX + F(J * WIDTH) * WIDTH
C 11 CONTINUE
CC [END TRAPEZOIDAL RULE APPROX]
C SVALUE = APPROX
CC [SVALUE = NEW SEQUENCE VALUE]
C CALL DEA(NEWFLG,SVALUE,LIMEXP,RESULT,ABSERR,EPSTAB,IERR)
CC [CALL DEA FOR BETTER ESTIMATE]
C WRITE(*,12) NPARTS,APPROX,RESULT,ABSERR
C 12 FORMAT(' ',I4,T20,F16.13,T40,F16.13,T60,D11.4)
C 10 CONTINUE
C STOP
C END
C
C DOUBLE PRECISION FUNCTION F(X)
C DOUBLE PRECISION X
C F = DSIN(X)
CC [INTEGRAND]
C RETURN
C END
C
C Output from the above program will be:
C
C NO. PANELS TRAP. APPROX APPROX W/EA ABSERR
C
C 1 .7853981633974 .7853981633974 .7854D+00
C 2 .9480594489685 .9480594489685 .9760D+00
C 4 .9871158009728 .9994567212570 .2141D+00
C 8 .9967851718862 .9999667417647 .3060D-02
C 16 .9991966804851 .9999998781041 .6094D-03
C 32 .9997991943200 .9999999981026 .5767D-03
C 64 .9999498000921 .9999999999982 .3338D-04
C 128 .9999874501175 1.0000000000000 .1238D-06
C
C-----------------------------------------------------------------------
C***REFERENCES "Acceleration de la convergence en analyse numerique",
C C. Brezinski, "Lecture Notes in Math.", vol. 584,
C Springer-Verlag, New York, 1977.
C***ROUTINES CALLED D1MACH,XERROR
C***END PROLOGUE DEA
double precision, dimension(*), intent(inout) :: EPSTAB
double precision, intent(out) :: RESULT !, ABSERR
double precision, intent(inout) :: ABSERR
double precision, intent(in) :: SVALUE
INTEGER, INTENT(IN) :: LIMEXP
INTEGER, INTENT(OUT) :: IERR
LOGICAL, intent(INOUT) :: NEWFLG
DOUBLE PRECISION :: DELTA1,DELTA2,DELTA3,DRELPR,DEPRN,
1 ERROR,ERR1,ERR2,ERR3,E0,E1,E2,E3,RES,
2 SS,TOL1,TOL2,TOL3
double precision, dimension(3) :: RES3LA
INTEGER I,IB,IB2,IE,IN,K1,K2,K3,N,NEWELM,NUM,NRES
C
C
C LIMEXP is the maximum number of elements the
C epsilon table data can contain. The epsilon table
C is stored in the first (LIMEXP+2) entries of EPSTAB.
C
C
C LIST OF MAJOR VARIABLES
C -----------------------
C E0,E1,E2,E3 - DOUBLE PRECISION
C The 4 elements on which the computation of
C a new element in the epsilon table is based.
C NRES - INTEGER
C Number of extrapolation results actually
C generated by the epsilon algorithm in prior
C calls to the routine.
C NEWELM - INTEGER
C Number of elements to be computed in the
C new diagonal of the epsilon table. The
C condensed epsilon table is computed. Only
C those elements needed for the computation of
C the next diagonal are preserved.
C RES - DOUBLE PRECISION
C New element in the new diagonal of the
C epsilon table.
C ERROR - DOUBLE PRECISION
C An estimate of the absolute error of RES.
C Routine decides whether RESULT=RES or
C RESULT=SVALUE by comparing ERROR with
C ABSERR from the previous call.
C RES3LA - DOUBLE PRECISION
C Vector of DIMENSION 3 containing at most
C the last 3 results.
C
C
C MACHINE DEPENDENT CONSTANTS
C ---------------------------
C DRELPR is the largest relative spacing.
C
C***FIRST EXECUTABLE STATEMENT DEA
IF(LIMEXP.LT.3) THEN
IERR = 1
! CALL XERROR('LIMEXP IS LESS THAN 3',21,1,1)
GO TO 110
ENDIF
IERR = 0
RES3LA(1)=EPSTAB(LIMEXP+5)
RES3LA(2)=EPSTAB(LIMEXP+6)
RES3LA(3)=EPSTAB(LIMEXP+7)
RESULT=SVALUE
IF(NEWFLG) THEN
N=1
NRES=0
NEWFLG=.FALSE.
EPSTAB(N)=SVALUE
ABSERR=ABS(RESULT)
GO TO 100
ELSE
N=INT(EPSTAB(LIMEXP+3))
NRES=INT(EPSTAB(LIMEXP+4))
IF(N.EQ.2) THEN
EPSTAB(N)=SVALUE
ABSERR=.6D+01*ABS(RESULT-EPSTAB(1))
GO TO 100
ENDIF
ENDIF
EPSTAB(N)=SVALUE
DRELPR=D1MACH(4)
DEPRN=1.0D+01*DRELPR
EPSTAB(N+2)=EPSTAB(N)
NEWELM=(N-1)/2
NUM=N
K1=N
DO 40 I=1,NEWELM
K2=K1-1
K3=K1-2
RES=EPSTAB(K1+2)
E0=EPSTAB(K3)
E1=EPSTAB(K2)
E2=RES
DELTA2=E2-E1
ERR2=ABS(DELTA2)
TOL2=MAX(ABS(E2),ABS(E1))*DRELPR
DELTA3=E1-E0
ERR3=ABS(DELTA3)
TOL3=MAX(ABS(E1),ABS(E0))*DRELPR
IF(ERR2.GT.TOL2.OR.ERR3.GT.TOL3) GO TO 10
C
C IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
C ACCURACY, CONVERGENCE IS ASSUMED.
C RESULT=E2
C ABSERR=ABS(E1-E0)+ABS(E2-E1)
C
RESULT=RES
ABSERR=ERR2+ERR3
GO TO 50
10 IF(I.NE.1) THEN
E3=EPSTAB(K1)
EPSTAB(K1)=E1
DELTA1=E1-E3
ERR1=ABS(DELTA1)
TOL1=MAX(ABS(E1),ABS(E3))*DRELPR
C
C IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
C A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
C
IF(ERR1.LE.TOL1.OR.ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20
SS=0.1D+01/DELTA1+0.1D+01/DELTA2-0.1D+01/DELTA3
ELSE
EPSTAB(K1)=E1
IF(ERR2.LE.TOL2.OR.ERR3.LE.TOL3) GO TO 20
SS=0.1D+01/DELTA2-0.1D+01/DELTA3
ENDIF
C
C TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
C EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
C OF N
C
IF(ABS(SS*E1).GT.0.1D-03) GO TO 30
20 N=I+I-1
IF(NRES.EQ.0) THEN
ABSERR=ERR2+ERR3
RESULT=RES
ELSE IF(NRES.EQ.1) THEN
RESULT=RES3LA(1)
ELSE IF(NRES.EQ.2) THEN
RESULT=RES3LA(2)
ELSE
RESULT=RES3LA(3)
ENDIF
GO TO 50
C
C COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
C THE VALUE OF RESULT
C
30 RES=E1+0.1D+01/SS
EPSTAB(K1)=RES
K1=K1-2
IF(NRES.EQ.0) THEN
ABSERR=ERR2+ABS(RES-E2)+ERR3
RESULT=RES
GO TO 40
ELSE IF(NRES.EQ.1) THEN
ERROR=.6D+01*(ABS(RES-RES3LA(1)))
ELSE IF(NRES.EQ.2) THEN
ERROR=.2D+01*(ABS(RES-RES3LA(2))+ABS(RES-RES3LA(1)))
ELSE
ERROR=ABS(RES-RES3LA(3))+ABS(RES-RES3LA(2))
1 +ABS(RES-RES3LA(1))
ENDIF
IF(ERROR.GT.1.0D+01*ABSERR) GO TO 40
ABSERR=ERROR
RESULT=RES
40 CONTINUE
C
C COMPUTE ERROR ESTIMATE
C
IF(NRES.EQ.1) THEN
ABSERR=.6D+01*(ABS(RESULT-RES3LA(1)))
ELSE IF(NRES.EQ.2) THEN
ABSERR=.2D+01*ABS(RESULT-RES3LA(2))+ABS(RESULT-RES3LA(1))
ELSE IF(NRES.GT.2) THEN
ABSERR=ABS(RESULT-RES3LA(3))+ABS(RESULT-RES3LA(2))
1 +ABS(RESULT-RES3LA(1))
ENDIF
C
C SHIFT THE TABLE
C
50 IF(N.EQ.LIMEXP) N=2*(LIMEXP/2)-1
IB=1
IF((NUM/2)*2.EQ.NUM) IB=2
IE=NEWELM+1
DO 60 I=1,IE
IB2=IB+2
EPSTAB(IB)=EPSTAB(IB2)
IB=IB2
60 CONTINUE
IF(NUM.EQ.N) GO TO 80
IN=NUM-N+1
DO 70 I=1,N
EPSTAB(I)=EPSTAB(IN)
IN=IN+1
70 CONTINUE
C
C UPDATE RES3LA
C
80 IF(NRES.EQ.0) THEN
RES3LA(1)=RESULT
ELSE IF(NRES.EQ.1) THEN
RES3LA(2)=RESULT
ELSE IF(NRES.EQ.2) THEN
RES3LA(3)=RESULT
ELSE
RES3LA(1)=RES3LA(2)
RES3LA(2)=RES3LA(3)
RES3LA(3)=RESULT
ENDIF
90 ABSERR=MAX(ABSERR,DEPRN*ABS(RESULT))
NRES=NRES+1
100 N=N+1
EPSTAB(LIMEXP+3)=DBLE(N)
EPSTAB(LIMEXP+4)=DBLE(NRES)
EPSTAB(LIMEXP+5)=RES3LA(1)
EPSTAB(LIMEXP+6)=RES3LA(2)
EPSTAB(LIMEXP+7)=RES3LA(3)
110 RETURN
END subroutine DEA
subroutine AdaptiveIntWithBreaks(f,a,b,N,brks,epsi,iflg
$ ,abserr, val)
use AdaptiveGaussKronrod
implicit none
double precision :: f
integer, intent(in) :: N
double precision, intent(in) :: a,b,epsi
double precision, dimension(:), intent(in) :: brks
double precision, intent(out) :: abserr, val
integer, intent(out) :: iflg
external f
! Locals
double precision, dimension(N+2) :: pts
double precision :: LTol,tol, error, valk, excess, errorEstimate
double precision :: delta, deltaK
integer :: kflg, k, limit,neval
limit = 30
pts(1) = a
pts(N+2) = b
delta = b - a
do k = 2,N+1
pts(k) = minval(brks(k-1:N)) !add user supplied break points
enddo
LTol = epsi / delta
abserr = 0.0d0
val = 0.0D0
iflg = 0
do k = 1, N + 1
deltaK = pts(k+1) - pts(k)
tol = LTol * deltaK
if (deltaK < 0.5D0) then
call AdaptiveSimpson(f,pts(k),pts(k+1),tol, kflg,error,valk)
! call romberg(f,pts(k),pts(k+1),20,tol,kflg,error, valk)
else
! call AdaptiveSimpson3(f,pts(k),pts(k+1),tol,kflg,error,valk)
call dqagp(f,pts(k),pts(k+1),0,pts,tol,0.0D0,limit,valk,
* error,neval,kflg)
endif
abserr = abserr + abs(error)
errorEstimate = abserr + (b - pts(k+1)) * LTol
excess = epsi - errorEstimate
if (excess < 0.0D0 ) then
LTol = 0.1D0*LTol
elseif ( epsi < 2.0D0 * excess ) then
LTol = (epsi + excess*0.5D0) / delta
endif
val = val + valk
if (kflg>0) iflg = IOR(iflg, kflg)
end do
if (epsi<abserr) iflg = IOR(iflg, 3)
end subroutine AdaptiveIntWithBreaks
subroutine AdaptiveSimpsonWithBreaks(f,a,b,N,brks,epsi,iflg
$ ,abserr, val)
implicit none
double precision :: f
integer, intent(in) :: N
double precision, intent(in) :: a,b,epsi
double precision, dimension(:), intent(in) :: brks
double precision, intent(out) :: abserr, val
integer, intent(out) :: iflg
external f
! Locals
double precision, dimension(N+2) :: pts
double precision :: error, valk, excess, errorEstimate
double precision :: Lepsi,LTol, tol, delta, deltaK,small
integer :: kflg, k
pts(1) = a
pts(N+2) = b
delta = (b - a)
do k = 2, N+1
pts(k) = minval(brks(k-1:N)) !add user supplied break points
enddo
small = spacing(1.0D0)
Lepsi = max(epsi,small)
LTol = Lepsi / delta
abserr = 0.0d0
val = 0.0D0
iflg = 0
do k = 1, N + 1
deltaK = pts(k+1)-pts(k)
tol = LTol * deltaK
call AdaptiveSimpson(f,pts(k),pts(k+1),tol, kflg,error, valk)
abserr = abserr + abs(error)
deltaK = (b-pts(k+1))
errorEstimate = abserr + deltaK * LTol
excess = Lepsi - errorEstimate
if (excess < 0.0D0 ) then
if (deltaK>0.0d0 .and. Lepsi > abserr) then
LTol = (Lepsi - abserr) / deltaK
else
LTol = 0.1D0 * LTol
endif
elseif ( Lepsi < 5D0 * excess ) then
LTol = (Lepsi + excess) / delta
endif
val = val + valk
if (kflg>0) iflg = IOR(iflg, kflg)
end do
if (epsi<abserr) iflg = IOR(iflg, 4)
end subroutine AdaptiveSimpsonWithBreaks
subroutine AdaptiveSimpson1(f,a,b,epsi, iflg,abserr, val)
implicit none
c Numerical Analysis:
c The Mathematics of Scientific Computing
c D.R. Kincaid & E.W. Cheney
c Brooks/Cole Publ., 1990
C
! Revised by Per A. Brodtkorb 4 June 2003
! Added check on stepsize, i.e., hMin and hMax
double precision :: f
double precision, intent(in) :: a,b,epsi
double precision, intent(out) :: abserr, val
integer, intent(out) :: iflg
integer, parameter :: stackLimit = 30
double precision, dimension(6,stackLimit) :: v
external f
! Locals
double precision, parameter :: small = 1.0D-16
double precision, parameter :: hmin = 1.0D-9
double precision, parameter :: zero = 0.0D0
double precision, parameter :: zpz66666 = 0.06666666666666666666D0 !1/15
double precision, parameter :: onethird = 0.33333333333333333333D0
double precision, parameter :: half = 0.5D0
double precision, parameter :: one = 1.0D0
double precision, parameter :: two = 2.0D0
double precision, parameter :: three = 3.0D0
double precision, parameter :: four = 4.0D0
double precision :: delta,h,c,y,z,localError, correction
double precision :: abar, bbar, cbar,ybar,vbar,zbar
double precision :: hmax, S, Star, SStar, SSStar,Sdiff
integer :: k
logical :: acceptError, lastInStack
logical :: stepSizeTooSmall, stepSizeOK
hmax = 0.24D0
c
c initialize everything,
c particularly the first column vector in the stack.
c
val = zero
abserr = zero
iflg = 0
delta = b - a
h = half * delta
c = half * ( a + b )
k = 1
abar = f(a)
cbar = f(c)
bbar = f(b)
S = (abar + four * cbar + bbar) * h * onethird
v(1,1) = a
v(2,1) = h
v(3,1) = abar
v(4,1) = cbar
v(5,1) = bbar
v(6,1) = S
c
do while ((1<=k) .and. (k <= stackLimit))
c
c take the last column off the stack and process it.
c
h = half * v(2,k)
y = v(1,k) + h
z = v(1,k) + three * h
ybar = f(y)
zbar = f(z)
Star = ( v(3,k) + four * ybar + v(4,k) ) * h * onethird
SStar = ( v(4,k) + four * zbar + v(5,k) ) * h * onethird
SSStar = Star + SStar
Sdiff = (SSStar - v(6,k))
correction = Sdiff * zpz66666 !=0.066666... = 1/15.0D0
localError = abs(Sdiff) * two
! acceptError is made conservative in order to avoid premature termination
acceptError = (localError * delta <= two* epsi * h
& .or. localError < small)
lastInStack = ( stackLimit <= k)
stepSizeOK = ( h < hMax )
stepSizeTooSmall = ( h < hMin)
if (lastInStack .or. (stepSizeOK.and.acceptError)
& .or. stepSizeTooSmall ) then
! Stop subdividing interval when
! 1) accuracy is sufficient, or
! 2) interval too narrow, or
! 3) subdivided too often. (stack limit reached)
! Add partial integral and take a new vector from the bottom of the stack.
abserr = abserr + localError
val = val + SSStar + correction
k = k - 1
c
if (.not.acceptError) then
if (lastInStack) iflg = IOR(iflg,1) ! stack limit reached
if (stepSizeTooSmall) iflg = IOR(iflg,2) ! stepSize limit reached
endif
if (k <= 0) then
return
endif
else
c Subdivide the interval and create two new vectors in the stack,
c one of which overwrites the vector just processed.
vbar = v(5,k)
v(2,k) = h
v(5,k) = v(4,k)
v(4,k) = ybar
v(6,k) = Star
c
k = k + 1
v(1,k) = v(1,k-1) + two * h
v(2,k) = h
v(3,k) = v(5,k-1)
v(4,k) = zbar
v(5,k) = vbar
v(6,k) = SStar
endif
enddo ! while
end subroutine AdaptiveSimpson1
subroutine AdaptiveSimpson2(f,a,b,epsi, iflg,abserr, val)
implicit none
! by Per A. Brodtkorb 4 June 2003
! based on psudo code in chapter 7, Kincaid and Cheney (1991).
! Added check on stepsize, i.e., hMin and hMax
! Added an alternitive check on termination: this is more robust
! Reference:
! D.R. Kincaid & E.W. Cheney (1991)
! "Numerical Analysis"
! Brooks/Cole Publ., 1991
!
! C. Brezinski (1977)
! "Acceleration de la convergence en analyse numerique",
! , "Lecture Notes in Math.", vol. 584,
! Springer-Verlag, New York, 1977.
double precision :: f
double precision, intent(in) :: a,b,epsi
double precision, intent(out) :: abserr, val
integer, intent(out) :: iflg
integer, parameter :: stackLimit = 300
double precision, dimension(10,stackLimit) :: v
external f
! Locals
double precision, parameter :: zero = 0.0D0
double precision, parameter :: zpz66666 = 0.06666666666666666666D0!1/15
double precision, parameter :: zpz588 = 0.05882352941176D0 !1/17
double precision, parameter :: onethird = 0.33333333333333333333D0
double precision, parameter :: half = 0.5D0
double precision, parameter :: one = 1.0D0
double precision, parameter :: two = 2.0D0
double precision, parameter :: three = 3.0D0
double precision, parameter :: four = 4.0D0
double precision, parameter :: six = 6.0D0
double precision, parameter :: eight = 8.0D0
double precision, parameter :: ten = 10.0D0
double precision, dimension(4) :: x, fx, Sn
double precision, dimension(5) :: d4fx
double precision, dimension(55) :: EPSTAB
double precision :: small
double precision :: delta, h, h8, localError, correction
double precision :: Sn1, Sn2, Sn4, Sn1e, Sn2e, Sn4e
double precision :: Sn12, Sn24, Sn124, Sn12e, Sn24e
double precision :: hmax, hmin, dhmin, val0
double precision :: Lepsi,Ltol, excess, deltaK, errorEstimate
integer :: k, kp1, i, j,ix, numExtrap, IERR
integer, parameter :: LIMEXP = 5
logical :: acceptError, lastInStack
logical :: stepSizeTooSmall, stepSizeOK
logical :: NEWFLG !, useDEA
small = spacing(one)
! useDEA = .TRUE.
Lepsi = max(epsi,small/ten)
if (Lepsi < 1.0D-7) then
numExtrap = 1
else
numExtrap = 0
endif
numExtrap = 1
hmax = one
hmin = 1.0D-9
dhmin = 1.0D-1
c
c initialize everything,
c particularly the first column vector in the stack.
c
val = zero
abserr = zero
iflg = 0
delta = b - a
h = half * delta
Ltol = Lepsi / delta
x(1) = a
x(3) = half * ( a + b )
x(2) = half * ( a + x(3) )
x(4) = half * ( x(3) + b )
k = 1
do I = 1,4
v(I,1) = f(x(I))
enddo
v(5,1) = f(b)
Sn(1) = ( v(1,1) + four * v(3,1) + v(5,1) ) * h * onethird
h = h * half
Sn(2) = ( v(1,1) + four * v(2,1) + v(3,1) ) * h * onethird
Sn(3) = ( v(3,1) + four * v(4,1) + v(5,1) ) * h * onethird
v(6 ,1) = x(1)
v(7 ,1) = h
v(8:10,1) = Sn(1:3);
do while ((1<=k) .and. (k <= stackLimit))
!
! take the last column off the stack and process it.
!
h = half * v(7,k)
do i = 1,4
x(i) = v(6,k) + dble(2*i-1)*h
fx(i) = f(x(i))
Sn(i) = ( v(i,k) + four * fx(i) + v(i+1,k) ) * h * onethird
enddo
stepSizeOK = ( h < hMax )
lastInStack = ( stackLimit <= k)
if (lastInStack .OR. stepSizeOK) then
Sn1 = v(8,k)
Sn2 = ( v(9,k) + v(10,k) )
Sn4 = Sn(1) + Sn(2) + Sn(3) + Sn(4)
if (numExtrap>0) then
Sn12 = (Sn1 - Sn2)
Sn24 = (Sn2 - Sn4)
! Extrapolate Sn1 and Sn2:
Sn1e = Sn2 - Sn12 * zpz66666
Sn2e = Sn4 - Sn24 * zpz66666
Sn12e = ( Sn1e - Sn2e )
Sn24e = (Sn2e - Sn4)
! Sn1e = Sn2e - Sn12e * zpz66666
! Sn12e = (Sn1e - Sn2e)
Sn124 = (Sn12e - Sn24)
if ((abs(Sn124)<= hmin) .or.
& .false..and.(Sn24*Sn12e < zero)) then
! Correction based on the assumption of slowly varying fourth derivative
correction = -Sn24 * zpz588 !
else
! Correction based on assumption that the termination error
! is of the form: C*h^q
correction = -Sn24 * Sn24 / Sn124
endif
Sn4e = Sn4 + correction
! NEWFLG = .TRUE.
! CALL DEA(NEWFLG,Sn1,LIMEXP,val0,localError,EPSTAB,IERR)
! CALL DEA(NEWFLG,Sn2,LIMEXP,val0,localError,EPSTAB,IERR)
! CALL DEA(NEWFLG,Sn1e,LIMEXP,val0,localError,EPSTAB,IERR)
! CALL DEA(NEWFLG,Sn4,LIMEXP,val0,localError,EPSTAB,IERR)
! CALL DEA(NEWFLG,Sn4e,LIMEXP,val0,localError,EPSTAB,IERR)
! localError is made conservative in order to avoid premature
! termination
CALL DEA3(Sn1e,Sn2e,Sn4e,localError,val0)
!if (h>dhMin) then
!localError = max(localError,abs(correction))
!else
!val0 = Sn4e
!localError = abs(correction)*two
!endif
else
CALL DEA3(Sn1,Sn2,Sn4,localError,val0)
endif
acceptError = ( localError <= Ltol * h * eight
& .or. localError < small)
else
acceptError = .FALSE.
endif
stepSizeTooSmall = ( h < hMin)
if (lastInStack .or.
& ( stepSizeOK .and. acceptError ) .or.
& stepSizeTooSmall) then
! Stop subdividing interval when
! 1) accuracy is sufficient, or
! 2) interval too narrow, or
! 3) subdivided too often. (stack limit reached)
! Add partial integral and take a new vector from the bottom of the stack.
abserr = abserr + max(localError, ten*small*val0)
val = val + val0
k = k - 1
if (.not.acceptError) then
if (lastInStack) iflg = IOR(iflg,1) !stack limit reached
if (stepSizeTooSmall) iflg = IOR(iflg,2) !stepSize limit reached
endif
if (k <= 0) then
exit ! while loop
endif
deltaK = (v(6,k+1)-a)
errorEstimate = abserr + deltaK * Ltol
excess = Lepsi - errorEstimate
if (excess < zero ) then
if (deltaK > zero .and. Lepsi > abserr) then
LTol = (Lepsi - abserr) / deltaK
else
LTol = 0.1D0 * LTol
endif
elseif (.true..or. Lepsi < four * excess ) then
LTol = (Lepsi + 0.9D0 * excess) / delta
endif
else
! Subdivide the interval and create two new vectors in the stack,
! one of which overwrites the vector just processed.
!
! v(:,k) = [fx1,fx2,fx3,fx4,fx5,x1,h,S,SL,SR]
kp1 = k + 1;
! Process right interval
v(1,kp1) = v(3,k); !fx1R
v(2,kp1) = fx(3); !fx2R
v(3,kp1) = v(4,k); !fx3R
v(4,kp1) = fx(4); !fx4R
v(5,kp1) = v(5,k); !fx5R
v(6,kp1) = v(6,k) + four * h; ! x1R
v(7,kp1) = h;
v(8,kp1) = v(10,k); ! S
v(9:10,kp1) = Sn(3:4); ! SL, SR
! Process left interval
v(5,k) = v(3,k); ! fx5L
v(4,k) = fx(2); ! fx4L
v(3,k) = v(2,k); ! fx3L
v(2,k) = fx(1); ! fx2L
! v(1,k) unchanged fx1L
! v(6,k) unchanged x1L
v(7,k) = h;
v(8,k) = v(9,k); ! S
v(9:10,k) = Sn(1:2); ! SL, SR
k = kp1;
endif
enddo ! while
if (epsi<abserr) iflg = IOR(iflg,4)
end subroutine AdaptiveSimpson2
subroutine AdaptiveSimpson3(f,a,b,epsi, iflg,abserr, val)
implicit none
! by Per A. Brodtkorb 4 June 2003
! based on psudo code in chapter 7, Kincaid and Cheney (1991).
! Added check on stepsize, i.e., hMin and hMax
! Added an alternitive check on termination: this is more robust
! Reference:
! D.R. Kincaid & E.W. Cheney (1991)
! "Numerical Analysis"
! Brooks/Cole Publ., 1991
!
! C. Brezinski (1977)
! "Acceleration de la convergence en analyse numerique",
! , "Lecture Notes in Math.", vol. 584,
! Springer-Verlag, New York, 1977.
double precision :: f
double precision, intent(in) :: a,b,epsi
double precision, intent(out) :: abserr, val
integer, intent(out) :: iflg
integer, parameter :: stackLimit = 10
double precision, dimension(18,stackLimit) :: v
external f
! Locals
double precision, parameter :: zero = 0.0D0
double precision, parameter :: zp125 = 0.125D0 ! 1/8
double precision, parameter :: zpz66666 = 0.06666666666666666666D0!1/15
double precision, parameter :: zpz588 = 0.05882352941176D0 !1/17
double precision, parameter :: onethird = 0.33333333333333333333D0
double precision, parameter :: half = 0.5D0
double precision, parameter :: one = 1.0D0
double precision, parameter :: two = 2.0D0
double precision, parameter :: three = 3.0D0
double precision, parameter :: four = 4.0D0
double precision, parameter :: six = 6.0D0
double precision, parameter :: eight = 8.0D0
double precision, parameter :: ten = 10.0D0
double precision, parameter :: sixteen = 16.0D0
double precision, dimension(8) :: fx, Sn
double precision, dimension(55) :: EPSTAB
double precision :: small, x, a0
double precision :: delta, h, h8, localError, correction
double precision :: Sn1, Sn2, Sn4, Sn8, Sn1e, Sn2e, Sn4e
double precision :: Sn12, Sn24,Sn48,Sn124, Sn12e, Sn24e
double precision :: hmax, hmin, dhmin, val0
double precision :: Lepsi,Ltol, excess, deltaK, errorEstimate
integer :: k, kp1, i, j,ix, numExtrap, IERR, Nrule, step
integer, parameter :: LIMEXP = 5
logical :: acceptError, lastInStack
logical :: stepSizeTooSmall, stepSizeOK
logical :: NEWFLG !, useDEA
Nrule = 9
small = spacing(one)
! useDEA = .TRUE.
Lepsi = max(epsi,small/ten)
if (Lepsi < 1.0D-5) then
numExtrap = 1
else
numExtrap = 0
endif
numExtrap = 0
hmax = one
hmin = 1.0D-9
dhmin = 1.0D-1
c
c initialize everything,
c particularly the first column vector in the stack.
c
localError = one
val = zero
abserr = zero
iflg = 0
delta = b - a
h = delta * zp125
Ltol = Lepsi / delta
k = 1
do I = 1,Nrule-1
x = a + dble(i-1) * h
v(I,1) = f(x)
enddo
v(Nrule, 1) = f(b)
v(Nrule+1, 1) = a
v(Nrule+2, 1) = h
step = 8
ix = Nrule + 2
do I = 1,3
step = step / 2
do j = 1,Nrule-2,2*step
ix = ix + 1
v(ix,1) = ( v(j,1) + four * v(j+step,1) + v(j + 2*step,1) )*
& h * onethird * dble(step)
enddo
enddo
do while ((1<=k) .and. (k <= stackLimit))
!
! take the last column off the stack and process it.
!
h = half * v(Nrule+2,k)
a0 = v(Nrule + 1,k)
Sn8 = zero
do i = 1, Nrule - 1
x = a0 + dble(2*i-1)*h
fx(i) = f(x)
Sn(i) = ( v(i,k) + four * fx(i) + v(i+1,k) ) * h * onethird
Sn8 = Sn8 + Sn(i)
enddo
ix = Nrule + 3
Sn1 = v(ix,k)
Sn2 = v(ix+1,k) + v(ix+2,k)
Sn4 = v(ix+3,k) + v(ix+4,k) + v(ix+5,k) + v(ix+6,k)
stepSizeOK = ( h < hMax )
lastInStack = ( stackLimit <= k )
if (lastInStack .OR. stepSizeOK) then
if (numExtrap>0) then
Sn12 = (Sn1 - Sn2)
Sn24 = (Sn2 - Sn4)
Sn48 = (Sn4 - Sn8)
! Extrapolate Sn1 and Sn2:
Sn1e = Sn2 - Sn12 * zpz66666
Sn2e = Sn4 - Sn24 * zpz66666
Sn4e = Sn8 - Sn48 * zpz588
Sn12e = (Sn1e - Sn2e)
Sn24e = (Sn2e - Sn4e)
Sn124 = (Sn12e - Sn24e)
if ((abs(Sn124)<= hmin) .or.
& (Sn12e*Sn24e < zero)) then
! Correction based on the assumption of slowly varying fourth derivative
correction = -Sn48*zpz588 !
else
! Correction based on assumption that the termination error
! is of the form: C*h^q
correction = -Sn24e * Sn24e / Sn124
!Sn4e = Sn4e + correction
endif
CALL DEA3(Sn1e,Sn2e,Sn4e,localError,val0)
! localError is made conservative in order to avoid premature
! termination
! localError = max(localError,abs(correction)*three)
! localError = abs(correction)*three
else
!CALL DEA3(Sn1,Sn2,Sn4,localError,val0)
NEWFLG = .TRUE.
CALL DEA(NEWFLG,Sn1,LIMEXP,val0,localError,EPSTAB,IERR)
CALL DEA(NEWFLG,Sn2,LIMEXP,val0,localError,EPSTAB,IERR)
CALL DEA(NEWFLG,Sn4,LIMEXP,val0,localError,EPSTAB,IERR)
CALL DEA(NEWFLG,Sn8,LIMEXP,val0,localError,EPSTAB,IERR)
endif
acceptError = ( localError <= Ltol * h * sixteen
& .or. localError < small)
else
acceptError = .FALSE.
endif
stepSizeTooSmall = ( h < hMin)
if (lastInStack .or.
& ( stepSizeOK .and. acceptError ) .or.
& stepSizeTooSmall) then
! Stop subdividing interval when
! 1) accuracy is sufficient, or
! 2) interval too narrow, or
! 3) subdivided too often. (stack limit reached)
! Add partial integral and take a new vector from the bottom of the stack.
abserr = abserr + max(localError, ten*small*val0)
val = val + val0
k = k - 1
if (.not.acceptError) then
if (lastInStack) iflg = IOR(iflg,1) !stack limit reached
if (stepSizeTooSmall) iflg = IOR(iflg,2) !stepSize limit reached
endif
if (k <= 0) then
exit ! while loop
endif
deltaK = (v(Nrule+1,k+1)-a)
errorEstimate = abserr + deltaK * Ltol
excess = Lepsi - errorEstimate
if (excess < zero ) then
if (deltaK > zero .and. Lepsi > abserr) then
LTol = (Lepsi - abserr) / deltaK
else
LTol = 0.1D0 * LTol
endif
elseif (.TRUE..or. Lepsi < four * excess ) then
LTol = (Lepsi + 0.9D0 * excess) / delta
endif
else
! Subdivide the interval and create two new vectors in the stack,
! one of which overwrites the vector just processed.
!
! v(:,k) = [fx1,fx2,..,fx8,fx9,x1,h,S,SL,SR,SL1,SL2 SR1,SR2]
kp1 = k + 1;
! Process right interval
v(1,kp1) = v(5,k); !fx1R
v(2,kp1) = fx(5); !fx2R
v(3,kp1) = v(6,k); !fx3R
v(4,kp1) = fx(6); !fx4R
v(5,kp1) = v(7,k); !fx5R
v(6,kp1) = fx(7); !fx6R
v(7,kp1) = v(8,k); !fx7R
v(8,kp1) = fx(8); !fx8R
v(9,kp1) = v(9,k); !fx9R
v(Nrule+1,kp1) = v(Nrule+1,k) + eight * h ! x1R
v(Nrule+2,kp1) = h;
v(Nrule+3,kp1) = v(Nrule+5,k); ! S
v(Nrule+4,kp1) = v(Nrule+8,k); ! SL
v(Nrule+5,kp1) = v(Nrule+9,k); ! SR
v(Nrule+6:Nrule+9,kp1) = Sn(5:8); ! SL1,SL2,SR1, SR2
! Process left interval
v(9,k) = v(5,k); ! fx9L
v(8,k) = fx(4); ! fx8L
v(7,k) = v(4,k); ! fx7L
v(6,k) = fx(3); ! fx6L
v(5,k) = v(3,k); ! fx5L
v(4,k) = fx(2); ! fx4L
v(3,k) = v(2,k); ! fx3L
v(2,k) = fx(1); ! fx2L
! v(1,k) = v(1,k); ! fx1L
! v(Nrule+1,k) unchanged x1L
v(Nrule+2,k) = h;
v(Nrule+3,k) = v(Nrule + 4,k); ! S
v(Nrule+4,k) = v(Nrule+6,k); ! SL
v(Nrule+5,k) = v(Nrule+7,k); ! SR
v(Nrule+6:Nrule+9,k) = Sn(1:4); ! SL1,SL2,SR1, SR2
k = kp1;
endif
enddo ! while
if (epsi<abserr) iflg = IOR(iflg,4)
end subroutine AdaptiveSimpson3
subroutine AdaptiveTrapzWithBreaks(f,a,b,N,brks,epsi,iflg
$ ,abserr, val)
implicit none
double precision :: f
integer, intent(in) :: N
double precision, intent(in) :: a,b,epsi
double precision, dimension(:), intent(in) :: brks
double precision, intent(out) :: abserr, val
integer, intent(out) :: iflg
external f
! Locals
double precision, dimension(N+2) :: pts
double precision :: LTol,tol, error, valk, excess, errorEstimate
integer :: kflg, k
pts(1) = a
pts(N+2) = b
do k = 2,N+1
pts(k) = minval(brks(k-1:N)) !add user supplied break points
enddo
LTol = epsi / dble( N+2 )
tol = LTol
abserr = 0.0d0
val = 0.0D0
iflg = 0
do k = 1, N + 1
call AdaptiveTrapz(f,pts(k),pts(k+1),tol, kflg,error, valk)
abserr = abserr + abs(error)
excess = LTol - abs(error)
errorEstimate = abserr + dble(N-k+1)*LTol
if (epsi < errorEstimate ) then
tol = max(0.1D0*LTol,2.0D-16)
! elseif ( LTol < excess / 10.0D0 ) then
! tol = LTol + excess*0.5D0
else
tol = LTol
endif
val = val + valk
if (kflg>0) iflg = IOR(iflg, kflg)
end do
if (epsi<abserr) iflg = IOR(iflg, 3)
end subroutine AdaptiveTrapzWithBreaks
subroutine AdaptiveTrapz1(f,a,b,epsi, iflg,abserr, val)
implicit none
! by Per A. Brodtkorb 4 June 2003
! based on psudo code in chapter 7, Kincaid and Cheney (1991).
! Added check on stepsize, i.e., hMin and hMax
! Added an alternitive check on termination: this is more robust
! Reference:
! D.R. Kincaid & E.W. Cheney (1991)
! "Numerical Analysis"
! Brooks/Cole Publ., 1991
!
double precision :: f
double precision, intent(in) :: a,b,epsi
double precision, intent(out) :: abserr, val
integer, intent(out) :: iflg
integer, parameter :: stackLimit = 30
double precision, dimension(8,stackLimit) :: v
external f
! Locals
double precision, parameter :: small = 1.0D-16
double precision, parameter :: hmin = 1.0D-10
double precision, parameter :: zero = 0.0D0
double precision, parameter :: zpz66666 = 0.06666666666666666666D0 !1/15
double precision, parameter :: onethird = 0.33333333333333333333D0
double precision, parameter :: half = 0.5D0
double precision, parameter :: one = 1.0D0
double precision, parameter :: two = 2.0D0
double precision, parameter :: three = 3.0D0
double precision, parameter :: four = 4.0D0
double precision, dimension(4) :: x, fx, Sn
double precision :: delta,h,localError, correction
double precision :: hmax, Sn1, Sn2, Sn4, Sn12, Sn24, Sn124
integer :: k, kp1, i
logical :: acceptError, lastInStack
logical :: stepSizeTooSmall, stepSizeOK
hmax = 0.24D0
c
c initialize everything,
c particularly the first column vector in the stack.
c
val = zero
abserr = zero
iflg = 0
delta = b - a
h = delta
x(1) = a
x(2) = half * ( a + b )
x(3) = b
k = 1
do I = 1,3
v(I,1) = f(x(I))
enddo
Sn(1) = ( v(1,1) + v(3,1) ) * h * half
h = h * half
Sn(2) = ( v(1,1) + v(2,1) ) * h * half
Sn(3) = ( v(2,1) + v(3,1) ) * h * half
v(4 ,1) = x(1)
v(5 ,1) = h
v(6:8 ,1) = Sn(1:3);
do while ((1<=k) .and. (k <= stackLimit))
!
! take the last column off the stack and process it.
!
h = half * v(5,k)
Sn1 = v(6,k)
Sn2 = (v(7,k) + v(8,k))
Sn4 = zero
do i = 1,2
x(i) = v(4,k) + dble(2*i-1)*h
fx(i) = f(x(i))
Sn(2*i-1) = ( v(i ,k) + fx(i) ) * h * half
Sn(2*i ) = ( v(i+1,k) + fx(i) ) * h * half
Sn4 = Sn4 + Sn(2*i-1) + Sn(2*i)
enddo
Sn12 = (Sn1 - Sn2)
Sn24 = (Sn2 - Sn4);
Sn124 = (Sn12 - Sn24)
! Correction based on assumption that the termination error
! is of the form: C*h^q
if (Sn124==zero) then
correction = zero !-Sn24 * Sn24 / sign(small,Sn124)
if (Sn12 == zero) then
!round off error?
endif
else
correction = -Sn24 * Sn24 / Sn124
endif
localError = max(abs(correction),abs(Sn24)*half)
! acceptError is made conservative in order to avoid premature termination
acceptError = (localError * delta <= two * epsi * h
& .or. localError < small)
lastInStack = ( stackLimit <= k)
stepSizeOK = ( h < hMax )
stepSizeTooSmall = ( h < hMin)
if (lastInStack .or. (stepSizeOK.and.acceptError)
& .or. stepSizeTooSmall) then
! Stop subdividing interval when
! 1) accuracy is sufficient, or
! 2) interval too narrow, or
! 3) subdivided too often. (stack limit reached)
! Add partial integral and take a new vector from the bottom of the stack.
abserr = abserr + localError
val = val + Sn4 + correction
k = k - 1
if (.not.acceptError) then
if (lastInStack) iflg = IOR(iflg,1) ! stack limit reached
if (stepSizeTooSmall) iflg = IOR(iflg,2) ! stepSize limit reached
endif
if (k <= 0) then
return
endif
else
! Subdivide the interval and create two new vectors in the stack,
! one of which overwrites the vector just processed.
!
! v(:,k) = [fx1,fx2,fx3,x1,h,S,SL,SR]
kp1 = k + 1;
! Process right interval
v(1,kp1) = v(2,k); !fx1R
v(2,kp1) = fx(2); !fx2R
v(3,kp1) = v(3,k); !fx3R
v(4,kp1) = v(4,k) + two * h; ! x1R
v(5,kp1) = h;
v(6,kp1) = v(8,k); ! S
v(7:8,kp1) = Sn(3:4); ! SL, SR
! Process left interval
v(3,k) = v(2,k); ! fx5L
v(2,k) = fx(1); ! fx4L
! v(1,k) unchanged fx1L
! v(4,k) unchanged x1L
v(5,k) = h;
v(6,k) = v(7,k); ! S
v(7:8,k) = Sn(1:2); ! SL, SR
k = kp1;
endif
enddo ! while
end subroutine AdaptiveTrapz1
subroutine RombergWithBreaks(f,a,b,N,brks,epsi,iflg
$ ,abserr, val)
implicit none
double precision :: f
integer, intent(in) :: N
double precision, intent(in) :: a,b,epsi
double precision, dimension(:), intent(in) :: brks
double precision, intent(out) :: abserr, val
integer, intent(out) :: iflg
external f
! Locals
double precision, dimension(N+2) :: pts
double precision :: tol, LTol, error, valk
double precision :: excess, errorEstimate
double precision :: delta, deltaK
integer :: kflg, k, deciDig
pts(1) = a
pts(N+2) = b
delta = (b - a) + 0.1D0
do k = 2,N + 1
pts(k) = minval(brks(k-1:N)) !add user supplied break points
enddo
LTol = epsi / delta
decidig = max(abs(NINT(log10(epsi)))+5,3)
abserr = 0.0d0
val = 0.0D0
iflg = 0
do k = 1, N+1
deltaK = pts( k + 1 ) - pts( k )
tol = LTol * deltaK
call romberg(f,pts(k),pts(k+1),decidig,tol,kflg,error, valk)
abserr = abserr + abs(error)
errorEstimate = abserr + (b - pts( k + 1 ) )* LTol
excess = epsi - errorEstimate
if (excess < 0.0D0 ) then
LTol = 0.1D0 * LTol
elseif ( epsi < 2.0D0 * excess ) then
LTol = (epsi + excess*0.5D0) / delta
endif
val = val + valk
if (kflg>0) iflg = ior(iflg,kflg)
end do
if (epsi<abserr) iflg = IOR(iflg, 3)
end subroutine RombergWithBreaks
subroutine Romberg1(f,a,b,decdigs,abseps,errFlg,abserr,VAL)
implicit none
double precision :: f
double precision, intent(in) :: a,b, abseps
integer, intent(in) :: decdigs ! Relative number of decimal digits
integer, intent(out) :: errFlg
double precision, intent(out) :: abserr, val
external f
! Locals
double precision, dimension(decdigs) :: rom1,rom2,fp
double precision, dimension(decdigs + 7) :: EPSTAB
integer :: LIMEXP
double precision :: h, fourk, Un5, T1n, T2n, T4n, T12n, T24n
double precision :: correction
integer :: ipower, k, i, IERR
logical :: stepSizeTooSmall, NEWFLG
double precision, parameter :: four = 4.0D0
double precision, parameter :: one = 1.0D0
double precision, parameter :: half = 0.5D0
double precision, parameter :: zero = 0.0D0
double precision, parameter :: hmin = 1.0d-10
LIMEXP = decdigs
val = zero
errFlg = 0
rom1(:) = zero
rom2(:) = zero
fp(1) = four;
h = ( b - a )
rom1(1) = h * ( f(a) + f(b) ) * half
ipower = 1
T1n = zero
T2n = zero
T4n = rom1(1)
NEWFLG = .TRUE.
CALL DEA(NEWFLG,T4n,LIMEXP,val,abserr,EPSTAB,IERR)
stepSizeTooSmall = ( h < hMin)
do i = 2, decdigs
h = h * half
Un5 = zero
do k = 1, ipower
Un5 = Un5 + f(a + DBLE(2*k-1)*h)
enddo
! trapezoidal approximations
rom2(1) = half * rom1(1) + h * Un5
! Richardson extrapolation
do k = 1, i-1
rom2(k+1) = ( fp(k)*rom2(k)-rom1(k) ) / ( fp(k) - one )
enddo
T1n = T2n
T2n = T4n
T4n = rom2(i)
CALL DEA(NEWFLG,T4n,LIMEXP,val,abserr,EPSTAB,IERR)
if (3<=i) then
! T12n = ( T1n - T2n )
! T24n = ( T2n - T4n )
! correction = -T24n * T24n / ( T12n - T24n )
! abserr = abs( correction )
stepSizeTooSmall = ( h < hMin)
if ( abserr < abseps .or. stepSizeTooSmall) then
exit ! exit do loop
endif
endif
rom1(1:i) = rom2(1:i)
ipower = ipower * 2
fp(i) = four * fp(i-1)
enddo
if (decdigs < 3) then
! val = T4n
abserr = min(abs(T4n-T2n) * half, abserr)
endif
if (abseps < abserr) then
if (stepSizeTooSmall) errFlg = ior(errFlg,2) ! step size limit reached
endif
end subroutine Romberg1
end module Integration1DModule
module mvnProdCorrPrbMod
implicit none
private
public :: mvnprodcorrprb
double precision, parameter :: mINFINITY = 8.25D0 !
! Inputs to integrand
integer mNdim ! # of mRho/=0 and mRho/=+/-1 and -inf<a or b<inf
double precision, allocatable, dimension(:) :: mRho, mDen
double precision, allocatable, dimension(:) :: mA,mB
INTERFACE mvnprodcorrprb
MODULE PROCEDURE mvnprodcorrprb
END INTERFACE
INTERFACE FI
MODULE PROCEDURE FI
END INTERFACE
INTERFACE FI2
MODULE PROCEDURE FI2
END INTERFACE
INTERFACE FIINV
MODULE PROCEDURE FIINV
END INTERFACE
INTERFACE GetBreakPoints
MODULE PROCEDURE GetBreakPoints
END INTERFACE
INTERFACE NarrowLimits
MODULE PROCEDURE NarrowLimits
END INTERFACE
INTERFACE GetTruncationError
MODULE PROCEDURE GetTruncationError
END INTERFACE
INTERFACE integrand
MODULE PROCEDURE integrand
END INTERFACE
INTERFACE integrand1
MODULE PROCEDURE integrand1
END INTERFACE
contains
SUBROUTINE SORTRE(rarray,indices)
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(:), INTENT(inout) :: rarray
INTEGER, DIMENSION(:), OPTIONAL, INTENT(inout) :: indices
! local variables
double precision :: tmpR
INTEGER :: i,im,j,k,m,n, tmpI
! diminishing increment sort as described by
! Donald E. Knuth (1973) "The art of computer programming,",
! Vol. 3, pp 84- (sorting and searching)
n = size(rarray)
! if (present(indices)) then
! if the below is commented out then assume indices are already initialized
! forall(i=1,n) indices(i) = i
! endif
100 continue
if (n.le.1) goto 800
m=1
200 continue
m=m+m
if (m.lt.n) goto 200
m=m-1
300 continue
m=m/2
if (m.eq.0) goto 800
k=n-m
j=1
400 continue
i=j
500 continue
im=i+m
if (rarray(i).gt.rarray(im)) goto 700
600 continue
j=j+1
if (j.gt.k) goto 300
goto 400
700 continue
tmpR = rarray(i)
rarray(i) = rarray(im)
rarray(im) = tmpR
if (present(indices)) then
tmpI = indices(i)
indices(i) = indices(im)
indices(im) = tmpI
endif
i=i-m
if (i.lt.1) goto 600
goto 500
800 continue
RETURN
END SUBROUTINE SORTRE
subroutine mvnprodcorrprb(rho,a,b,abseps,releps,useBreakPoints,
& useSimpson,abserr,errFlg,prb)
use AdaptiveGaussKronrod
use Integration1DModule
! use numerical_libraries
implicit none
double precision,dimension(:),intent(in) :: rho,a,b
double precision,intent(in) :: abseps,releps
logical, intent(in) :: useBreakPoints,useSimpson
double precision,intent(out) :: abserr,prb
integer, intent(out) :: errFlg
! Locals
double precision, parameter :: ZERO = 0.0D0
double precision, parameter :: ZPT1 = 0.1D0
double precision, parameter :: ZPTZ5 = 0.05D0
double precision, parameter :: ZPTZZ1 = 0.001D0
double precision, parameter :: ZPTZZZ1 = 0.0001D0
double precision, parameter :: ONE = 1.d0
double precision :: small, LTol, val0,val, truncError
double precision :: zCutOff, zlo, zup, As, Bs
double precision, dimension(1000) :: breakPoints
integer :: n, k , limit, Npts, neval
logical :: isSingular, isLimitsNarrowed
small = MAX(spacing(one),1.0D-16)
isSingular = .FALSE.
n = size(a,DIM=1)
LTol = max(abseps,small)
errFlg = 0
prb = ZERO
abserr = small
if ( any(b(:)<=a(:)).or.
& any(b(:)<=-mINFINITY) .or.
& any(mINFINITY<=a(:))) then
goto 999 ! end program
endif
As = - mInfinity
Bs = mInfinity
zCutOff = abs(max(FIINV(ZPTZ5*LTol),-mINFINITY));
zlo = - zCutOff
zup = zCutOff
allocate(mA(n),mB(n),mRho(n),mDen(n))
do k = 1,n
if (one <= abs(rho(k)) ) then
mRho(k) = sign(one,rho(k))
mDen(k) = zero
else
mRho(k) = rho(k)
mDen(k) = sqrt(one - rho(k))*sqrt(one + rho(k))
endif
end do
! See if we may narrow down the integration region: zlo, zup
CALL NarrowLimits(zlo,zup,As,Bs,zCutOff,n,a,b,mRho,mDen)
if (zup <= zlo) goto 999 ! end program
! Move only significant variables to mA,mB, and mRho
! (Note: If you scale them with mDen, the integrand must also be changed)
mNdim = 0
val0 = one
do k = 1, n
if (small < abs(mRho(k))) then
if ( ONE <= abs(mRho(k))) then
! rho(k) == 1
isSingular = .TRUE.
elseif ((-mINFINITY < a(k)) .OR. (b(k) < mINFINITY)) then
mNdim = mNdim + 1
mA(mNdim) = a(k) / mDen(k)
mB(mNdim) = b(k) / mDen(k)
mRho(mNdim) = mRho(k) / mDen(k)
mDen(mNdim) = mDen(k)
endif
else ! independent variables which are evaluated separately
val0 = val0 * ( FI( b(k) ) - FI( a(k) ) )
endif
enddo
CALL GetTruncationError(zlo, zup, As, Bs, truncError)
select case(mNdim)
case (0)
if (isSingular) then
prb = ( FI( zup ) - FI( zlo ) ) * val0
abserr = sqrt(small) + truncError
else
prb = val0;
abserr = small+truncError;
endif
goto 999 ! end program
case (1)
if (.not.isSingular) then
prb = (FI(mB(1)*mDen(1))-FI(mA(1)*mDen(1))) * val0
abserr = small + truncError
goto 999 ! end program
endif
end select
if (small < val0) then
isLimitsNarrowed = ((-7.D0 < zlo) .or. (zup < 7.D0))
Npts = 0
if (isLimitsNarrowed.AND.useBreakPoints) then
! Provide abscissas for break points
CALL GetBreakPoints(zlo,zup,mNdim,mA,mB,mRho,mDen,
& breakPoints,Npts)
endif
LTol = LTol - truncError
!
if (useSimpson) then
call AdaptiveSimpson(integrand,zlo,zup,Npts,breakPoints,LTol
& ,errFlg,abserr, val)
else
limit = 100
call dqagp(integrand,zlo,zup,Npts,breakPoints,LTol,zero,
& limit,val,abserr,neval,errFlg)
endif
prb = val * val0;
abserr = (abserr + truncError)* val0;
else
prb = zero
abserr = small + truncError
endif
999 continue
if (allocated(mDen)) deallocate(mDen)
if (allocated(mA)) deallocate(mA,mB,mRho)
return
end subroutine mvnprodcorrprb
subroutine GetTruncationError(zlo, zup, As, Bs, truncError)
double precision, intent(in) :: zlo, zup, As, Bs
double precision, intent(out) :: truncError
double precision :: upError,loError
! Computes the upper bound for the truncation error
upError = integrand1(zup) * abs( FI( Bs ) - FI( zup ) )
loError = integrand1(zlo) * abs( FI( zlo ) - FI( As ) )
truncError = loError + upError
end subroutine GetTruncationError
subroutine GetBreakPoints(xlo,xup,n,a,b,rho,den,
& breakPoints,Npts)
implicit none
double precision, intent(in) :: xlo, xup
double precision,dimension(:), intent(in) :: a,b, rho,den
integer, intent(in) :: n
double precision,dimension(:), intent(inout) :: breakPoints
integer, intent(inout) :: Npts
! Locals
integer, dimension(2*n) :: indices
integer, dimension(4*n) :: indices2
double precision, dimension(2*n) :: brkPts
double precision, dimension(4*n) :: brkPtsVal
double precision, parameter :: zero = 0.0D0, brkSplit = 2.5D0
double precision, parameter :: stepSize = 0.24
double precision :: brk,brk1,hMin,distance, xLow, dx
double precision :: z1, z2, val1,val2
integer :: j,k, kL,kU , Nprev, Nk
hMin = 1.0D-5
kL = 0
Npts = 0
if (.false.) then
if (xup-xlo>stepSize) then
Nk = floor((xup-xlo)/stepSize) + 1
dx = (xup-xlo)/dble(Nk)
do j=1, Nk -1
Npts = Npts + 1
breakPoints(Npts) = xlo + dx * dble( j )
enddo
endif
else
! Compute candidates for the breakpoints
brkPts(1:2*n) = xup
forall(k=1:n,rho(k) .ne. zero)
indices(2*k-1) = k
indices(2*k ) = k
brkPts(2*k-1) = a(k)/rho(k)
brkPts(2*k ) = b(k)/rho(k)
end forall
! Sort the candidates
call sortre(brkPts,indices)
! Make unique list of breakpoints
do k = 1,2*n
brk = brkPts(k)
if (xlo < brk) then
if ( xup <= brk ) exit ! terminate do loop
! if (Npts>0) then
! xLow = max(xlo, breakPoints(Npts))
! else
! xLow = xlo
! endif
! if (brk-xLow>stepSize) then
! Nk = floor((brk-xLow)/stepSize)
! dx = (brk-xLow)/dble(Nk)
! do j=1, Nk -1
! Npts = Npts + 1
! breakPoints(Npts) = brk + dx * dble( j )
! enddo
! endif
kU = indices(k)
!if ( xlo + distance < brk .and. brk + distance < xup )
!then
if ( den(kU) < 0.2) then
distance = max(brkSplit*den(kU),hMin)
z1 = brk + distance
z2 = brk - distance
if (Npts <= 0) then
if (xlo + distance < z1) then
Npts = Npts + 1
breakPoints(Npts) = z1
brkPtsVal(Npts) = integrand(z1)
indices2(Npts) = kU
endif
! Nprev = Nprev + 1
! breakPoints(Npts + Nprev) = brk
if ( z2 + distance < xup) then
Npts = Npts + 1
breakPoints(Npts) = z2
brkPtsVal(Npts) = integrand(z2)
indices2(Npts) = kU
endif
kL = kU
elseif (breakPoints(Npts)+ max(distance
& ,brkSplit*den(kL)) < z1) then
if (breakPoints(Npts) + distance < z1) then
Npts = Npts + 1
breakPoints(Npts) = z1
brkPtsVal(Npts) = integrand(z1)
indices2(Npts) = kU
kL = kU
endif
! Nprev = Nprev + 1
! breakPoints(Npts + Nprev) = brk
if ( z2 + distance < xup) then
Npts = Npts + 1
breakPoints(Npts) = z2
brkPtsVal(Npts) = integrand(z2)
indices2(Npts) = kU
kL = kU
endif
else
val1 = 0.0d0
val2 = 0.0d0
brkPts(Npts+1) = integrand(z1)
brkPts(Npts+2) = integrand(z2)
if ((xlo+ distance < z1) .and. (z1 + distance < xup))
& val2 = brkPts(Npts +1)
if ((xlo+ distance < z2) .and. (z2 + distance < xup))
& val2 = max(val2,brkPts(Npts +2))
val1 = breakPoints(Npts)
Nprev = 1
if (Npts>1) then
if (indices2(Npts-1)==kL) then
Nprev = 2
val1 = max(val1,breakPoints(Npts-1))
endif
endif
if (val1 < val2) then
!overwrite previous candidate
Npts = Npts - Nprev
if (Npts>0) then
val1 = breakPoints(Npts)+ distance
else
val1 = xlo+ distance
endif
if (val1 < z1) then
Npts = Npts + 1
breakPoints(Npts) = z1
brkPtsVal(Npts) = brkPtsVal(Npts+Nprev)
indices2(Npts) = kU
endif
! Nprev = Nprev + 1
! breakPoints(Npts + Nprev) = brk
if ((val1< z2) .and. (z2 + distance < xup)) then
Npts = Npts + 1
breakPoints(Npts) = z2
brkPtsVal(Npts) = integrand(z2)
indices2(Npts) = kU
endif
if (Npts>0) kL = indices2(Npts)
endif
endif
endif
endif
enddo
endif
end subroutine GetBreakPoints
subroutine NarrowLimits(zMin,zMax,As,Bs,zCutOff,n,a,b,rho,den)
implicit none
double precision, intent(inout) :: zMin, zMax, As, Bs
double precision,dimension(*),intent(in) :: rho,a,b,den
double precision, intent(in) :: zCutOff
integer, intent(in) :: n
! Locals
double precision, parameter :: zero = 0.0D0, one = 1.0D0
integer :: k
! Uses the regression equation to limit the
! integration limits zMin and zMax
do k = 1,n
if (ZERO < rho(k)) then
zMax = max(zMin, min(zMax,(b(k)+den(k)*zCutOff)/rho(k)))
zMin = min(zMax, max(zMin,(a(k)-den(k)*zCutOff)/rho(k)))
if ( one <= rho(k) ) then
if ( b(k) < Bs ) Bs = b(k)
if ( As < a(k) ) As = a(k)
endif
elseif (rho(k)< ZERO) then
zMax = max(zMin,min(zMax,(a(k)-den(k)*zCutOff)/rho(k)))
zMin = min(zMax,max(zMin,(b(k)+den(k)*zCutOff)/rho(k)))
if ( rho(k) <= -one ) then
if ( -a(k) < Bs ) Bs = -a(k)
if ( As < -b(k) ) As = -b(k)
endif
endif
enddo
As = min(As,Bs)
end subroutine NarrowLimits
function integrand(z) result (val)
implicit none
DOUBLE PRECISION, INTENT(IN) :: Z
DOUBLE PRECISION :: VAL
double precision, parameter :: sqtwopi1 = 0.39894228040143D0
double precision, parameter :: half = 0.5D0
val = sqtwopi1 * exp(-half * z * z) * integrand1(z)
return
end function integrand
function integrand1(z) result (val)
implicit none
double precision, intent(in) :: z
double precision :: val
double precision :: xUp,xLo,zRho
double precision, parameter :: one = 1.0D0, zero = 0.0D0
integer :: I
val = one
do I = 1, mNdim
zRho = z * mRho(I)
! Uncomment / mDen below if mRho, mA, mB is not scaled
xUp = ( mB(I) - zRho ) !/ mDen(I)
xLo = ( mA(I) - zRho ) !/ mDen(I)
if (zero<xLo) then
val = val * ( FI( -xLo ) - FI( -xUp ) )
else
val = val * ( FI( xUp ) - FI( xLo ) )
endif
enddo
end function integrand1
FUNCTION FIINV(P) RESULT (VAL)
IMPLICIT NONE
*
* ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
*
* Produces the normal deviate Z corresponding to a given lower
* tail area of P.
* Absolute error less than 1e-13
* Relative error less than 1e-15 for abs(VAL)>0.1
*
* The hash sums below are the sums of the mantissas of the
* coefficients. They are included for use in checking
* transcription.
*
DOUBLE PRECISION, INTENT(in) :: P
DOUBLE PRECISION :: VAL
!local variables
DOUBLE PRECISION SPLIT1, SPLIT2, CONST1, CONST2, ONE, ZERO, HALF,
& A0, A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7,
& C0, C1, C2, C3, C4, C5, C6, C7, D1, D2, D3, D4, D5, D6, D7,
& E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, F4, F5, F6, F7,
& Q, R
PARAMETER ( SPLIT1 = 0.425D0, SPLIT2 = 5.D0,
& CONST1 = 0.180625D0, CONST2 = 1.6D0,
& ONE = 1.D0, ZERO = 0.D0, HALF = 0.5D0 )
*
* Coefficients for P close to 0.5
*
PARAMETER (
* A0 = 3.38713 28727 96366 6080D0,
* A1 = 1.33141 66789 17843 7745D+2,
* A2 = 1.97159 09503 06551 4427D+3,
* A3 = 1.37316 93765 50946 1125D+4,
* A4 = 4.59219 53931 54987 1457D+4,
* A5 = 6.72657 70927 00870 0853D+4,
* A6 = 3.34305 75583 58812 8105D+4,
* A7 = 2.50908 09287 30122 6727D+3,
* B1 = 4.23133 30701 60091 1252D+1,
* B2 = 6.87187 00749 20579 0830D+2,
* B3 = 5.39419 60214 24751 1077D+3,
* B4 = 2.12137 94301 58659 5867D+4,
* B5 = 3.93078 95800 09271 0610D+4,
* B6 = 2.87290 85735 72194 2674D+4,
* B7 = 5.22649 52788 52854 5610D+3 )
* HASH SUM AB 55.88319 28806 14901 4439
*
* Coefficients for P not close to 0, 0.5 or 1.
*
PARAMETER (
* C0 = 1.42343 71107 49683 57734D0,
* C1 = 4.63033 78461 56545 29590D0,
* C2 = 5.76949 72214 60691 40550D0,
* C3 = 3.64784 83247 63204 60504D0,
* C4 = 1.27045 82524 52368 38258D0,
* C5 = 2.41780 72517 74506 11770D-1,
* C6 = 2.27238 44989 26918 45833D-2,
* C7 = 7.74545 01427 83414 07640D-4,
* D1 = 2.05319 16266 37758 82187D0,
* D2 = 1.67638 48301 83803 84940D0,
* D3 = 6.89767 33498 51000 04550D-1,
* D4 = 1.48103 97642 74800 74590D-1,
* D5 = 1.51986 66563 61645 71966D-2,
* D6 = 5.47593 80849 95344 94600D-4,
* D7 = 1.05075 00716 44416 84324D-9 )
* HASH SUM CD 49.33206 50330 16102 89036
*
* Coefficients for P near 0 or 1.
*
PARAMETER (
* E0 = 6.65790 46435 01103 77720D0,
* E1 = 5.46378 49111 64114 36990D0,
* E2 = 1.78482 65399 17291 33580D0,
* E3 = 2.96560 57182 85048 91230D-1,
* E4 = 2.65321 89526 57612 30930D-2,
* E5 = 1.24266 09473 88078 43860D-3,
* E6 = 2.71155 55687 43487 57815D-5,
* E7 = 2.01033 43992 92288 13265D-7,
* F1 = 5.99832 20655 58879 37690D-1,
* F2 = 1.36929 88092 27358 05310D-1,
* F3 = 1.48753 61290 85061 48525D-2,
* F4 = 7.86869 13114 56132 59100D-4,
* F5 = 1.84631 83175 10054 68180D-5,
* F6 = 1.42151 17583 16445 88870D-7,
* F7 = 2.04426 31033 89939 78564D-15 )
* HASH SUM EF 47.52583 31754 92896 71629
*
Q = ( P - HALF)
IF ( ABS(Q) .LE. SPLIT1 ) THEN ! Central range.
R = CONST1 - Q*Q
VAL = Q*( ( ( ((((A7*R + A6)*R + A5)*R + A4)*R + A3)
* *R + A2 )*R + A1 )*R + A0 )
* /( ( ( ((((B7*R + B6)*R + B5)*R + B4)*R + B3)
* *R + B2 )*R + B1 )*R + ONE)
ELSE ! near the endpoints
R = MIN( P, ONE - P )
IF (R .GT.ZERO) THEN ! ( 2.d0*R .GT. CFxCutOff) THEN ! R .GT.0.d0
R = SQRT( -LOG(R) )
IF ( R .LE. SPLIT2 ) THEN
R = R - CONST2
VAL = ( ( ( ((((C7*R + C6)*R + C5)*R + C4)*R + C3)
* *R + C2 )*R + C1 )*R + C0 )
* /( ( ( ((((D7*R + D6)*R + D5)*R + D4)*R + D3)
* *R + D2 )*R + D1 )*R + ONE )
ELSE
R = R - SPLIT2
VAL = ( ( ( ((((E7*R + E6)*R + E5)*R + E4)*R + E3)
* *R + E2 )*R + E1 )*R + E0 )
* /( ( ( ((((F7*R + F6)*R + F5)*R + F4)*R + F3)
* *R + F2 )*R + F1 )*R + ONE )
END IF
ELSE
VAL = 37.D0 !XMAX 9.d0
END IF
IF ( Q < ZERO ) VAL = - VAL
END IF
RETURN
END FUNCTION FIINV
FUNCTION FI2( Z ) RESULT (VALUE)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(in) :: Z
DOUBLE PRECISION :: VALUE
*
* Normal distribution probabilities accurate to 1.e-15.
* relative error less than 1e-8;
* Z = no. of standard deviations from the mean.
*
* Based upon algorithm 5666 for the error function, from:
* Hart, J.F. et al, 'Computer Approximations', Wiley 1968
*
* Programmer: Alan Miller
*
* Latest revision - 30 March 1986
*
DOUBLE PRECISION :: P0, P1, P2, P3, P4, P5, P6,
* Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7,XMAX,
* P, EXPNTL, CUTOFF, ROOTPI, ZABS, Z2
PARAMETER(
* P0 = 220.20 68679 12376 1D0,
* P1 = 221.21 35961 69931 1D0,
* P2 = 112.07 92914 97870 9D0,
* P3 = 33.912 86607 83830 0D0,
* P4 = 6.3739 62203 53165 0D0,
* P5 = 0.70038 30644 43688 1D0,
* P6 = 0.035262 49659 98910 9D0 )
PARAMETER(
* Q0 = 440.41 37358 24752 2D0,
* Q1 = 793.82 65125 19948 4D0,
* Q2 = 637.33 36333 78831 1D0,
* Q3 = 296.56 42487 79673 7D0,
* Q4 = 86.780 73220 29460 8D0,
* Q5 = 16.064 17757 92069 5D0,
* Q6 = 1.7556 67163 18264 2D0,
* Q7 = 0.088388 34764 83184 4D0 )
PARAMETER( ROOTPI = 2.5066 28274 63100 1D0 )
PARAMETER( CUTOFF = 7.0710 67811 86547 5D0 )
PARAMETER( XMAX = 8.25D0 )
*
ZABS = ABS(Z)
*
* |Z| > 37 (or XMAX)
*
IF ( Z .GT. XMAX .OR. ZABS .GT. 37) THEN
P = 0.d0
ELSE
*
* |Z| <= 37
*
Z2 = ZABS * ZABS
EXPNTL = EXP( -Z2 * 0.5D0 )
*
* |Z| < CUTOFF = 10/SQRT(2)
*
IF ( ZABS < CUTOFF ) THEN
P = EXPNTL*( (((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS
* + P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS
* + Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS
* + Q0 )
*
* |Z| >= CUTOFF.
*
ELSE
P = EXPNTL/( ZABS + 1.d0/( ZABS + 2.d0/( ZABS + 3.d0/( ZABS
* + 4.d0/( ZABS + 0.65D0 ) ) ) ) )/ROOTPI
END IF
END IF
IF ( Z .GT. 0.d0 ) P = 1.d0 - P
VALUE = P
RETURN
END FUNCTION FI2
FUNCTION FI( Z ) RESULT (VALUE)
USE ERFCOREMOD
IMPLICIT NONE
DOUBLE PRECISION, INTENT(in) :: Z
DOUBLE PRECISION :: VALUE
! Local variables
DOUBLE PRECISION, PARAMETER:: SQ2M1 = 0.70710678118655D0 ! 1/SQRT(2)
DOUBLE PRECISION, PARAMETER:: HALF = 0.5D0
VALUE = DERFC(-Z*SQ2M1)*HALF
RETURN
END FUNCTION FI
end module mvnProdCorrPrbMod