dbfqad Subroutine

private subroutine dbfqad(f, t, bcoef, n, k, id, x1, x2, tol, quad, iflag, work)

dbfqad computes the integral on (x1,x2) of a product of a function f and the id-th derivative of a k-th order b-spline, using the b-representation (t,bcoef,n,k). (x1,x2) must be a subinterval of t(k) <= x <= t(n+1). an integration routine, dbsgq8 (a modification of gaus8), integrates the product on subintervals of (x1,x2) formed by included (distinct) knots

Reference

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

History

  • 800901 Amos, D. E., (SNLA)
  • 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. Some changes.

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.

Note

Extrapolation is not enabled for this routine.

Arguments

Type IntentOptional Attributes Name
procedure(b1fqad_func) :: f

external function of one argument for the integrand bf(x)=f(x)*dbvalu(t,bcoef,n,k,id,x,inbv,work)
real(kind=wp), intent(in), dimension(n+k) :: t

knot array
real(kind=wp), intent(in), dimension(n) :: bcoef

coefficient array
integer(kind=ip), intent(in) :: n

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

order of b-spline, k >= 1
integer(kind=ip), intent(in) :: id

order of the spline derivative, 0 <= id <= k-1 id=0 gives the spline function
real(kind=wp), intent(in) :: x1

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

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

desired accuracy for the quadrature, suggest 10*dtol < tol <= 0.1 where dtol is the maximum of 1.0e-300 and real(wp) unit roundoff for the machine
real(kind=wp), intent(out) :: quad

integral of bf(x) on (x1,x2)
integer(kind=ip), intent(out) :: iflag

status flag:

  • 0: no errors
  • 1001: k does not satisfy k>=1
  • 1002: n does not satisfy n>=k
  • 1003: d does not satisfy 0<=id<k
  • 1004: x1 or x2 or both do not satisfy t(k)<=x<=t(n+1)
  • 1005: tol is less than dtol or greater than 0.1
real(kind=wp), intent(inout), dimension(:) :: work

work vector of length 3*k

Calls

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

Called by

proc~~dbfqad~~CalledByGraph proc~dbfqad bspline_sub_module::dbfqad 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 dbfqad(f,t,bcoef,n,k,id,x1,x2,tol,quad,iflag,work)

    implicit none

    procedure(b1fqad_func)              :: f      !! external function of one argument for the
                                                  !! integrand `bf(x)=f(x)*dbvalu(t,bcoef,n,k,id,x,inbv,work)`
    integer(ip),intent(in)              :: n      !! length of coefficient array
    integer(ip),intent(in)              :: k      !! order of b-spline, `k >= 1`
    real(wp),dimension(n+k),intent(in)  :: t      !! knot array
    real(wp),dimension(n),intent(in)    :: bcoef  !! coefficient array
    integer(ip),intent(in)              :: id     !! order of the spline derivative, `0 <= id <= k-1`
                                                  !! `id=0` gives the spline function
    real(wp),intent(in)                 :: x1     !! left point of quadrature interval in `t(k) <= x <= t(n+1)`
    real(wp),intent(in)                 :: x2     !! right point of quadrature interval in `t(k) <= x <= t(n+1)`
    real(wp),intent(in)                 :: tol    !! desired accuracy for the quadrature, suggest
                                                  !! `10*dtol < tol <= 0.1` where `dtol` is the maximum
                                                  !! of `1.0e-300` and real(wp) unit roundoff for
                                                  !! the machine
    real(wp),intent(out)                :: quad   !! integral of `bf(x)` on `(x1,x2)`
    real(wp),dimension(:),intent(inout) :: work   !! work vector of length `3*k`
    integer(ip),intent(out)             :: iflag  !! status flag:
                                                  !!
                                                  !! * 0: no errors
                                                  !! * 1001: `k` does not satisfy `k>=1`
                                                  !! * 1002: `n` does not satisfy `n>=k`
                                                  !! * 1003: `d` does not satisfy `0<=id<k`
                                                  !! * 1004: `x1` or `x2` or both do not
                                                  !!   satisfy `t(k)<=x<=t(n+1)`
                                                  !! * 1005: `tol` is less than `dtol`
                                                  !!   or greater than 0.1

    integer(ip) :: inbv,ilo,il1,il2,left,mflag,npk,np1
    real(wp) :: a,aa,ans,b,bb,q,ta,tb,err

    real(wp),parameter :: min_tol = max(epsilon(1.0_wp),1.0e-300_wp) !! minimum allowed `tol`

    iflag = 0_ip
    quad = 0.0_wp
    err = tol
    if ( k<1_ip ) then
        iflag = 1001_ip     ! error
    elseif ( n<k ) then
        iflag = 1002_ip     ! error
    elseif ( id<0_ip .or. id>=k ) then
        iflag = 1003_ip     ! error
    else
        if ( tol>=min_tol .and. tol<=0.1_wp ) then
            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
                    ilo = 1_ip
                    call dintrv(t,npk,aa,ilo,il1,mflag)
                    call dintrv(t,npk,bb,ilo,il2,mflag)
                    if ( il2>=np1 ) il2 = n
                    inbv = 1_ip
                    q = 0.0_wp
                    do left = il1 , il2
                        ta = t(left)
                        tb = t(left+1_ip)
                        if ( ta/=tb ) then
                            a = max(aa,ta)
                            b = min(bb,tb)
                            call dbsgq8(f,t,bcoef,n,k,id,a,b,inbv,err,ans,iflag,work)
                            if ( iflag/=0_ip .and. iflag/=1101_ip ) return
                            q = q + ans
                        end if
                    end do
                    if ( x1>x2 ) q = -q
                    quad = q
                end if
            else
                iflag = 1004_ip  ! error
            end if
        else
            iflag = 1005_ip  ! error
        end if
    end if

    end subroutine dbfqad