# include # include # include # include using namespace std; # include "local_min_rc.hpp" //****************************************************************************80 double local_min_rc ( double &a, double &b, int &status, double value ) //****************************************************************************80 // // Purpose: // // 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 GNU LGPL license. // // Modified: // // 29 May 2021 // // 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, double &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. // // Input/output, int &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, double VALUE, the function value at ARG, as requested // by the routine on the previous call. // // Output, double LOCAL_MIN_RC, the currently considered point. // On return with STATUS positive, the user is requested to evaluate the // function at this point, and return the value in VALUE. On return with // STATUS zero, this is the routine's estimate for the function minimizer. // // Local: // // double C: the squared inverse of the golden ratio. // // double EPS: the square root of the relative machine precision. // { static double arg; static double c; static double d; static double e; static double eps; static double fu; static double fv; static double fw; static double fx; static double midpoint; static double p; static double q; static double r; static double tol; static double tol1; static double tol2; static double u; static double v; static double w; static double x; // // STATUS (INPUT) = 0, startup. // if ( status == 0 ) { if ( b <= a ) { cerr << "\n"; cerr << "local_min_rc(): Fatal error!\n"; cerr << " A < B is required, but\n"; cerr << " A = " << a << "\n"; cerr << " B = " << b << "\n"; status = -1; exit ( 1 ); } c = 0.5 * ( 3.0 - sqrt ( 5.0 ) ); eps = sqrt ( DBL_EPSILON ); tol = DBL_EPSILON; v = a + c * ( b - a ); w = v; x = v; e = 0.0; status = 1; arg = x; return arg; } // // STATUS (INPUT) = 1, return with initial function value of FX. // else if ( status == 1 ) { fx = value; fv = fx; fw = fx; } // // STATUS (INPUT) = 2 or more, update the data. // else if ( 2 <= status ) { fu = value; if ( fu <= fx ) { if ( x <= u ) { a = x; } else { b = x; } v = w; fv = fw; w = x; fw = fx; x = u; fx = fu; } else { if ( u < x ) { a = u; } else { b = u; } if ( fu <= fw || w == x ) { v = w; fv = fw; w = u; fw = fu; } else if ( fu <= fv || v == x || v == w ) { v = u; fv = fu; } } } // // Take the next step. // midpoint = 0.5 * ( a + b ); tol1 = eps * fabs ( x ) + tol / 3.0; tol2 = 2.0 * tol1; // // If the stopping criterion is satisfied, we can exit. // if ( fabs ( x - midpoint ) <= ( tol2 - 0.5 * ( b - a ) ) ) { status = 0; return arg; } // // Is golden-section necessary? // if ( fabs ( e ) <= tol1 ) { if ( midpoint <= x ) { e = a - x; } else { e = b - x; } 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.0 * ( q - r ); if ( 0.0 < q ) { p = - p; } q = fabs ( q ); r = e; e = d; // // Choose a golden-section step if the parabola is not advised. // if ( ( fabs ( 0.5 * q * r ) <= fabs ( p ) ) || ( p <= q * ( a - x ) ) || ( q * ( b - x ) <= p ) ) { if ( midpoint <= x ) { e = a - x; } else { e = b - x; } d = c * e; } // // Choose a parabolic interpolation step. // else { d = p / q; u = x + d; if ( ( u - a ) < tol2 ) { d = copysign ( tol1, midpoint - x ); } if ( ( b - u ) < tol2 ) { d = copysign ( tol1, midpoint - x ); } } } // // F must not be evaluated too close to X. // if ( tol1 <= fabs ( d ) ) { u = x + d; } if ( fabs ( d ) < tol1 ) { u = x + copysign ( tol1, d ); } // // Request value of F(U). // arg = u; status = status + 1; return arg; }