/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:41 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dfmin.h" #include #include #include /* PARAMETER translations */ #define FIVE 5.0e0 #define HALF 0.5e0 #define PT95 0.95e0 #define THREE 3.0e0 #define TWO 2.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dfmin( double *x, double *xorf, long *mode, double tol) { static long int _l0, next; static double a, asave, b, bsave, c, c3, d, e, fu, fv, fw, fx, p, q, r, tol1, tol2, toli, u, v, w, xi, xm; static double eps = ZERO; static long np = 0; static long ic = 0; /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2004-11-11 DFMIN Krogh Removed if (..) C3 = statement doing nil. *>> 2000-12-01 DFMIN Krogh Removed unused parameter ONE. *>> 1996-03-30 DFMIN Krogh Added external statement. *>> 1994-10-20 DFMIN Krogh Changes to use M77CON *>> 1994-04-20 DFMIN CLL Edited to make DP & SP files similar. *>> 1987-12-09 DFMIN Lawson Initial code. * * Finds a local minimum of a function f(x) in the closed interval * between A and B. * The function f(x) is defined by code in the calling program. * The user must initially provide A, B, and a tolerance, TOL. * The function will not be evaluated at A or B unless the * minimization search leads to one of these endpoints. * * This subroutine uses reverse communication, i.e., it returns to * the calling program each time it needs to have f() evaluated at a * new value of x. * ------------------------------------------------------------------ * Subroutine Arguments * * X, XORF, MODE [all are inout] * On the initial call to this subroutine to solve a * new problem, the user must set MODE = 0, and must set * X = A * XORF = B * TOL = desired tolerance * A and B denote endpoints defining a closed * interval in which a local minimum is to be found. * Permit A < B or A > B or A = B. If A = B the solution, x = A, * will be found immediately. * On any call the user can set MODE negative to start or stop * detail printing. On such an entry we set NP = -1 - MODE, and * the normal algorithm will not be executed. On the next NP calls * with MODE = 0 or 1, a detail line will be printed. * * TOL [in] is an absolute tolerance on the uncertainty in the final * estimate of the minimizing abcissa, X1. Recommend TOL .ge. 0. * This subr will set TOLI = TOL if TOL > 0., and otherwise sets * TOLI = 0.. The operational tolerance at any trial abcissa, X, * will be * TOL2 = 2 * abs(X) * sqrt(D1MACH(4)) + (2/3) * C3 * where C3 = TOLI/3.D0 + D1MACH(4)**2 where TOLI = max(TOL,0.D0). * * On each return this subr will set MODE to a value in the range * [1:3] to indicate the action needed from the calling program or * the status on termination. * * = 1 means the calling program must evaluate * f(X), store the value in XORF, and then call this * subr again. * * = 2 means normal termination. X contains the point of the * minimum and XORF contains f(X). * * = 3 as for 2 but the requested accuracy could not be obtained. * * = 4 means termination on an error condition: * MODE on entry not in [0:1]. * ------------------------------------------------------------------ * This algorithm is due to Richard. P. Brent. It is presented as * an ALGOL procedure, LOCALMIN, in his 1973 book, * "Algorithms for minimization without derivatives", Prentice-Hall. * Published as subroutine FMIN in Forsythe, Malcolm, and Moler, * "Computer Methods for Mathematical Computations", Prentice-Hall, * 1977. * The current subroutine adapted from F.,M.& M. by C. L. Lawson, and * F. T. Krogh, JPL, Oct 1987, for use in Fortran 77 in the JPL * MATH77 library. The changes improve performance when a minimum is * at an endpoint, and change the user interface to reverse * communication. Unlike the original, the changed algorithm may * evaluate the function at one of the end points. * ------------------------------------------------------------------ * Subprograms referenced: D1MACH, IERM1 * ------------------------------------------------------------------ * 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 AX OR * BX), THEN CONVERGENCE IS SUPERLINEAR, AND USUALLY OF THE ORDER OF * ABOUT 1.324.... * THE FUNCTION F IS NEVER EVALUATED AT TWO POINTS CLOSER TOGETHER * THAN EPS*ABS(FMIN) + (TOLI/3), WHERE EPS IS APPROXIMATELY THE SQUARE * ROOT OF THE RELATIVE MACHINE PRECISION. IF F IS A UNIMODAL * FUNCTION AND THE COMPUTED VALUES OF F ARE ALWAYS UNIMODAL WHEN * SEPARATED BY AT LEAST EPS*ABS(X) + (TOLI/3), THEN FMIN APPROXIMATES * THE ABCISSA OF THE GLOBAL MINIMUM OF F ON THE INTERVAL [AX,BX] WITH * AN ERROR LESS THAN 3*EPS*ABS(FMIN) + TOLI. IF F IS NOT UNIMODAL, * THEN FMIN MAY APPROXIMATE A LOCAL, BUT PERHAPS NON-GLOBAL, MINIMUM TO * THE SAME ACCURACY. (Comments from F., M. & M.) * ------------------------------------------------------------------ *--D replaces "?": ?FMIN * Both versions use IERM1 * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ if (*mode < 0) { np = -1 - *mode; return; } if (np > 0) { np -= 1; printf(" In DFMIN: X= %g , XORF= %g , IC= %ld \n", *x, *xorf, ic); } if (*mode == 1) { switch (next) { case 1: goto L_301; case 2: goto L_302; } } if (*mode != 0) { /* Error: MODE > 1 on entry. */ ierm1( "DFMIN", *mode, 0, "Input value of MODE exceeds 1." , "MODE", *mode, '.' ); *mode = 4; return; } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * C is the squared inverse of the "Golden Ratio", 1.618... * C = 0.381966 * */ c = HALF*(THREE - sqrt( FIVE )); /* EPS IS THE SQUARE ROOT OF THE RELATIVE MACHINE PRECISION. * */ if (eps == ZERO) eps = sqrt( DBL_EPSILON ); /* INITIALIZATION * */ a = *x; b = *xorf; if (a > b) { xi = a; a = b; b = xi; } /* Now we have A .le. B */ asave = a; bsave = b; toli = tol; if (toli <= ZERO) toli = ZERO; c3 = toli/THREE + powi(eps,4); v = a + c*(b - a); w = v; xi = v; ic = 0; e = ZERO; /* Return to calling prog for computation of FX = f(XI) */ *x = xi; *mode = 1; next = 1; return; /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ L_301: ; next = 2; fx = *xorf; fv = fx; fw = fx; /* MAIN LOOP STARTS HERE */ L_20: xm = HALF*(a + b); tol1 = eps*fabs( xi ) + c3; tol2 = TWO*tol1; /* CHECK STOPPING CRITERION * */ if (fabs( xi - xm ) <= (tol2 - HALF*(b - a))) goto L_90; /* IS GOLDEN-SECTION NECESSARY ? * */ if (fabs( e ) <= tol1) goto L_40; /* FIT PARABOLA * */ r = (xi - w)*(fx - fv); q = (xi - v)*(fx - fw); p = (xi - v)*q - (xi - w)*r; q = TWO*(q - r); if (q > ZERO) p = -p; q = fabs( q ); r = e; e = d; /* IS PARABOLA ACCEPTABLE * */ if (fabs( p ) >= fabs( HALF*q*r )) goto L_40; if (p <= q*(a - xi)) goto L_40; if (p >= q*(b - xi)) goto L_40; /* A PARABOLIC INTERPOLATION STEP * */ d = p/q; u = xi + d; /* F MUST NOT BE EVALUATED TOO CLOSE TO A OR B * */ if ((u - a) < tol2 || (b - u) < tol2) d = sign( tol1, xm - xi ); goto L_50; /* A GOLDEN-SECTION STEP * */ L_40: ic += 1; if (ic > 3) { if (a == asave) { if (ic == 4) { e = a + (eps*fabs( a ) + c3) - xi; } else { if (b != w) goto L_45; e = a - xi; } } else if (b == bsave) { if (ic == 4) { e = b - (eps*fabs( b ) + c3) - xi; } else { if (a != w) goto L_45; e = b - xi; } } else { ic = -99; goto L_45; } d = e; goto L_50; } L_45: if (xi >= xm) { e = a - xi; } else { e = b - xi; } d = c*e; /* F MUST NOT BE EVALUATED TOO CLOSE TO XI * */ L_50: if (fabs( d ) >= tol1) { u = xi + d; } else { if (b - a > THREE*tol1) tol1 = tol2; u = fmin( b, fmax( a, xi + sign( PT95*tol1, d ) ) ); } /* Return to calling prog for computation of FU = f(U) * Returning with MODE = 1 and NEXT = 2 */ *x = u; return; /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ L_302: ; fu = *xorf; /* UPDATE A, B, V, W, AND XI * */ if (fu <= fx) { if (fu == fx) { if (b - fmin( xi, u ) > fmax( xi, u ) - a) { b = fmax( xi, u ); } else { a = fmin( xi, u ); } } else { if (u >= xi) { a = xi; } else { b = xi; } } v = w; fv = fw; w = xi; fw = fx; xi = u; fx = fu; goto L_20; } if (u < xi) { a = u; } else { b = u; } if (fu <= fw || w == xi) { v = w; fv = fw; w = u; fw = fu; } else if ((fu <= fv || v == xi) || v == w) { v = u; fv = fu; } goto L_20; /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ L_90: ; *x = xi; *xorf = fx; *mode = 2; if ((b - a > THREE*toli) && (toli != ZERO)) *mode = 3; return; } /* end of function */