dqng – quadpack

public subroutine dqng(f, a, b, Epsabs, Epsrel, Result, Abserr, Neval, Ier)

1D non-adaptive automatic integrator
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)<=max(epsabs,epsrel*abs(i)).
History

QUADPACK: date written 800101, revision date 810101 (yymmdd), kahaner,david,nbs - modified (2/82)
Arguments

Type	Intent		Name
procedure(func)		::	f	function subprogram defining the integrand function `f(x)`.
real(kind=wp),	intent(in)	::	a	lower limit of integration
real(kind=wp),	intent(in)	::	b	upper limit of integration
real(kind=wp),	intent(in)	::	Epsabs	absolute accuracy requested
real(kind=wp),	intent(in)	::	Epsrel	relative accuracy requested if `epsabs<=0` and `epsrel<max(50*rel.mach.acc.,0.5e-28)`, the routine will end with ier = 6.
real(kind=wp),	intent(out)	::	Result	approximation to the integral i result is obtained by applying the 21-point gauss-kronrod rule (res21) obtained by optimal addition of abscissae to the 10-point gauss rule (res10), or by applying the 43-point rule (res43) obtained by optimal addition of abscissae to the 21-point gauss-kronrod rule, or by applying the 87-point rule (res87) obtained by optimal addition of abscissae to the 43-point rule.
real(kind=wp),	intent(out)	::	Abserr	estimate of the modulus of the absolute error, which should equal or exceed `abs(i-result)`
integer,	intent(out)	::	Neval	number of integrand evaluations
integer,	intent(out)	::	Ier	ier = 0 normal and reliable termination of the routine. it is assumed that the requested accuracy has been achieved. ier>0 abnormal termination of the routine. it is assumed that the requested accuracy has not been achieved. error messages: ier = 1 the maximum number of steps has been executed. the integral is probably too difficult to be calculated by dqng. ier = 6 the input is invalid, because `epsabs<=0` and `epsrel<max(50*rel.mach.acc.,0.5e-28)`. `result`, `abserr` and `neval` are set to zero.
Source Code

    subroutine dqng(f, a, b, Epsabs, Epsrel, Result, Abserr, Neval, Ier)
        implicit none

        procedure(func) :: f !! function subprogram defining the integrand function `f(x)`.
        real(wp), intent(in) :: a !! lower limit of integration
        real(wp), intent(in) :: b !! upper limit of integration
        real(wp), intent(in) :: Epsabs !! absolute accuracy requested
        real(wp), intent(in) :: Epsrel !! relative accuracy requested
                                      !! if `epsabs<=0`
                                      !! and `epsrel<max(50*rel.mach.acc.,0.5e-28)`,
                                      !! the routine will end with ier = 6.
        real(wp), intent(out) :: Result !! approximation to the integral i
                                        !! result is obtained by applying the 21-point
                                        !! gauss-kronrod rule (res21) obtained by optimal
                                        !! addition of abscissae to the 10-point gauss rule
                                        !! (res10), or by applying the 43-point rule (res43)
                                        !! obtained by optimal addition of abscissae to the
                                        !! 21-point gauss-kronrod rule, or by applying the
                                        !! 87-point rule (res87) obtained by optimal addition
                                        !! of abscissae to the 43-point rule.
        real(wp), intent(out) :: Abserr !! estimate of the modulus of the absolute error,
                                        !! which should equal or exceed `abs(i-result)`
        integer, intent(out) :: Neval !! number of integrand evaluations
        integer, intent(out) :: Ier !! * ier = 0 normal and reliable termination of the
                                    !!   routine. it is assumed that the requested
                                    !!   accuracy has been achieved.
                                    !! * ier>0 abnormal termination of the routine. it is
                                    !!   assumed that the requested accuracy has
                                    !!   not been achieved.
                                    !!
                                    !! error messages:
                                    !!
                                    !! * ier = 1 the maximum number of steps has been
                                    !!   executed. the integral is probably too
                                    !!   difficult to be calculated by dqng.
                                    !! * ier = 6 the input is invalid, because
                                    !!   `epsabs<=0` and
                                    !!   `epsrel<max(50*rel.mach.acc.,0.5e-28)`.
                                    !!   `result`, `abserr` and `neval` are set to zero.

        real(wp) :: dhlgth, fval1, fval2, fv1(5), fv2(5), fv3(5), fv4(5), reskh
        integer :: ipx, k, l
        real(wp) :: centr !! mid point of the integration interval
        real(wp) :: hlgth !! half-length of the integration interval
        real(wp) :: fcentr !! function value at mid point
        real(wp) :: absc !! abscissa
        real(wp) :: fval !! function value
        real(wp) :: savfun(21) !! array of function values which have already been computed
        real(wp) :: res10 !! 10-point gauss result
        real(wp) :: res21 !! 21-point kronrod result
        real(wp) :: res43 !! 43-point result
        real(wp) :: res87 !! 87-point result
        real(wp) :: resabs !! approximation to the integral of `abs(f)`
        real(wp) :: resasc !! approximation to the integral of `abs(f-i/(b-a))`

        ! the following data statements contain the
        ! abscissae and weights of the integration rules used.

        real(wp), dimension(5), parameter :: x1 = [ &
                                             9.73906528517171720077964012084452053428e-1_wp, &
                                             8.65063366688984510732096688423493048528e-1_wp, &
                                             6.79409568299024406234327365114873575769e-1_wp, &
                                             4.33395394129247190799265943165784162200e-1_wp, &
                                             1.48874338981631210884826001129719984618e-1_wp] !! abscissae common to the 10-, 21-, 43- and 87-point rule

        real(wp), dimension(5), parameter :: w10 = [ &
                                             6.66713443086881375935688098933317928579e-2_wp, &
                                             1.49451349150580593145776339657697332403e-1_wp, &
                                             2.19086362515982043995534934228163192459e-1_wp, &
                                             2.69266719309996355091226921569469352860e-1_wp, &
                                             2.95524224714752870173892994651338329421e-1_wp] !! weights of the 10-point formula

        real(wp), dimension(5), parameter :: x2 = [ &
                                             9.95657163025808080735527280689002847921e-1_wp, &
                                             9.30157491355708226001207180059508346225e-1_wp, &
                                             7.80817726586416897063717578345042377163e-1_wp, &
                                             5.62757134668604683339000099272694140843e-1_wp, &
                                             2.94392862701460198131126603103865566163e-1_wp] !! abscissae common to the 21-, 43- and 87-point rule

        real(wp), dimension(5), parameter :: w21a = [ &
                                             3.25581623079647274788189724593897606174e-2_wp, &
                                             7.50396748109199527670431409161900093952e-2_wp, &
                                             1.09387158802297641899210590325804960272e-1_wp, &
                                             1.34709217311473325928054001771706832761e-1_wp, &
                                             1.47739104901338491374841515972068045524e-1_wp] !! weights of the 21-point formula for abscissae x1

        real(wp), dimension(6), parameter :: w21b = [ &
                                             1.16946388673718742780643960621920483962e-2_wp, &
                                             5.47558965743519960313813002445801763737e-2_wp, &
                                             9.31254545836976055350654650833663443900e-2_wp, &
                                             1.23491976262065851077958109831074159512e-1_wp, &
                                             1.42775938577060080797094273138717060886e-1_wp, &
                                             1.49445554002916905664936468389821203745e-1_wp] !! weights of the 21-point formula for abscissae x2

        ! 43 and 87 coefficients are computed via the algorithm in the quadpack
        ! manual, section 2.2.2.
        !TODO: They need to be regenerated with the same precision as the others.
        real(wp), dimension(11), parameter :: x3 = [ &
                                              0.999333360901932081394099323919911_wp, &
                                              0.987433402908088869795961478381209_wp, &
                                              0.954807934814266299257919200290473_wp, &
                                              0.900148695748328293625099494069092_wp, &
                                              0.825198314983114150847066732588520_wp, &
                                              0.732148388989304982612354848755461_wp, &
                                              0.622847970537725238641159120344323_wp, &
                                              0.499479574071056499952214885499755_wp, &
                                              0.364901661346580768043989548502644_wp, &
                                              0.222254919776601296498260928066212_wp, &
                                              0.074650617461383322043914435796506_wp] !! abscissae common to the 43- and 87-point rule

        real(wp), dimension(10), parameter :: w43a = [ &
                                              0.016296734289666564924281974617663_wp, &
                                              0.037522876120869501461613795898115_wp, &
                                              0.054694902058255442147212685465005_wp, &
                                              0.067355414609478086075553166302174_wp, &
                                              0.073870199632393953432140695251367_wp, &
                                              0.005768556059769796184184327908655_wp, &
                                              0.027371890593248842081276069289151_wp, &
                                              0.046560826910428830743339154433824_wp, &
                                              0.061744995201442564496240336030883_wp, &
                                              0.071387267268693397768559114425516_wp] !! weights of the 43-point formula for abscissae x1, x3

        real(wp), dimension(12), parameter :: w43b = [ &
                                              0.001844477640212414100389106552965_wp, &
                                              0.010798689585891651740465406741293_wp, &
                                              0.021895363867795428102523123075149_wp, &
                                              0.032597463975345689443882222526137_wp, &
                                              0.042163137935191811847627924327955_wp, &
                                              0.050741939600184577780189020092084_wp, &
                                              0.058379395542619248375475369330206_wp, &
                                              0.064746404951445885544689259517511_wp, &
                                              0.069566197912356484528633315038405_wp, &
                                              0.072824441471833208150939535192842_wp, &
                                              0.074507751014175118273571813842889_wp, &
                                              0.074722147517403005594425168280423_wp] !! weights of the 43-point formula for abscissae x3

        real(wp), dimension(22), parameter :: x4 = [ &
                                              0.999902977262729234490529830591582_wp, &
                                              0.997989895986678745427496322365960_wp, &
                                              0.992175497860687222808523352251425_wp, &
                                              0.981358163572712773571916941623894_wp, &
                                              0.965057623858384619128284110607926_wp, &
                                              0.943167613133670596816416634507426_wp, &
                                              0.915806414685507209591826430720050_wp, &
                                              0.883221657771316501372117548744163_wp, &
                                              0.845710748462415666605902011504855_wp, &
                                              0.803557658035230982788739474980964_wp, &
                                              0.757005730685495558328942793432020_wp, &
                                              0.706273209787321819824094274740840_wp, &
                                              0.651589466501177922534422205016736_wp, &
                                              0.593223374057961088875273770349144_wp, &
                                              0.531493605970831932285268948562671_wp, &
                                              0.466763623042022844871966781659270_wp, &
                                              0.399424847859218804732101665817923_wp, &
                                              0.329874877106188288265053371824597_wp, &
                                              0.258503559202161551802280975429025_wp, &
                                              0.185695396568346652015917141167606_wp, &
                                              0.111842213179907468172398359241362_wp, &
                                              0.037352123394619870814998165437704_wp] !! abscissae of the 87-point rule

        real(wp), dimension(21), parameter :: w87a = [ &
                                              0.008148377384149172900002878448190_wp, &
                                              0.018761438201562822243935059003794_wp, &
                                              0.027347451050052286161582829741283_wp, &
                                              0.033677707311637930046581056957588_wp, &
                                              0.036935099820427907614589586742499_wp, &
                                              0.002884872430211530501334156248695_wp, &
                                              0.013685946022712701888950035273128_wp, &
                                              0.023280413502888311123409291030404_wp, &
                                              0.030872497611713358675466394126442_wp, &
                                              0.035693633639418770719351355457044_wp, &
                                              0.000915283345202241360843392549948_wp, &
                                              0.005399280219300471367738743391053_wp, &
                                              0.010947679601118931134327826856808_wp, &
                                              0.016298731696787335262665703223280_wp, &
                                              0.021081568889203835112433060188190_wp, &
                                              0.025370969769253827243467999831710_wp, &
                                              0.029189697756475752501446154084920_wp, &
                                              0.032373202467202789685788194889595_wp, &
                                              0.034783098950365142750781997949596_wp, &
                                              0.036412220731351787562801163687577_wp, &
                                              0.037253875503047708539592001191226_wp] !! weights of the 87-point formula for abscissae x1, x2, x3

        real(wp), dimension(23), parameter :: w87b = [ &
                                              0.000274145563762072350016527092881_wp, &
                                              0.001807124155057942948341311753254_wp, &
                                              0.004096869282759164864458070683480_wp, &
                                              0.006758290051847378699816577897424_wp, &
                                              0.009549957672201646536053581325377_wp, &
                                              0.012329447652244853694626639963780_wp, &
                                              0.015010447346388952376697286041943_wp, &
                                              0.017548967986243191099665352925900_wp, &
                                              0.019938037786440888202278192730714_wp, &
                                              0.022194935961012286796332102959499_wp, &
                                              0.024339147126000805470360647041454_wp, &
                                              0.026374505414839207241503786552615_wp, &
                                              0.028286910788771200659968002987960_wp, &
                                              0.030052581128092695322521110347341_wp, &
                                              0.031646751371439929404586051078883_wp, &
                                              0.033050413419978503290785944862689_wp, &
                                              0.034255099704226061787082821046821_wp, &
                                              0.035262412660156681033782717998428_wp, &
                                              0.036076989622888701185500318003895_wp, &
                                              0.036698604498456094498018047441094_wp, &
                                              0.037120549269832576114119958413599_wp, &
                                              0.037334228751935040321235449094698_wp, &
                                              0.037361073762679023410321241766599_wp] !! weights of the 87-point formula for abscissae x4

        ! test on validity of parameters

        Result = 0.0_wp
        Abserr = 0.0_wp
        Neval = 0
        Ier = 6
        if (Epsabs > 0.0_wp .or. Epsrel >= max(50.0_wp*epmach, 0.5e-28_wp)) &
            then
            hlgth = 0.5_wp*(b - a)
            dhlgth = abs(hlgth)
            centr = 0.5_wp*(b + a)
            fcentr = f(centr)
            Neval = 21
            Ier = 1

            ! compute the integral using the 10- and 21-point formula.

            do l = 1, 3
                select case (l)
                case (2)

                    ! compute the integral using the 43-point formula.

                    res43 = w43b(12)*fcentr
                    Neval = 43
                    do k = 1, 10
                        res43 = res43 + savfun(k)*w43a(k)
                    end do
                    do k = 1, 11
                        ipx = ipx + 1
                        absc = hlgth*x3(k)
                        fval = f(absc + centr) + f(centr - absc)
                        res43 = res43 + fval*w43b(k)
                        savfun(ipx) = fval
                    end do

                    ! test for convergence.

                    Result = res43*hlgth
                    Abserr = abs((res43 - res21)*hlgth)
                case (3)

                    ! compute the integral using the 87-point formula.

                    res87 = w87b(23)*fcentr
                    Neval = 87
                    do k = 1, 21
                        res87 = res87 + savfun(k)*w87a(k)
                    end do
                    do k = 1, 22
                        absc = hlgth*x4(k)
                        res87 = res87 + w87b(k)*(f(absc + centr) + f(centr - absc))
                    end do
                    Result = res87*hlgth
                    Abserr = abs((res87 - res43)*hlgth)
                case default
                    res10 = 0.0_wp
                    res21 = w21b(6)*fcentr
                    resabs = w21b(6)*abs(fcentr)
                    do k = 1, 5
                        absc = hlgth*x1(k)
                        fval1 = f(centr + absc)
                        fval2 = f(centr - absc)
                        fval = fval1 + fval2
                        res10 = res10 + w10(k)*fval
                        res21 = res21 + w21a(k)*fval
                        resabs = resabs + w21a(k)*(abs(fval1) + abs(fval2))
                        savfun(k) = fval
                        fv1(k) = fval1
                        fv2(k) = fval2
                    end do
                    ipx = 5
                    do k = 1, 5
                        ipx = ipx + 1
                        absc = hlgth*x2(k)
                        fval1 = f(centr + absc)
                        fval2 = f(centr - absc)
                        fval = fval1 + fval2
                        res21 = res21 + w21b(k)*fval
                        resabs = resabs + w21b(k)*(abs(fval1) + abs(fval2))
                        savfun(ipx) = fval
                        fv3(k) = fval1
                        fv4(k) = fval2
                    end do

                    ! test for convergence.

                    Result = res21*hlgth
                    resabs = resabs*dhlgth
                    reskh = 0.5_wp*res21
                    resasc = w21b(6)*abs(fcentr - reskh)
                    do k = 1, 5
                        resasc = resasc + w21a(k) &
                                 *(abs(fv1(k) - reskh) + abs(fv2(k) - reskh)) &
                                 + w21b(k) &
                                 *(abs(fv3(k) - reskh) + abs(fv4(k) - reskh))
                    end do
                    Abserr = abs((res21 - res10)*hlgth)
                    resasc = resasc*dhlgth
                end select
                if (resasc /= 0.0_wp .and. Abserr /= 0.0_wp) &
                    Abserr = resasc*min(1.0_wp, (200.0_wp*Abserr/resasc)**1.5_wp)
                if (resabs > uflow/(50.0_wp*epmach)) &
                    Abserr = max((epmach*50.0_wp)*resabs, Abserr)
                if (Abserr <= max(Epsabs, Epsrel*abs(Result))) Ier = 0
                ! ***jump out of do-loop
                if (Ier == 0) return
            end do
        end if
        call xerror('abnormal return from dqng ', Ier, 0)

    end subroutine dqng
dqng Subroutine

Contents

Source Code

public subroutine dqng(f, a, b, Epsabs, Epsrel, Result, Abserr, Neval, Ier)

History

Arguments

Source Code
