dintrv Subroutine

private pure subroutine dintrv(xt, x, ilo, ileft, iright, mflag, inearest)

Returns the indices in xt that bound x, to use for interpolation. If outside the range, then the indices are returned that can be used for extrapolation. Precisely,

         if            x < xt(1)   then ileft=1,   iright=2,    mflag=-1
         if   xt(i) <= x < xt(i+1) then ileft=i,   iright=i+1,  mflag=0
         if   xt(n) <= x           then ileft=n-1, iright=n,    mflag=1

History

  • interv written by carl de boor [5]
  • dintrv author: amos, d. e., (snla) : date written 800901
  • revision date 820801
  • Jacob Williams, 2/24/2015 : updated to free-form Fortran.
  • Jacob Williams, 2/17/2016 : additional refactoring (eliminated GOTOs).
  • Jacob Williams, 2/22/2016 : modified bspline-fortran dintrv routine for linear interpolation/extrapolation use.
  • Jacob Williams, 10/9/2019 : added optional inearest output.

Arguments

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

a knot or break point vector
real(kind=wp), intent(in) :: x

argument
integer, intent(inout) :: ilo

an initialization parameter which must be set to 1 the first time the array xt is processed by dintrv. ilo contains information for efficient processing after the initial call and ilo must not be changed by the user. each dimension requires a distinct ilo parameter.
integer, intent(out) :: ileft

left index
integer, intent(out) :: iright

right index
integer, intent(out) :: mflag

signals when x lies out of bounds
integer, intent(out), optional :: inearest

nearest index

Called by

proc~~dintrv~~CalledByGraph proc~dintrv linear_interpolation_module::dintrv proc~interp_1d linear_interpolation_module::linear_interp_1d%interp_1d proc~interp_1d->proc~dintrv proc~interp_2d linear_interpolation_module::linear_interp_2d%interp_2d proc~interp_2d->proc~dintrv proc~interp_3d linear_interpolation_module::linear_interp_3d%interp_3d proc~interp_3d->proc~dintrv proc~interp_4d linear_interpolation_module::linear_interp_4d%interp_4d proc~interp_4d->proc~dintrv proc~interp_5d linear_interpolation_module::linear_interp_5d%interp_5d proc~interp_5d->proc~dintrv proc~interp_6d linear_interpolation_module::linear_interp_6d%interp_6d proc~interp_6d->proc~dintrv proc~nearest_1d linear_interpolation_module::nearest_interp_1d%nearest_1d proc~nearest_1d->proc~dintrv proc~nearest_2d linear_interpolation_module::nearest_interp_2d%nearest_2d proc~nearest_2d->proc~dintrv proc~nearest_3d linear_interpolation_module::nearest_interp_3d%nearest_3d proc~nearest_3d->proc~dintrv proc~nearest_4d linear_interpolation_module::nearest_interp_4d%nearest_4d proc~nearest_4d->proc~dintrv proc~nearest_5d linear_interpolation_module::nearest_interp_5d%nearest_5d proc~nearest_5d->proc~dintrv proc~nearest_6d linear_interpolation_module::nearest_interp_6d%nearest_6d proc~nearest_6d->proc~dintrv

Source Code

    pure subroutine dintrv(xt,x,ilo,ileft,iright,mflag,inearest)

    implicit none

    real(wp),dimension(:),intent(in) :: xt       !! a knot or break point vector
    real(wp),intent(in)              :: x        !! argument
    integer,intent(inout)            :: ilo      !! an initialization parameter which must be set
                                                 !! to 1 the first time the array `xt` is
                                                 !! processed by dintrv. `ilo` contains information for
                                                 !! efficient processing after the initial call and `ilo`
                                                 !! must not be changed by the user.  each dimension
                                                 !! requires a distinct `ilo` parameter.
    integer,intent(out)              :: ileft    !! left index
    integer,intent(out)              :: iright   !! right index
    integer,intent(out)              :: mflag    !! signals when `x` lies out of bounds
    integer,intent(out),optional     :: inearest !! nearest index

    integer :: ihi, istep, imid, n

    n = size(xt)

    if (n==1) then
        ! this is only allowed for nearest interpolation
        if (present(inearest)) then
            inearest = 1
            return
        end if
    end if

    ihi = ilo + 1
    if ( ihi>=n ) then
        if ( x>=xt(n) ) then
            mflag = 1
            ileft = n-1
            iright= n
            if (present(inearest)) inearest = n
            return
        end if
        if ( n<=1 ) then
            mflag = -1
            ileft = 1
            iright= 2
            if (present(inearest)) inearest = 1
            return
        end if
        ilo = n - 1
        ihi = n
    endif

    if ( x>=xt(ihi) ) then

        ! now x >= xt(ilo). find upper bound
        istep = 1
        do
            ilo = ihi
            ihi = ilo + istep
            if ( ihi>=n ) then
                if ( x>=xt(n) ) then
                    mflag = 1
                    ileft = n-1
                    iright= n
                    if (present(inearest)) inearest = n
                    return
                end if
                ihi = n
            elseif ( x>=xt(ihi) ) then
                istep = istep*2
                cycle
            endif
            exit
        end do

    else

        if ( x>=xt(ilo) ) then
            mflag = 0
            ileft = ilo
            iright= ilo+1
            if (present(inearest)) then
                if ( abs(x-xt(ileft)) <= abs(x-xt(iright)) ) then
                    inearest = ileft
                else
                    inearest = iright
                end if
            end if
            return
        end if
        ! now x <= xt(ihi). find lower bound
        istep = 1
        do
            ihi = ilo
            ilo = ihi - istep
            if ( ilo<=1 ) then
                ilo = 1
                if ( x<xt(1) ) then
                    mflag = -1
                    ileft = 1
                    iright= 2
                    if (present(inearest)) inearest = 1
                    return
                end if
            elseif ( x<xt(ilo) ) then
                istep = istep*2
                cycle
            endif
            exit
        end do

    endif

    ! now xt(ilo) <= x < xt(ihi). narrow the interval
    do
        imid = (ilo+ihi)/2
        if ( imid==ilo ) then
            mflag = 0
            ileft = ilo
            iright= ilo+1
            if (present(inearest)) then
                if ( abs(x-xt(ileft)) <= abs(x-xt(iright)) ) then
                    inearest = ileft
                else
                    inearest = iright
                end if
            end if
            return
        end if
        ! note. it is assumed that imid = ilo in case ihi = ilo+1
        if ( x<xt(imid) ) then
            ihi = imid
        else
            ilo = imid
        endif
    end do

    end subroutine dintrv