db1extremum Subroutine

    pure subroutine db1extremum(tx,nx,kx,bcoef,xa,xb,find_max,xopt,fopt,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)                  :: xa       !! one end of the interval
    real(wp),intent(in)                  :: xb       !! other end of the interval
    logical,intent(in)                   :: find_max !! find the maximum (otherwise the minimum)
    real(wp),intent(out)                 :: xopt     !! location of the extremum
    real(wp),intent(out)                 :: fopt     !! value of the spline at `xopt`
    integer(ip),intent(out)              :: iflag    !! status flag (see [[db1min]])

    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(:),allocatable :: xc    !! candidate points
    real(wp),dimension(3*kx) :: work  !! work array for [[dbvalu]]
    real(wp),dimension(kx,1) :: vnikx  !! basis function values from [[dbspvd]]
    real(wp),dimension((kx+1)*(kx+2)/2) :: work2  !! work array for [[dbspvd]]
    integer(ip) :: i       !! counter
    integer(ip) :: m       !! multiplicity of an interior knot
    integer(ip) :: inbv    !! initialization parameter for [[dbvalu]]
    integer(ip) :: iflag2  !! status flag from [[spline_roots]] for \( s'(x) \) (not returned)
    real(wp) :: x1         !! left end of the interval: `min(xa,xb)`
    real(wp) :: x2         !! right end of the interval: `max(xa,xb)`
    real(wp) :: f          !! spline value at a candidate point

    xopt = 0.0_wp
    fopt = 0.0_wp

    call check_root_inputs(nx,kx,iflag); if (iflag/=0_ip) return

    x1 = min(xa,xb)
    x2 = max(xa,xb)
    if (x1<tx(kx) .or. x2>tx(nx+1_ip)) then
        iflag = 4003_ip
        return
    end if

    ! candidates: the end points, the interior knots (the derivative
    ! may be discontinuous there), and the critical points
    xc = [x1, x2, pack(tx(kx+1_ip:nx), tx(kx+1_ip:nx)>x1 .and. tx(kx+1_ip:nx)<x2)]
    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)  ! flat intervals are covered by the knots
        xc = [xc, pack(crit, crit>x1 .and. crit<x2)]
    end if

    inbv = 1_ip
    do i = 1_ip, size(xc,kind=ip)
        call dbvalu(tx,bcoef,nx,kx,0_ip,xc(i),inbv,work,iflag,f)
        if (iflag/=0_ip) return
        if (i==1_ip .or. (find_max .and. f>fopt) .or. (.not. find_max .and. f<fopt)) then
            xopt = xc(i)
            fopt = f
        end if
    end do

    ! the limits from the left at interior knots where the spline is
    ! discontinuous (multiplicity >= kx), in (x1,x2]. Each such knot is
    ! processed at the first of its copies, i, so interval i-1 is the
    ! nonempty interval to its left.
    do i = kx+1_ip, nx
        if (tx(i)<=x1 .or. tx(i)>x2 .or. tx(i)==tx(i-1_ip)) cycle
        m = 1_ip
        do while (i+m<=nx+kx)
            if (tx(i+m)/=tx(i)) exit
            m = m + 1_ip
        end do
        if (m<kx .or. tx(i)>=tx(nx+1_ip)) cycle  ! continuous, or the right end (where
                                                 ! s(x) is already the limit from the left)
        call dbspvd(tx,kx,1_ip,tx(i),i-1_ip,kx,vnikx,work2,iflag2)
        if (iflag2/=0_ip) cycle
        f = dot_product(vnikx(:,1), bcoef(i-kx:i-1_ip))
        if ((find_max .and. f>fopt) .or. (.not. find_max .and. f<fopt)) then
            xopt = tx(i)
            fopt = f
            iflag = 4007_ip
        end if
    end do

    end subroutine db1extremum