Homec4science
Diffusion SPClibs (master) Browse

spclibs/pppack/interv.f90
03540012fedemaster

Last Change
Raw File

interv.f90
View Options

!>
!> @file interv.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> spclibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine interv ( xt, lxt, x, left, mflag )
!*******************************************************************************
!
!! INTERV brackets a real value in an ascending vector of values.
!
! Discussion:
!
! The XT array is a set of increasing values. The goal of the routine
! is to determine the largest index I so that XT(I) <= X.
!
! The routine is designed to be efficient in the common situation
! that it is called repeatedly, with X taken from an increasing
! or decreasing sequence.
!
! This will happen when a piecewise polynomial is to be graphed.
! The first guess for LEFT is therefore taken to be the value
! returned at the previous call and stored in the local variable ILO.
!
! A first check ascertains that ILO < LXT. This is necessary
! since the present call may have nothing to do with the previous
! call. Then, if XT(ILO) < = X < XT(ILO+1), we set LEFT=ILO
! and are done after just three comparisons.
!
! Otherwise, we repeatedly double the difference ISTEP=IHI-ILO
! while also moving ILO and IHI in the direction of X, until
! XT(ILO) < = X < XT(IHI)
! after which we use bisection to get, in addition, ILO+1=IHI.
! LEFT=ILO is then returned.
!
! Modified:
!
! 05 February 2004
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) XT(LXT), a nondecreasing sequence of values.
!
! Input, integer LXT, the dimension of XT.
!
! Input, real ( kind = 8 ) X, the point whose location with
! respect to the sequence XT is to be determined.
!
! Output, integer LEFT, the index of the bracketing value:
! 1 if X < XT(1)
! I if XT(I) <= X < XT(I+1)
! LXT if XT(LXT) <= X
!
! Output, integer MFLAG, indicates whether X lies within the
! range of the data.
! -1: X < XT(1)
! 0: XT(I) <= X < XT(I+1)
! +1: XT(LXT) <= X
!
implicit none
integer lxt
integer left
integer mflag
integer ihi
integer, save :: ilo = 1
integer istep
integer middle
real ( kind = 8 ) x
real ( kind = 8 ) xt(lxt)
!$omp threadprivate(ilo)
ihi = ilo + 1
if ( lxt <= ihi ) then
if ( xt(lxt) <= x ) then
go to 110
end if
if ( lxt <= 1 ) then
mflag = -1
left = 1
return
end if
ilo = lxt - 1
ihi = lxt
end if
if ( xt(ihi) <= x ) then
go to 40
end if
if ( xt(ilo) <= x ) then
mflag = 0
left = ilo
return
end if
!
! Now X < XT(ILO). Decrease ILO to capture X.
!
istep = 1
31 continue
ihi = ilo
ilo = ihi - istep
if ( 1 < ilo ) then
if ( xt(ilo) <= x ) then
go to 50
end if
istep = istep * 2
go to 31
end if
ilo = 1
if ( x < xt(1) ) then
mflag = -1
left = 1
return
end if
go to 50
!
! Now XT(IHI) <= X. Increase IHI to capture X.
!
40 continue
istep = 1
41 continue
ilo = ihi
ihi = ilo + istep
if ( ihi < lxt ) then
if ( x < xt(ihi) ) then
go to 50
end if
istep = istep * 2
go to 41
end if
if ( xt(lxt) <= x ) then
go to 110
end if
!
! Now XT(ILO) < = X < XT(IHI). Narrow the interval.
!
ihi = lxt
50 continue
do
middle = ( ilo + ihi ) / 2
if ( middle == ilo ) then
mflag = 0
left = ilo
return
end if
!
! It is assumed that MIDDLE = ILO in case IHI = ILO+1.
!
if ( xt(middle) <= x ) then
ilo = middle
else
ihi = middle
end if
end do
!
! Set output and return.
!
110 continue
mflag = 1
if ( x == xt(lxt) ) then
mflag = 0
end if
do left = lxt, 1, -1
if ( xt(left) < xt(lxt) ) then
return
end if
end do
return
end
c4science · Help