dbsqad Subroutine

private pure subroutine dbsqad(t, bcoef, n, k, x1, x2, bquad, work, iflag)

DBSQAD computes the integral on (x1,x2) of a k-th order b-spline using the b-representation (t,bcoef,n,k). orders k as high as 20 are permitted by applying a 2, 6, or 10 point gauss formula on subintervals of (x1,x2) which are formed by included (distinct) knots.

If orders k greater than 20 are needed, use dbfqad with f(x) = 1.

Note

  • The maximum number of significant digits obtainable in DBSQAD is the smaller of ~300 and the number of digits carried in real(wp) arithmetic.

References

  • D. E. Amos, "Quadrature subroutines for splines and B-splines", Report SAND79-1825, Sandia Laboratories, December 1979.

History

  • Author: Amos, D. E., (SNLA)
  • 800901 DATE WRITTEN
  • 890531 Changed all specific intrinsics to generic. (WRB)
  • 890531 REVISION DATE from Version 3.2
  • 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)
  • 920501 Reformatted the REFERENCES section. (WRB)
  • Jacob Williams, 9/6/2017 : refactored to modern Fortran. Added higher precision coefficients.

Note

Extrapolation is not enabled for this routine.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in), dimension(:) :: t

knot array of length n+k
real(kind=wp), intent(in), dimension(:) :: bcoef

b-spline coefficient array of length n
integer(kind=ip), intent(in) :: n

length of coefficient array
integer(kind=ip), intent(in) :: k

order of b-spline, 1 <= k <= 20
real(kind=wp), intent(in) :: x1

end point of quadrature interval in t(k) <= x <= t(n+1)
real(kind=wp), intent(in) :: x2

end point of quadrature interval in t(k) <= x <= t(n+1)
real(kind=wp), intent(out) :: bquad

integral of the b-spline over (x1,x2)
real(kind=wp), intent(inout), dimension(:) :: work

work vector of length 3*k
integer(kind=ip), intent(out) :: iflag

status flag:

  • 0: no errors
  • 901: k does not satisfy 1<=k<=20
  • 902: n does not satisfy n>=k
  • 903: x1 or x2 or both do not satisfy t(k)<=x<=t(n+1)

Calls

proc~~dbsqad~~CallsGraph proc~dbsqad bspline_sub_module::dbsqad proc~dbvalu bspline_sub_module::dbvalu proc~dbsqad->proc~dbvalu proc~dintrv bspline_sub_module::dintrv proc~dbsqad->proc~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~~dbsqad~~CalledByGraph proc~dbsqad bspline_sub_module::dbsqad proc~db1sqad bspline_sub_module::db1sqad proc~db1sqad->proc~dbsqad proc~integral_1d bspline_oo_module::bspline_1d%integral_1d proc~integral_1d->proc~db1sqad

Source Code

    pure subroutine dbsqad(t,bcoef,n,k,x1,x2,bquad,work,iflag)

    implicit none

    real(wp),dimension(:),intent(in)    :: t       !! knot array of length `n+k`
    real(wp),dimension(:),intent(in)    :: bcoef   !! b-spline coefficient array of length `n`
    integer(ip),intent(in)              :: n       !! length of coefficient array
    integer(ip),intent(in)              :: k       !! order of b-spline, `1 <= k <= 20`
    real(wp),intent(in)                 :: x1      !! end point of quadrature interval
                                                   !! in `t(k) <= x <= t(n+1)`
    real(wp),intent(in)                 :: x2      !! end point of quadrature interval
                                                   !! in `t(k) <= x <= t(n+1)`
    real(wp),intent(out)                :: bquad   !! integral of the b-spline over (`x1`,`x2`)
    real(wp),dimension(:),intent(inout) :: work    !! work vector of length `3*k`
    integer(ip),intent(out)             :: iflag   !! status flag:
                                                   !!
                                                   !! * 0: no errors
                                                   !! * 901: `k` does not satisfy `1<=k<=20`
                                                   !! * 902: `n` does not satisfy `n>=k`
                                                   !! * 903: `x1` or `x2` or both do
                                                   !!   not satisfy `t(k)<=x<=t(n+1)`

    integer(ip) :: i,il1,il2,ilo,inbv,jf,left,m,mf,mflag,npk,np1
    real(wp) :: a,aa,b,bb,bma,bpa,c1,gx,q,ta,tb,y1,y2
    real(wp),dimension(5) :: s  !! sum

    real(wp),dimension(9),parameter :: gpts = [ &
        &0.577350269189625764509148780501957455647601751270126876018602326483977&
        &67230293334569371539558574952522520871380513556767665664836499965082627&
        &05518373647912161760310773007685273559916067003615583077550051041144223&
        &01107628883557418222973945990409015710553455953862673016662179126619796&
        &4892168_wp,&
        &0.238619186083196908630501721680711935418610630140021350181395164574274&
        &93427563984224922442725734913160907222309701068720295545303507720513526&
        &28872175189982985139866216812636229030578298770859440976999298617585739&
        &46921613621659222233462641640013936777894532787145324672151888999339900&
        &0945406150514997832_wp,&
        &0.661209386466264513661399595019905347006448564395170070814526705852183&
        &49660714310094428640374646145642988837163927514667955734677222538043817&
        &23198010093367423918538864300079016299442625145884902455718821970386303&
        &22362011735232135702218793618906974301231555871064213101639896769013566&
        &1651261150514997832_wp,&
        &0.932469514203152027812301554493994609134765737712289824872549616526613&
        &50084420019627628873992192598504786367972657283410658797137951163840419&
        &21786180750210169211578452038930846310372961174632524612619760497437974&
        &07422632089671621172178385230505104744277222209386367655366917903888025&
        &2326771150514997832_wp,&
        &0.148874338981631210884826001129719984617564859420691695707989253515903&
        &61735566852137117762979946369123003116080525533882610289018186437654023&
        &16761969968090913050737827720371059070942475859422743249837177174247346&
        &21691485290294292900319346665908243383809435507599683357023000500383728&
        &0634351_wp,&
        &0.433395394129247190799265943165784162200071837656246496502701513143766&
        &98907770350122510275795011772122368293504099893794727422475772324920512&
        &67741032822086200952319270933462032011328320387691584063411149801129823&
        &14148878744320432476641442157678880770848387945248811854979703928792696&
        &4254222_wp,&
        &0.679409568299024406234327365114873575769294711834809467664817188952558&
        &57539507492461507857357048037949983390204739931506083674084257663009076&
        &82741718202923543197852846977409718369143712013552962837733153108679126&
        &93254495485472934132472721168027426848661712101171203022718105101071880&
        &4444161_wp,&
        &0.865063366688984510732096688423493048527543014965330452521959731845374&
        &75513805556135679072894604577069440463108641176516867830016149345356373&
        &92729396890950011571349689893051612072435760480900979725923317923795535&
        &73929059587977695683242770223694276591148364371481692378170157259728913&
        &9322313_wp,&
        &0.973906528517171720077964012084452053428269946692382119231212066696595&
        &20323463615962572356495626855625823304251877421121502216860143447777992&
        &05409587259942436704413695764881258799146633143510758737119877875210567&
        &06745243536871368303386090938831164665358170712568697066873725922944928&
        &4383797_wp]

    real(wp),dimension(9),parameter :: gwts = [ &
        &1.0_wp,&
        &0.467913934572691047389870343989550994811655605769210535311625319963914&
        &20162039812703111009258479198230476626878975479710092836255417350295459&
        &35635592733866593364825926382559018030281273563502536241704619318259000&
        &99756987095900533474080074634376824431808173206369174103416261765346292&
        &7888917150514997832_wp,&
        &0.360761573048138607569833513837716111661521892746745482289739240237140&
        &03783726171832096220198881934794311720914037079858987989027836432107077&
        &67872114085818922114502722525757771126000732368828591631602895111800517&
        &40813685547074482472486101183259931449817216402425586777526768199930950&
        &3106873150514997832_wp,&
        &0.171324492379170345040296142172732893526822501484043982398635439798945&
        &76054234015464792770542638866975211652206987440430919174716746217597462&
        &96492293180314484520671351091683210843717994067668872126692485569940481&
        &59429327357024984053433824182363244118374610391205239119044219703570297&
        &7497812150514997832_wp,&
        &0.295524224714752870173892994651338329421046717026853601354308029755995&
        &93821715232927035659579375421672271716440125255838681849078955200582600&
        &19363424941869666095627186488841680432313050615358674090830512706638652&
        &87483901746874726597515954450775158914556548308329986393605934912382356&
        &670244_wp,&
        &0.269266719309996355091226921569469352859759938460883795800563276242153&
        &43231917927676422663670925276075559581145036869830869292346938114524155&
        &64658846634423711656014432259960141729044528030344411297902977067142537&
        &53480628460839927657500691168674984281408628886853320804215041950888191&
        &6391898_wp,&
        &0.219086362515982043995534934228163192458771870522677089880956543635199&
        &91065295128124268399317720219278659121687281288763476662690806694756883&
        &09211843316656677105269915322077536772652826671027878246851010208832173&
        &32006427348325475625066841588534942071161341022729156547776892831330068&
        &8702802_wp,&
        &0.149451349150580593145776339657697332402556639669427367835477268753238&
        &65472663001094594726463473195191400575256104543633823445170674549760147&
        &13716011937109528798134828865118770953566439639333773939909201690204649&
        &08381561877915752257830034342778536175692764212879241228297015017259084&
        &2897331_wp,&
        &0.066671344308688137593568809893331792857864834320158145128694881613412&
        &06408408710177678550968505887782109005471452041933148750712625440376213&
        &93049873169940416344953637064001870112423155043935262424506298327181987&
        &18647480566044117862086478449236378557180717569208295026105115288152794&
        &421677_wp]

    iflag = 0_ip
    bquad = 0.0_wp

    if ( k<1_ip .or. k>20_ip ) then

        iflag = 901_ip ! error return

    else if ( n<k ) then

        iflag = 902_ip ! error return

    else

        aa = min(x1,x2)
        bb = max(x1,x2)
        if ( aa>=t(k) ) then
            np1 = n + 1_ip
            if ( bb<=t(np1) ) then
            if ( aa==bb ) return
            npk = n + k
            ! selection of 2, 6, or 10 point gauss formula
            jf = 0_ip
            mf = 1_ip
            if ( k>4_ip ) then
                jf = 1_ip
                mf = 3_ip
                if ( k>12_ip ) then
                    jf = 4_ip
                    mf = 5_ip
                end if
            end if
            do i = 1_ip , mf
                s(i) = 0.0_wp
            end do
            ilo = 1_ip
            inbv = 1_ip
            call dintrv(t,npk,aa,ilo,il1,mflag)
            call dintrv(t,npk,bb,ilo,il2,mflag)
            if ( il2>=np1 ) il2 = n
            do left = il1 , il2
                ta = t(left)
                tb = t(left+1_ip)
                if ( ta/=tb ) then
                    a = max(aa,ta)
                    b = min(bb,tb)
                    bma = 0.5_wp*(b-a)
                    bpa = 0.5_wp*(b+a)
                    do m = 1_ip , mf
                        c1 = bma*gpts(jf+m)
                        gx = -c1 + bpa
                        call dbvalu(t,bcoef,n,k,0_ip,gx,inbv,work,iflag,y2)
                        if (iflag/=0_ip) return
                        gx = c1 + bpa
                        call dbvalu(t,bcoef,n,k,0_ip,gx,inbv,work,iflag,y1)
                        if (iflag/=0_ip) return
                        s(m) = s(m) + (y1+y2)*bma
                    end do
                end if
            end do
            q = 0.0_wp
            do m = 1_ip , mf
                q = q + gwts(jf+m)*s(m)
            end do
            if ( x1>x2 ) q = -q
                bquad = q
                return
            end if
        end if

        iflag = 903_ip ! error return

    end if

    end subroutine dbsqad