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