dqagse Subroutine

    subroutine dqagse(f, a, b, Epsabs, Epsrel, Limit, Result, Abserr, Neval, &
                      Ier, Alist, Blist, Rlist, Elist, Iord, Last)
        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.
        integer, intent(in) :: Limit !! gives an upperbound on the number of subintervals
                                     !! in the partition of `(a,b)`
        real(wp), intent(out) :: Result !! approximation to the integral
        integer, intent(out) :: Neval !! number of integrand evaluations
        real(wp), intent(out) :: Abserr !! estimate of the modulus of the absolute error,
                                        !! which should equal or exceed `abs(i-result)`
        real(wp), intent(out) :: Alist(Limit) !! 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)
        real(wp), intent(out) :: Blist(Limit) !! 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)
        real(wp), intent(out) :: Elist(Limit) !! vector of dimension at least `limit`, the first
                                              !! `last` elements of which are the moduli of the
                                              !! absolute error estimates on the subintervals
        real(wp), intent(out) :: Rlist(Limit) !! vector of dimension at least `limit`, the first
                                              !! `last` elements of which are the integral
                                              !! approximations on the subintervals
        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
                                    !!   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 sub-
                                    !!   divisions 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 (e.g. singularity,
                                    !!   discontinuity within the interval) one
                                    !!   will probably gain from splitting up the
                                    !!   interval at this point and calling the
                                    !!   integrator on the subranges. if possible,
                                    !!   an appropriate special-purpose integrator
                                    !!   should be used, which is designed for
                                    !!   handling the type of difficulty involved.
                                    !! * ier = 2 the occurrence of roundoff error is
                                    !!   detected, which prevents the requested
                                    !!   tolerance from being achieved.
                                    !!   the error may be under-estimated.
                                    !! * ier = 3 extremely bad integrand behaviour
                                    !!   occurs at some points of the integration
                                    !!   interval.
                                    !! * ier = 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.
                                    !! * ier = 5 the integral is probably divergent, or
                                    !!   slowly convergent. it must be noted that
                                    !!   divergence can occur with any other value
                                    !!   of ier.
                                    !! * ier = 6 the input is invalid, because
                                    !!   `epsabs<=0` and
                                    !!   `epsrel<max(50*rel.mach.acc.,0.5e-28)`.
                                    !!   `result`, `abserr`, `neval`, `last`, `rlist(1)`,
                                    !!   `iord(1)` and `elist(1)` are set to zero.
                                    !!   `alist(1)` and `blist(1)` are set to a and b
                                    !!   respectively.
        integer, intent(out) :: Iord(Limit) !! 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<=(limit/2+2)`, and `k = limit+1-last`
                                            !! otherwise
        integer, intent(out) :: Last !! number of subintervals actually produced in the
                                     !! subdivision process

        real(wp) :: abseps, correc, defabs, dres, &
                    ertest, resabs, reseps, res3la(3)
        integer :: id, ierro, iroff1, iroff2, iroff3, &
                   jupbnd, k, ksgn, ktmin, nrmax
        real(wp) :: area12 !! `area1 + area2`
        real(wp) :: erro12 !! `error1 + error2`
        real(wp) :: area1, a1, b1, defab1, error1 !! variable for the left interval
        real(wp) :: area2, a2, b2, defab2, error2 !! variable for the right interval
        real(wp) :: rlist2(limexp + 2) !! array of dimension at least `limexp+2` containing
                                       !! the part of the epsilon table which is still
                                       !! needed for further computations.
        integer :: maxerr !! pointer to the interval with largest error estimate
        integer :: nres !! number of calls to the extrapolation routine
        integer :: numrl2 !! number of elements currently 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.
        real(wp) :: errmax !! elist(maxerr)
        real(wp) :: erlast !! error on the interval currently subdivided
                           !! (before that subdivision has taken place)
        real(wp) :: area !! sum of the integrals over the subintervals
        real(wp) :: errsum !! sum of the errors over the subintervals
        real(wp) :: errbnd !! requested accuracy `max(epsabs,epsrel*abs(result))`
        real(wp) :: small !! length of the smallest interval considered up
                          !! to now, multiplied by 1.5
        real(wp) :: erlarg !! sum of the errors over the intervals larger
                           !! than the smallest interval considered up to now
        logical :: 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`.
        logical :: noext !! logical variable denoting that extrapolation
                         !! is no longer allowed (true value)

        ! test on validity of parameters

        Ier = 0
        Neval = 0
        Last = 0
        Result = 0.0_wp
        Abserr = 0.0_wp
        Alist(1) = a
        Blist(1) = b
        Rlist(1) = 0.0_wp
        Elist(1) = 0.0_wp
        if (Epsabs <= 0.0_wp .and. Epsrel < max(50.0_wp*epmach, 0.5e-28_wp)) Ier = 6
        if (Ier == 6) return

        main : block

            ! first approximation to the integral

            ierro = 0
            call dqk21(f, a, b, Result, Abserr, defabs, resabs)

            ! test on accuracy.

            dres = abs(Result)
            errbnd = max(Epsabs, Epsrel*dres)
            Last = 1
            Rlist(1) = Result
            Elist(1) = Abserr
            Iord(1) = 1
            if (Abserr <= 100.0_wp*epmach*defabs .and. Abserr > errbnd) &
                Ier = 2
            if (Limit == 1) Ier = 1
            if (Ier /= 0 .or. (Abserr <= errbnd .and. Abserr /= resabs) .or. &
                Abserr == 0.0_wp) then
                Neval = 42*Last - 21
                return
            end if

            ! initialization

            rlist2(1) = Result
            errmax = Abserr
            maxerr = 1
            area = Result
            errsum = Abserr
            Abserr = oflow
            nrmax = 1
            nres = 0
            numrl2 = 2
            ktmin = 0
            extrap = .false.
            noext = .false.
            iroff1 = 0
            iroff2 = 0
            iroff3 = 0
            ksgn = -1
            if (dres >= (1.0_wp - 50.0_wp*epmach)*defabs) ksgn = 1

            ! main do-loop

            loop: do Last = 2, Limit

                ! bisect the subinterval with the nrmax-th largest error
                ! estimate.

                a1 = Alist(maxerr)
                b1 = 0.5_wp*(Alist(maxerr) + Blist(maxerr))
                a2 = b1
                b2 = Blist(maxerr)
                erlast = errmax
                call dqk21(f, a1, b1, area1, error1, resabs, defab1)
                call dqk21(f, a2, b2, area2, error2, resabs, defab2)

                ! improve previous approximations to integral
                ! and error and test for accuracy.

                area12 = area1 + area2
                erro12 = error1 + error2
                errsum = errsum + erro12 - errmax
                area = area + area12 - Rlist(maxerr)
                if (defab1 /= error1 .and. defab2 /= error2) then
                    if (abs(Rlist(maxerr) - area12) <= 0.1e-4_wp*abs(area12) &
                        .and. erro12 >= 0.99_wp*errmax) then
                        if (extrap) iroff2 = iroff2 + 1
                        if (.not. extrap) iroff1 = iroff1 + 1
                    end if
                    if (Last > 10 .and. erro12 > errmax) iroff3 = iroff3 + 1
                end if
                Rlist(maxerr) = area1
                Rlist(Last) = area2
                errbnd = max(Epsabs, Epsrel*abs(area))

                ! test for roundoff error and eventually set error flag.

                if (iroff1 + iroff2 >= 10 .or. iroff3 >= 20) Ier = 2
                if (iroff2 >= 5) ierro = 3

                ! set error flag in the case that the number of subintervals
                ! equals limit.

                if (Last == Limit) Ier = 1

                ! set error flag in the case of bad integrand behaviour
                ! at a point of the integration range.

                if (max(abs(a1), abs(b2)) <= (1.0_wp + 100.0_wp*epmach) &
                    *(abs(a2) + 1000.0_wp*uflow)) Ier = 4

                ! append the newly-created intervals to the list.

                if (error2 > error1) then
                    Alist(maxerr) = a2
                    Alist(Last) = a1
                    Blist(Last) = b1
                    Rlist(maxerr) = area2
                    Rlist(Last) = area1
                    Elist(maxerr) = error2
                    Elist(Last) = error1
                else
                    Alist(Last) = a2
                    Blist(maxerr) = b1
                    Blist(Last) = b2
                    Elist(maxerr) = error1
                    Elist(Last) = error2
                end if

                ! 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).

                call dqpsrt(Limit, Last, maxerr, errmax, Elist, Iord, nrmax)
                ! ***jump out of do-loop
                if (errsum <= errbnd) exit main
                ! ***jump out of do-loop
                if (Ier /= 0) exit loop
                if (Last == 2) then
                    small = abs(b - a)*0.375_wp
                    erlarg = errsum
                    ertest = errbnd
                    rlist2(2) = area
                elseif (.not. (noext)) then
                    erlarg = erlarg - erlast
                    if (abs(b1 - a1) > small) erlarg = erlarg + erro12
                    if (.not. (extrap)) then
                        ! test whether the interval to be bisected next is the
                        ! smallest interval.
                        if (abs(Blist(maxerr) - Alist(maxerr)) > small) cycle loop
                        extrap = .true.
                        nrmax = 2
                    end if
                    if (ierro /= 3 .and. erlarg > ertest) then
                        ! 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 > (2 + Limit/2)) jupbnd = Limit + 3 - Last
                        do k = id, jupbnd
                            maxerr = Iord(nrmax)
                            errmax = Elist(maxerr)
                            ! ***jump out of do-loop
                            if (abs(Blist(maxerr) - Alist(maxerr)) > small) cycle loop
                            nrmax = nrmax + 1
                        end do
                    end if

                    ! perform extrapolation.

                    numrl2 = numrl2 + 1
                    rlist2(numrl2) = area
                    call dqelg(numrl2, rlist2, reseps, abseps, res3la, nres)
                    ktmin = ktmin + 1
                    if (ktmin > 5 .and. Abserr < 0.1e-02_wp*errsum) Ier = 5
                    if (abseps < Abserr) then
                        ktmin = 0
                        Abserr = abseps
                        Result = reseps
                        correc = erlarg
                        ertest = max(Epsabs, Epsrel*abs(reseps))
                        ! ***jump out of do-loop
                        if (Abserr <= ertest) exit loop
                    end if

                    ! prepare bisection of the smallest interval.

                    if (numrl2 == 1) noext = .true.
                    if (Ier == 5) exit loop
                    maxerr = Iord(1)
                    errmax = Elist(maxerr)
                    nrmax = 1
                    extrap = .false.
                    small = small*0.5_wp
                    erlarg = errsum
                end if
            end do loop

            ! set final result and error estimate.

            if (Abserr /= oflow) then
                if (Ier + ierro /= 0) then
                    if (ierro == 3) Abserr = Abserr + correc
                    if (Ier == 0) Ier = 3
                    if (Result == 0.0_wp .or. area == 0.0_wp) then
                        if (Abserr > errsum) exit main
                        if (area == 0.0_wp) then
                            if (Ier > 2) Ier = Ier - 1
                            Neval = 42*Last - 21
                            return
                        end if
                    elseif (Abserr/abs(Result) > errsum/abs(area)) then
                        exit main
                    end if
                end if

                ! test on divergence.

                if (ksgn /= (-1) .or. max(abs(Result), abs(area)) &
                    > defabs*0.01_wp) then
                    if (0.01_wp > (Result/area) .or. (Result/area) &
                        > 100.0_wp .or. errsum > abs(area)) Ier = 6
                end if
                if (Ier > 2) Ier = Ier - 1
                Neval = 42*Last - 21
                return
            end if

        end block main

        ! compute global integral sum.

        Result = sum(Rlist(1:Last))
        Abserr = errsum
        if (Ier > 2) Ier = Ier - 1
        Neval = 42*Last - 21

    end subroutine dqagse