db1root Subroutine

public pure subroutine db1root(tx, nx, kx, bcoef, a, roots, iflag)

Finds all the real roots of on , where is a 1D b-spline of order kx (kx >= 2).

Each knot interval is examined separately:

  • By the convex hull property, if the kx coefficients that are nonzero on an interval are all above (or all below) a, the interval has no root and is skipped.
  • Otherwise, the polynomial piece is converted to Bernstein form (directly from the b-spline coefficients by blossoming, which is numerically stable for any order). Pieces of degree one or two are solved in closed form. Higher degree roots are isolated by Bernstein subdivision (Descartes' rule of signs) and refined with a safeguarded Newton iteration.
  • Roots of even multiplicity (where the spline touches a without crossing it) are found as the roots of at which is within a small tolerance.

If over a whole knot interval, the roots inside it are not returned (its end points are), and iflag is set to 4004.

At an interior knot where the spline is discontinuous, is the value from the right (as for db1val), so a root of the piece on the left at the knot is not a root of the spline.

Reference

  • K. Morken, M. Reimers, "An unconditionally convergent method for computing zeros of splines and polynomials", Mathematics of Computation, 76 (2007), 845-865.

See also

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in), dimension(nx+kx) :: tx

knot array
integer(kind=ip), intent(in) :: nx

number of b-spline coefficients
integer(kind=ip), intent(in) :: kx

order of the b-spline, kx >= 2
real(kind=wp), intent(in), dimension(nx) :: bcoef

b-spline coefficient array
real(kind=wp), intent(in) :: a

the value to solve for
real(kind=wp), intent(out), dimension(:), allocatable :: roots

the roots, in increasing order (size 0 if there are none)
integer(kind=ip), intent(out) :: iflag

status flag:

  • 0: no errors
  • 4001: kx does not satisfy kx >= 2
  • 4002: nx does not satisfy nx >= kx
  • 4004: over a whole knot interval (other roots are still returned)

Calls

proc~~db1root~~CallsGraph proc~db1root db1root proc~check_root_inputs check_root_inputs proc~db1root->proc~check_root_inputs proc~dbvalu dbvalu proc~db1root->proc~dbvalu proc~in_flat_region in_flat_region proc~db1root->proc~in_flat_region proc~root_dup_tol root_dup_tol proc~db1root->proc~root_dup_tol proc~sort_unique sort_unique proc~db1root->proc~sort_unique proc~spline_derivative spline_derivative proc~db1root->proc~spline_derivative proc~spline_roots spline_roots proc~db1root->proc~spline_roots proc~dintrv dintrv proc~dbvalu->proc~dintrv proc~in_flat_region->proc~dintrv proc~spline_roots->proc~dbvalu proc~spline_roots->proc~in_flat_region proc~spline_roots->proc~root_dup_tol proc~spline_roots->proc~sort_unique proc~get_temp_x_for_extrap get_temp_x_for_extrap proc~dintrv->proc~get_temp_x_for_extrap

Called by

proc~~db1root~~CalledByGraph proc~db1root db1root proc~roots_1d bspline_1d%roots_1d proc~roots_1d->proc~db1root proc~inverse_1d bspline_1d%inverse_1d proc~inverse_1d->proc~roots_1d

Source Code

    pure subroutine db1root(tx,nx,kx,bcoef,a,roots,iflag)

    implicit none

    integer(ip),intent(in)                        :: nx     !! number of b-spline coefficients
    integer(ip),intent(in)                        :: kx     !! order of the b-spline, `kx >= 2`
    real(wp),dimension(nx+kx),intent(in)          :: tx     !! knot array
    real(wp),dimension(nx),intent(in)             :: bcoef  !! b-spline coefficient array
    real(wp),intent(in)                           :: a      !! the value to solve for
    real(wp),dimension(:),allocatable,intent(out) :: roots  !! the roots, in increasing order
                                                            !! (size 0 if there are none)
    integer(ip),intent(out)                       :: iflag  !! status flag:
                                                            !!
                                                            !! * 0: no errors
                                                            !! * 4001: `kx` does not satisfy `kx >= 2`
                                                            !! * 4002: `nx` does not satisfy `nx >= kx`
                                                            !! * 4004: \( s(x) = a \) over a whole knot
                                                            !!   interval (other roots are still returned)

    real(wp),dimension(:),allocatable :: crit  !! critical points (roots of \( s'(x) \))
    real(wp),dimension(:),allocatable :: td    !! knots of \( s'(x) \)
    real(wp),dimension(:),allocatable :: d     !! coefficients of \( s'(x) \)
    real(wp),dimension(3*kx) :: work  !! work array for [[dbvalu]]
    logical,dimension(:),allocatable :: keep  !! roots to keep (those not replaced
                                              !! by a critical point)
    integer(ip) :: i       !! counter
    integer(ip) :: j       !! counter
    integer(ip) :: inbv    !! initialization parameter for [[dbvalu]]
    integer(ip) :: iflag2  !! status flag from [[dbvalu]] and [[spline_roots]]
                           !! for \( s'(x) \) (not returned)
    real(wp) :: f          !! spline value \( s(x) \)
    real(wp) :: f2         !! second derivative \( s''(x) \) at a critical point
    real(wp) :: ftol       !! absolute tolerance on \( |s(x)-a| \)
    real(wp) :: xtol       !! distance below which two roots may be merged
                           !! into a single root of even multiplicity
    real(wp) :: dmax       !! distance from a root of even multiplicity within which
                           !! other roots are replaced by it
    real(wp) :: band_max   !! upper bound on `dmax` (for a root of higher multiplicity,
                           !! where \( s''(x) = 0 \))

    call check_root_inputs(nx,kx,iflag)
    if (iflag/=0_ip) then
        allocate(roots(0))
        return
    end if

    call spline_roots(tx,bcoef,nx,kx,a,roots,iflag)

    ! roots of even multiplicity: critical points where s(x) is equal to a
    if (kx>=3_ip) then
        call spline_derivative(tx,bcoef,nx,kx,td,d)
        call spline_roots(td,d,nx-1_ip,kx-1_ip,0.0_wp,crit,iflag2)
        ftol = root_ftol * max(abs(a),maxval(abs(bcoef)))
        ! (the domain width is scaled before the subtraction, so it does not overflow)
        xtol = (sqrt(epsilon(1.0_wp))*tx(nx+1_ip) - sqrt(epsilon(1.0_wp))*tx(kx)) + &
               root_dup_tol(tx(kx),tx(nx+1_ip))
        band_max = sqrt(sqrt(epsilon(1.0_wp)))*tx(nx+1_ip) - sqrt(sqrt(epsilon(1.0_wp)))*tx(kx)
        inbv = 1_ip
        do i = 1_ip, size(crit,kind=ip)
            call dbvalu(tx,bcoef,nx,kx,0_ip,crit(i),inbv,work,iflag2,f)
            if (iflag2/=0_ip) cycle
            if (abs(f-a)<=ftol) then
                if (iflag==4004_ip) then
                    if (in_flat_region(tx,bcoef,nx,kx,a,ftol,crit(i))) cycle
                end if
                ! near a root of even multiplicity, the roots found above
                ! are only accurate to about sqrt(eps), but the critical
                ! point is accurate to about eps: replace any of them that
                ! are within the band around it where s(x) is within the
                ! noise level of `a`. For a double root, the half width
                ! of the band where |s(x)-a| <= ftol is sqrt(2*ftol/|s''|).
                ! Twice this is used, where |s(x)-a| <= 4*ftol.
                call dbvalu(tx,bcoef,nx,kx,2_ip,crit(i),inbv,work,iflag2,f2)
                dmax = band_max
                if (iflag2==0_ip .and. abs(f2)>0.0_wp) dmax = min(dmax, 2.0_wp*sqrt(2.0_wp*ftol/abs(f2)))
                dmax = max(dmax, xtol)
                allocate(keep(size(roots)))
                do j = 1_ip, size(roots,kind=ip)
                    keep(j) = .true.
                    if (abs(roots(j)-crit(i))<=dmax) then
                        ! check that the spline is within the band between
                        ! them (so they are not distinct roots)
                        call dbvalu(tx,bcoef,nx,kx,0_ip,0.5_wp*(roots(j)+crit(i)),inbv,work,iflag2,f)
                        keep(j) = .not. (iflag2==0_ip .and. abs(f-a)<=4.0_wp*ftol)
                    end if
                end do
                roots = [pack(roots, keep), crit(i)]
                deallocate(keep)
            end if
        end do
        call sort_unique(roots, root_dup_tol(tx(kx),tx(nx+1_ip)))

        ! a root of even multiplicity may be split by rounding into a
        ! close pair of simple roots: merge them if the spline is
        ! within the noise level of `a` between them
        i = 1_ip
        do while (i<size(roots,kind=ip))
            if (roots(i+1_ip)-roots(i)<=xtol) then
                call dbvalu(tx,bcoef,nx,kx,0_ip,0.5_wp*(roots(i)+roots(i+1_ip)),inbv,work,iflag2,f)
                if (iflag2==0_ip .and. abs(f-a)<=ftol) then
                    roots(i) = 0.5_wp*(roots(i)+roots(i+1_ip))
                    roots = [roots(1_ip:i), roots(i+2_ip:)]
                    cycle
                end if
            end if
            i = i + 1_ip
        end do
    end if

    end subroutine db1root