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
real(kind=wp), intent(inout), dimension(:) :: work

work vector of length 3*k

Calls

proc~~dbfqad~~CallsGraph proc~dbfqad dbfqad proc~dbsgq8 dbsgq8 proc~dbfqad->proc~dbsgq8 proc~dbspvd dbspvd proc~dbfqad->proc~dbspvd proc~dintrv dintrv proc~dbfqad->proc~dintrv proc~dbvalu dbvalu proc~dbsgq8->proc~dbvalu proc~dbspvn dbspvn proc~dbspvd->proc~dbspvn proc~get_temp_x_for_extrap get_temp_x_for_extrap proc~dintrv->proc~get_temp_x_for_extrap proc~dbvalu->proc~dintrv

Called by

proc~~dbfqad~~CalledByGraph proc~dbfqad dbfqad proc~db1fqad db1fqad proc~db1fqad->proc~dbfqad proc~fintegral_1d 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`
    real(wp),parameter :: narrow = 2.0_wp**(-40) !! a subinterval narrower than this relative to `|x|`
                                                 !! is integrated with a fixed Gauss rule

    iflag = 0_ip
    quad = 0.0_wp
    err = tol
    if ( k<1_ip ) then
        iflag = 1001_ip     ! error
    else if ( n<k ) then
        iflag = 1002_ip     ! error
    else if ( 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)
                            if ( abs(b-a) <= narrow*max(abs(a),abs(b)) ) then
                                ! too narrow relative to |x| for dbsgq8 (which returns
                                ! 1101, or 0 if a/b rounds to 1)
                                call gauss8(a,b,left,ans,iflag)
                            else
                                call dbsgq8(f,t,bcoef,n,k,id,a,b,inbv,err,ans,iflag,work)
                                if ( iflag==1101_ip ) call gauss8(a,b,left,ans,iflag)
                            end if
                            if ( iflag/=0_ip ) return
                            q = q + ans
                        end if
                    end do
                    if ( x1>x2 ) q = -q
                    quad = q
                else
                    iflag = 1004_ip  ! error
                end if
            else
                iflag = 1004_ip  ! error
            end if
        else
            iflag = 1005_ip  ! error
        end if
    end if

    contains

        subroutine gauss8(a,b,l,ans,iflag)

        !! 8-point Gauss-Legendre rule on `[a,b]`, for a subinterval that is too
        !! narrow (relative to `|x|`) for [[dbsgq8]]. `[a,b]` is inside knot
        !! interval `l`, where the spline is a polynomial, so the rule is exact
        !! for the spline part for orders up to 16 (and very accurate for a
        !! smooth `f` over such a narrow interval).
        !!
        !! The nodes may not be representable (the interval may be only a few
        !! ulps wide), so the spline is evaluated at the exact offsets from the
        !! midpoint, from its Taylor expansion there. `f` is evaluated at the
        !! nearest representable points.
        !!
        !!@note For orders above 16, the rule is not exact for the spline part
        !!      (but is still accurate to better than about `1e-8`). The
        !!      derivatives in the Taylor expansion can overflow if the knot
        !!      interval is extremely small (e.g., about `1e-200` wide).

        real(wp),intent(in)     :: a      !! left end of the subinterval
        real(wp),intent(in)     :: b      !! right end of the subinterval
        integer(ip),intent(in)  :: l      !! knot interval containing `[a,b]`: `t(l) <= a < b <= t(l+1)`
        real(wp),intent(out)    :: ans    !! the integral
        integer(ip),intent(out) :: iflag  !! status flag from [[dbspvd]]

        real(wp),dimension(4),parameter :: xg = [gl8_x1, gl8_x2, gl8_x3, gl8_x4]  !! nodes (positive half)
        real(wp),dimension(4),parameter :: wg = [gl8_w1, gl8_w2, gl8_w3, gl8_w4]  !! weights
        real(wp),dimension(k,k) :: vnikx  !! basis functions and their derivatives at the midpoint
        real(wp),dimension((k+1)*(k+2)/2) :: wrk  !! work array for [[dbspvd]]
        real(wp),dimension(0:k-1) :: dc   !! Taylor coefficients of the spline (or its derivative)
                                          !! about the midpoint
        real(wp) :: c   !! midpoint of the subinterval
        real(wp) :: h   !! half width of the subinterval
        real(wp) :: r   !! rounding error of `c` (the true midpoint is `c+r`)
        real(wp) :: dx  !! offset of a node from `c`
        real(wp) :: fs  !! spline (or derivative) value at a node
        integer(ip) :: i   !! counter
        integer(ip) :: j   !! side of the midpoint (-1 or 1)
        integer(ip) :: m   !! degree of the Taylor expansion: `k-1-id`
        integer(ip) :: jj  !! counter

        c = a + 0.5_wp*(b-a)
        h = 0.5_wp*(b-a)
        r = (a - c) + h  ! the true midpoint minus the rounded one
        m = k - 1_ip - id
        call dbspvd(t,k,k,c,l,k,vnikx,wrk,iflag)
        if (iflag/=0_ip) return
        do jj = 0_ip, m
            dc(jj) = dot_product(vnikx(:,id+jj+1_ip), bcoef(l-k+1_ip:l))
            if (jj>1_ip) dc(jj) = dc(jj) / product([(real(i,wp), i=2_ip,jj)])
        end do
        ans = 0.0_wp
        do i = 1_ip, 4_ip
            do j = -1_ip, 1_ip, 2_ip
                dx = r + real(j,wp)*h*xg(i)
                fs = dc(m)
                do jj = m-1_ip, 0_ip, -1_ip
                    fs = fs*dx + dc(jj)
                end do
                ans = ans + wg(i) * f(min(b, max(a, c + dx))) * fs
            end do
        end do
        ans = h * ans

        end subroutine gauss8

    end subroutine dbfqad