subroutine local_min_rc ( a, b, arg, status, value )
!*****************************************************************************80
!
!! local_min_rc() seeks a minimizer of a scalar function of a scalar variable.
!
! Discussion:
!
! This routine seeks an approximation to the point where a function
! F attains a minimum on the interval (A,B).
!
! The method used is a combination of golden section search and
! successive parabolic interpolation. Convergence is never much
! slower than that for a Fibonacci search. If F has a continuous
! second derivative which is positive at the minimum (which is not
! at A or B), then convergence is superlinear, and usually of the
! order of about 1.324...
!
! The routine is a revised version of the Brent local minimization
! algorithm, using reverse communication.
!
! It is worth stating explicitly that this routine will NOT be
! able to detect a minimizer that occurs at either initial endpoint
! A or B. If this is a concern to the user, then the user must
! either ensure that the initial interval is larger, or to check
! the function value at the returned minimizer against the values
! at either endpoint.
!
! Licensing:
!
! This code is distributed under the MIT license.
!
! Modified:
!
! 22 October 2004
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Richard Brent,
! Algorithms for Minimization Without Derivatives,
! Dover, 2002,
! ISBN: 0-486-41998-3,
! LC: QA402.5.B74.
!
! David Kahaner, Cleve Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1989,
! ISBN: 0-13-627258-4,
! LC: TA345.K34.
!
! Parameters
!
! Input/output, real ( kind = rk ) A, B. On input, the left and right
! endpoints of the initial interval. On output, the lower and upper
! bounds for an interval containing the minimizer. It is required
! that A < B.
!
! Output, real ( kind = rk ) ARG, the currently considered point. The user
! does not need to initialize this value. On return with STATUS positive,
! the user is requested to evaluate the function at ARG, and return
! the value in VALUE. On return with STATUS zero, ARG is the routine's
! estimate for the function minimizer.
!
! Input/output, integer ( kind = 4 ) STATUS, used to communicate between
! the user and the routine. The user only sets STATUS to zero on the first
! call, to indicate that this is a startup call. The routine returns STATUS
! positive to request that the function be evaluated at ARG, or returns
! STATUS as 0, to indicate that the iteration is complete and that
! ARG is the estimated minimizer.
!
! Input, real ( kind = rk ) VALUE, the function value at ARG, as requested
! by the routine on the previous call.
!
! Local parameters:
!
! C is the squared inverse of the golden ratio.
!
! EPS is the square root of the relative machine precision.
!
implicit none
integer, parameter :: rk = kind ( 1.0D+00 )
real ( kind = rk ) a
real ( kind = rk ) arg
real ( kind = rk ) b
real ( kind = rk ), save :: c
real ( kind = rk ), save :: d
real ( kind = rk ), save :: e
real ( kind = rk ), save :: eps
real ( kind = rk ), save :: fu
real ( kind = rk ), save :: fv
real ( kind = rk ), save :: fw
real ( kind = rk ), save :: fx
real ( kind = rk ), save :: midpoint
real ( kind = rk ), save :: p
real ( kind = rk ), save :: q
real ( kind = rk ), save :: r
integer ( kind = 4 ) status
real ( kind = rk ), save :: tol
real ( kind = rk ), save :: tol1
real ( kind = rk ), save :: tol2
real ( kind = rk ), save :: u
real ( kind = rk ), save :: v
real ( kind = rk ) value
real ( kind = rk ), save :: w
real ( kind = rk ), save :: x
!
! STATUS (INPUT) = 0, startup.
!
if ( status == 0 ) then
if ( b <= a ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'local_min_rc(): Fatal error!'
write ( *, '(a)' ) ' A < B is required, but'
write ( *, '(a,g14.6)' ) ' A = ', a
write ( *, '(a,g14.6)' ) ' B = ', b
status = -1
stop 1
end if
c = 0.5D+00 * ( 3.0D+00 - sqrt ( 5.0D+00 ) )
eps = sqrt ( epsilon ( eps ) )
tol = epsilon ( tol )
v = a + c * ( b - a )
w = v
x = v
e = 0.0D+00
status = 1
arg = x
return
!
! STATUS (INPUT) = 1, return with initial function value of FX.
!
else if ( status == 1 ) then
fx = value
fv = fx
fw = fx
!
! STATUS (INPUT) = 2 or more, update the data.
!
else if ( 2 <= status ) then
fu = value
if ( fu <= fx ) then
if ( x <= u ) then
a = x
else
b = x
end if
v = w
fv = fw
w = x
fw = fx
x = u
fx = fu
else
if ( u < x ) then
a = u
else
b = u
end if
if ( fu <= fw .or. w == x ) then
v = w
fv = fw
w = u
fw = fu
else if ( fu <= fv .or. v == x .or. v == w ) then
v = u
fv = fu
end if
end if
end if
!
! Take the next step.
!
midpoint = 0.5D+00 * ( a + b )
tol1 = eps * abs ( x ) + tol / 3.0D+00
tol2 = 2.0D+00 * tol1
!
! If the stopping criterion is satisfied, we can exit.
!
if ( abs ( x - midpoint ) <= ( tol2 - 0.5D+00 * ( b - a ) ) ) then
status = 0
return
end if
!
! Is golden-section necessary?
!
if ( abs ( e ) <= tol1 ) then
if ( midpoint <= x ) then
e = a - x
else
e = b - x
end if
d = c * e
!
! Consider fitting a parabola.
!
else
r = ( x - w ) * ( fx - fv )
q = ( x - v ) * ( fx - fw )
p = ( x - v ) * q - ( x - w ) * r
q = 2.0D+00 * ( q - r )
if ( 0.0D+00 < q ) then
p = - p
end if
q = abs ( q )
r = e
e = d
!
! Choose a golden-section step if the parabola is not advised.
!
if ( &
( abs ( 0.5D+00 * q * r ) <= abs ( p ) ) .or. &
( p <= q * ( a - x ) ) .or. &
( q * ( b - x ) <= p ) ) then
if ( midpoint <= x ) then
e = a - x
else
e = b - x
end if
d = c * e
!
! Choose a parabolic interpolation step.
!
else
d = p / q
u = x + d
if ( ( u - a ) < tol2 ) then
d = sign ( tol1, midpoint - x )
end if
if ( ( b - u ) < tol2 ) then
d = sign ( tol1, midpoint - x )
end if
end if
end if
!
! F must not be evaluated too close to X.
!
if ( tol1 <= abs ( d ) ) then
u = x + d
end if
if ( abs ( d ) < tol1 ) then
u = x + sign ( tol1, d )
end if
!
! Request value of F(U).
!
arg = u
status = status + 1
return
end