dbsgq8 Subroutine

private subroutine dbsgq8(fun, xt, bc, n, kk, id, a, b, inbv, err, ans, iflag, work)

DBSGQ8, a modification of gaus8, integrates the product of fun(x) by the id-th derivative of a spline dbvalu between limits a and b using an adaptive 8-point Legendre-Gauss algorithm.

See also

History

  • 800901 Jones, R. E., (SNLA)
  • 890531 Changed all specific intrinsics to generic. (WRB)
  • 890911 Removed unnecessary intrinsics. (WRB)
  • 891214 Prologue converted to Version 4.0 format. (BAB)
  • 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
  • 900326 Removed duplicate information from DESCRIPTION section. (WRB)
  • 900328 Added TYPE section. (WRB)
  • 910408 Updated the AUTHOR section. (WRB)
  • Jacob Williams, 9/6/2017 : refactored to modern Fortran. Some changes. Added higher precision coefficients.

Arguments

Type IntentOptional Attributes Name
procedure(b1fqad_func) :: fun

name of external function of one argument which multiplies dbvalu.
real(kind=wp), intent(in), dimension(:) :: xt

knot array for dbvalu
real(kind=wp), intent(in), dimension(n) :: bc

b-coefficient array for dbvalu
integer(kind=ip), intent(in) :: n

number of b-coefficients for dbvalu
integer(kind=ip), intent(in) :: kk

order of the spline, kk>=1
integer(kind=ip), intent(in) :: id

Order of the spline derivative, 0<=id<=kk-1
real(kind=wp), intent(in) :: a

lower limit of integral
real(kind=wp), intent(in) :: b

upper limit of integral (may be less than a)
integer(kind=ip), intent(inout) :: inbv

initialization parameter for dbvalu
real(kind=wp), intent(inout) :: err

IN: is a requested pseudorelative error tolerance. normally pick a value of abs(err)<1e-3. ans will normally have no more error than abs(err) times the integral of the absolute value of fun(x)*[[dbvalu]]().

OUT: will be an estimate of the absolute error in ans if the input value of err was negative. (err is unchanged if the input value of err was nonnegative.) the estimated error is solely for information to the user and should not be used as a correction to the computed integral.
real(kind=wp), intent(out) :: ans

computed value of integral
integer(kind=ip), intent(out) :: iflag

a status code:

  • 0: ans most likely meets requested error tolerance, or a=b.
  • 1101: a and b are too nearly equal to allow normal integration. ans is set to zero.
  • 1102: ans probably does not meet requested error tolerance.
real(kind=wp), intent(inout), dimension(:) :: work

work vector of length 3*k for dbvalu

Calls

proc~~dbsgq8~~CallsGraph proc~dbsgq8 bspline_sub_module::dbsgq8 proc~dbvalu bspline_sub_module::dbvalu proc~dbsgq8->proc~dbvalu proc~dintrv bspline_sub_module::dintrv proc~dbvalu->proc~dintrv proc~get_temp_x_for_extrap bspline_sub_module::get_temp_x_for_extrap proc~dintrv->proc~get_temp_x_for_extrap

Called by

proc~~dbsgq8~~CalledByGraph proc~dbsgq8 bspline_sub_module::dbsgq8 proc~dbfqad bspline_sub_module::dbfqad proc~dbfqad->proc~dbsgq8 proc~db1fqad bspline_sub_module::db1fqad proc~db1fqad->proc~dbfqad proc~fintegral_1d bspline_oo_module::bspline_1d%fintegral_1d proc~fintegral_1d->proc~db1fqad

Source Code

    subroutine dbsgq8(fun,xt,bc,n,kk,id,a,b,inbv,err,ans,iflag,work)

    implicit none

    procedure(b1fqad_func)              :: fun     !! name of external function of one
                                                   !! argument which multiplies [[dbvalu]].
    integer(ip),intent(in)              :: n       !! number of b-coefficients for [[dbvalu]]
    integer(ip),intent(in)              :: kk      !! order of the spline, `kk>=1`
    real(wp),dimension(:),intent(in)    :: xt      !! knot array for [[dbvalu]]
    real(wp),dimension(n),intent(in)    :: bc      !! b-coefficient array for [[dbvalu]]
    integer(ip),intent(in)              :: id      !! Order of the spline derivative, `0<=id<=kk-1`
    real(wp),intent(in)                 :: a       !! lower limit of integral
    real(wp),intent(in)                 :: b       !! upper limit of integral (may be less than `a`)
    integer(ip),intent(inout)           :: inbv    !! initialization parameter for [[dbvalu]]
    real(wp),intent(inout)              :: err     !! **IN:** is a requested pseudorelative error
                                                   !! tolerance.  normally pick a value of
                                                   !! `abs(err)<1e-3`.  `ans` will normally
                                                   !! have no more error than `abs(err)` times
                                                   !! the integral of the absolute value of
                                                   !! `fun(x)*[[dbvalu]]()`.
                                                   !!
                                                   !! **OUT:** will be an estimate of the absolute
                                                   !! error in ans if the input value of `err`
                                                   !! was negative.  (`err` is unchanged if
                                                   !! the input value of `err` was nonnegative.)
                                                   !! the estimated error is solely for information
                                                   !! to the user and should not be used as a
                                                   !! correction to the computed integral.
    real(wp),intent(out)                :: ans     !! computed value of integral
    integer(ip),intent(out)             :: iflag   !! a status code:
                                                   !!
                                                   !! * 0: `ans` most likely meets requested
                                                   !!   error tolerance, or `a=b`.
                                                   !! * 1101: `a` and `b` are too nearly equal
                                                   !!   to allow normal integration.
                                                   !!   `ans` is set to zero.
                                                   !! * 1102: `ans` probably does not meet
                                                   !!   requested error tolerance.
    real(wp),dimension(:),intent(inout) :: work    !! work vector of length `3*k` for [[dbvalu]]

    integer(ip) :: k,l,lmn,lmx,mxl,nbits,nib,nlmx
    real(wp) :: ae,anib,area,c,ce,ee,ef,eps,est,gl,glr,tol,vr,x
    integer(ip),dimension(60)  :: lr
    real(wp),dimension(60) :: aa,hh,vl,gr

    integer(ip),parameter  :: i1mach14 = digits(1.0_wp)            !! i1mach(14)
    real(wp),parameter     :: d1mach5  = log10(real(radix(x),wp))  !! d1mach(5)
    real(wp),parameter     :: ln2      = log(2.0_wp)               !! 0.69314718d0
    real(wp),parameter     :: sq2      = sqrt(2.0_wp)
    integer(ip),parameter  :: nlmn     = 1
    integer(ip),parameter  :: kmx      = 5000
    integer(ip),parameter  :: kml      = 6

    ! initialize
    inbv  = 1_ip
    iflag = 0_ip
    k     = i1mach14
    anib  = d1mach5*k/0.30102000_wp
    nbits = int(anib,ip)
    nlmx  = min((nbits*5_ip)/8_ip,60_ip)
    ans   = 0.0_wp
    ce    = 0.0_wp

    if ( a==b ) then
        if ( err<0.0_wp ) err = ce
    else
        lmx = nlmx
        lmn = nlmn
        if ( b/=0.0_wp ) then
            if ( sign(1.0_wp,b)*a>0.0_wp ) then
                c = abs(1.0_wp-a/b)
                if ( c<=0.1_wp ) then
                    if ( c<=0.0_wp ) then
                        if ( err<0.0_wp ) err = ce
                        return
                    else
                        anib = 0.5_wp - log(c)/ln2
                        nib = int(anib,ip)
                        lmx = min(nlmx,nbits-nib-7_ip)
                        if ( lmx<1_ip ) then
                            ! a and b are too nearly equal
                            ! to allow normal integration
                            iflag = 1101_ip
                            if ( err<0.0_wp ) err = ce
                            return
                        else
                            lmn = min(lmn,lmx)
                        end if
                    end if
                end if
            end if
        end if
        tol = max(abs(err),2.0_wp**(5-nbits))/2.0_wp
        if ( err==0.0_wp ) tol = sqrt(epsilon(1.0_wp))
        eps = tol
        hh(1_ip) = (b-a)/4.0_wp
        aa(1_ip) = a
        lr(1_ip) = 1_ip
        l = 1_ip
        call g8(aa(l)+2.0_wp*hh(l),2.0_wp*hh(l),est,iflag)
        if (iflag/=0_ip) return
        k = 8_ip
        area = abs(est)
        ef = 0.5_wp
        mxl = 0_ip
    end if

    do
        ! compute refined estimates, estimate the error, etc.
        call g8(aa(l)+hh(l),hh(l),gl,iflag)
        if (iflag/=0_ip) return
        call g8(aa(l)+3.0_wp*hh(l),hh(l),gr(l),iflag)
        if (iflag/=0_ip) return
        k = k + 16_ip
        area = area + (abs(gl)+abs(gr(l))-abs(est))
        glr = gl + gr(l)
        ee = abs(est-glr)*ef
        ae = max(eps*area,tol*abs(glr))
        if ( ee>ae ) then
            ! consider the left half of this level
            if ( k>kmx ) lmx = kml
            if ( l>=lmx ) then
                mxl = 1_ip
            else
                l = l + 1_ip
                eps = eps*0.5_wp
                ef = ef/sq2
                hh(l) = hh(l-1)*0.5_wp
                lr(l) = -1_ip
                aa(l) = aa(l-1_ip)
                est = gl
                cycle
            end if
        end if
        ce = ce + (est-glr)
        if ( lr(l)<=0_ip ) then
            ! proceed to right half at this level
            vl(l) = glr
        else
            ! return one level
            vr = glr
            do
                if ( l<=1_ip ) then
                    ! exit
                    ans = vr
                    if ( (mxl/=0_ip) .and. (abs(ce)>2.0_wp*tol*area) ) then
                        iflag = 1102_ip
                    end if
                    if ( err<0.0_wp ) err = ce
                    return
                else
                    l = l - 1_ip
                    eps = eps*2.0_wp
                    ef = ef*sq2
                    if ( lr(l)<=0 ) then
                        vl(l) = vl(l+1_ip) + vr
                        exit
                    else
                        vr = vl(l+1_ip) + vr
                    end if
                end if
            end do
        end if
        est = gr(l-1_ip)
        lr(l) = 1_ip
        aa(l) = aa(l) + 4.0_wp*hh(l)
    end do

    contains

        subroutine g8(x,h,res,iflag)

        !! 8-point formula.
        !!
        !!@note Replaced the original double precision abscissa and weight
        !!      coefficients with the higher precision versions from here:
        !!      http://pomax.github.io/bezierinfo/legendre-gauss.html
        !!      So, if `wp` is changed to say, `real128`, more precision
        !!      can be obtained. These coefficients have about 300 digits.

        implicit none

        real(wp),intent(in)     :: x
        real(wp),intent(in)     :: h
        real(wp),intent(out)    :: res
        integer(ip),intent(out) :: iflag

        real(wp),dimension(8) :: f
        real(wp),dimension(8) :: v

        ! abscissa and weight coefficients:
        real(wp),parameter :: x1 = &
        &0.1834346424956498049394761423601839806667578129129737823171884736992044&
        &742215421141160682237111233537452676587642867666089196012523876865683788&
        &569995160663568104475551617138501966385810764205532370882654749492812314&
        &961247764619363562770645716456613159405134052985058171969174306064445289&
        &638150514997832_wp
        real(wp),parameter :: x2 = &
        &0.5255324099163289858177390491892463490419642431203928577508570992724548&
        &207685612725239614001936319820619096829248252608507108793766638779939805&
        &395303668253631119018273032402360060717470006127901479587576756241288895&
        &336619643528330825624263470540184224603688817537938539658502113876953598&
        &879150514997832_wp
        real(wp),parameter :: x3 = &
        &0.7966664774136267395915539364758304368371717316159648320701702950392173&
        &056764730921471519272957259390191974534530973092653656494917010859602772&
        &562074621689676153935016290342325645582634205301545856060095727342603557&
        &415761265140428851957341933710803722783136113628137267630651413319993338&
        &002150514997832_wp
        real(wp),parameter :: x4 = &
        &0.9602898564975362316835608685694729904282352343014520382716397773724248&
        &977434192844394389592633122683104243928172941762102389581552171285479373&
        &642204909699700433982618326637346808781263553346927867359663480870597542&
        &547603929318533866568132868842613474896289232087639988952409772489387324&
        &25615051499783203_wp
        real(wp),parameter :: w1 = &
        &0.3626837833783619829651504492771956121941460398943305405248230675666867&
        &347239066773243660420848285095502587699262967065529258215569895173844995&
        &576007862076842778350382862546305771007553373269714714894268328780431822&
        &779077846722965535548199601402487767505928976560993309027632737537826127&
        &502150514997832_wp
        real(wp),parameter :: w2 = &
        &0.3137066458778872873379622019866013132603289990027349376902639450749562&
        &719421734969616980762339285560494275746410778086162472468322655616056890&
        &624276469758994622503118776562559463287222021520431626467794721603822601&
        &295276898652509723185157998353156062419751736972560423953923732838789657&
        &919150514997832_wp
        real(wp),parameter :: w3 = &
        &0.2223810344533744705443559944262408844301308700512495647259092892936168&
        &145704490408536531423771979278421592661012122181231114375798525722419381&
        &826674532090577908613289536840402789398648876004385697202157482063253247&
        &195590228631570651319965589733545440605952819880671616779621183704306688&
        &233150514997832_wp
        real(wp),parameter :: w4 = &
        &0.1012285362903762591525313543099621901153940910516849570590036980647401&
        &787634707848602827393040450065581543893314132667077154940308923487678731&
        &973041136073584690533208824050731976306575729205467961435779467552492328&
        &730055025992954089946676810510810729468366466585774650346143712142008566&
        &866150514997832_wp

        res = 0.0_wp

        v(1_ip) = x-x1*h
        v(2_ip) = x+x1*h
        v(3_ip) = x-x2*h
        v(4_ip) = x+x2*h
        v(5_ip) = x-x3*h
        v(6_ip) = x+x3*h
        v(7_ip) = x-x4*h
        v(8_ip) = x+x4*h

        call dbvalu(xt,bc,n,kk,id,v(1_ip),inbv,work,iflag,f(1_ip)); if (iflag/=0_ip) return
        call dbvalu(xt,bc,n,kk,id,v(2_ip),inbv,work,iflag,f(2_ip)); if (iflag/=0_ip) return
        call dbvalu(xt,bc,n,kk,id,v(3_ip),inbv,work,iflag,f(3_ip)); if (iflag/=0_ip) return
        call dbvalu(xt,bc,n,kk,id,v(4_ip),inbv,work,iflag,f(4_ip)); if (iflag/=0_ip) return
        call dbvalu(xt,bc,n,kk,id,v(5_ip),inbv,work,iflag,f(5_ip)); if (iflag/=0_ip) return
        call dbvalu(xt,bc,n,kk,id,v(6_ip),inbv,work,iflag,f(6_ip)); if (iflag/=0_ip) return
        call dbvalu(xt,bc,n,kk,id,v(7_ip),inbv,work,iflag,f(7_ip)); if (iflag/=0_ip) return
        call dbvalu(xt,bc,n,kk,id,v(8_ip),inbv,work,iflag,f(8_ip)); if (iflag/=0_ip) return

        res = h*((w1*(fun(v(1_ip))*f(1_ip) + fun(v(2_ip))*f(2_ip))  + &
                  w2*(fun(v(3_ip))*f(3_ip) + fun(v(4_ip))*f(4_ip))) + &
                 (w3*(fun(v(5_ip))*f(5_ip) + fun(v(6_ip))*f(6_ip))  + &
                  w4*(fun(v(7_ip))*f(7_ip) + fun(v(8_ip))*f(8_ip))))

        end subroutine g8

    end subroutine dbsgq8