subroutine SFMIN(X,XORF,MODE,TOL) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. c>> 2004-11-11 SFMIN Krogh Removed if (..) C3 = statement doing nil. c>> 2000-12-01 SFMIN Krogh Removed unused parameter ONE. c>> 1996-03-30 SFMIN Krogh Added external statement. C>> 1994-10-20 SFMIN Krogh Changes to use M77CON C>> 1994-04-20 SFMIN CLL Edited to make DP & SP files similar. C>> 1987-12-09 SFMIN Lawson Initial code. c C Finds a local minimum of a function f(x) in the closed interval c between A and B. c The function f(x) is defined by code in the calling program. c The user must initially provide A, B, and a tolerance, TOL. c The function will not be evaluated at A or B unless the c minimization search leads to one of these endpoints. c c This subroutine uses reverse communication, i.e., it returns to c the calling program each time it needs to have f() evaluated at a c new value of x. c ------------------------------------------------------------------ c Subroutine Arguments c c X, XORF, MODE [all are inout] c On the initial call to this subroutine to solve a c new problem, the user must set MODE = 0, and must set c X = A c XORF = B c TOL = desired tolerance c A and B denote endpoints defining a closed c interval in which a local minimum is to be found. c Permit A < B or A > B or A = B. If A = B the solution, x = A, c will be found immediately. c On any call the user can set MODE negative to start or stop c detail printing. On such an entry we set NP = -1 - MODE, and c the normal algorithm will not be executed. On the next NP calls c with MODE = 0 or 1, a detail line will be printed. c c TOL [in] is an absolute tolerance on the uncertainty in the final c estimate of the minimizing abcissa, X1. Recommend TOL .ge. 0. c This subr will set TOLI = TOL if TOL > 0., and otherwise sets c TOLI = 0.. The operational tolerance at any trial abcissa, X, c will be c TOL2 = 2 * abs(X) * sqrt(R1MACH(4)) + (2/3) * C3 c where C3 = TOLI/3.E0 + R1MACH(4)**2 where TOLI = max(TOL,0.E0). c c On each return this subr will set MODE to a value in the range c [1:3] to indicate the action needed from the calling program or c the status on termination. c c = 1 means the calling program must evaluate c f(X), store the value in XORF, and then call this c subr again. c c = 2 means normal termination. X contains the point of the c minimum and XORF contains f(X). c c = 3 as for 2 but the requested accuracy could not be obtained. c c = 4 means termination on an error condition: c MODE on entry not in [0:1]. c ------------------------------------------------------------------ c This algorithm is due to Richard. P. Brent. It is presented as c an ALGOL procedure, LOCALMIN, in his 1973 book, c "Algorithms for minimization without derivatives", Prentice-Hall. c Published as subroutine FMIN in Forsythe, Malcolm, and Moler, c "Computer Methods for Mathematical Computations", Prentice-Hall, c 1977. c The current subroutine adapted from F.,M.& M. by C. L. Lawson, and c F. T. Krogh, JPL, Oct 1987, for use in Fortran 77 in the JPL c MATH77 library. The changes improve performance when a minimum is c at an endpoint, and change the user interface to reverse c communication. Unlike the original, the changed algorithm may c evaluate the function at one of the end points. c ------------------------------------------------------------------ c Subprograms referenced: R1MACH, IERM1 c ------------------------------------------------------------------ c THE METHOD USED IS A COMBINATION OF GOLDEN SECTION SEARCH AND C SUCCESSIVE PARABOLIC INTERPOLATION. CONVERGENCE IS NEVER MUCH SLOWER C THAN THAT FOR A FIBONACCI SEARCH. IF F HAS A CONTINUOUS SECOND C DERIVATIVE WHICH IS POSITIVE AT THE MINIMUM (WHICH IS NOT AT AX OR C BX), THEN CONVERGENCE IS SUPERLINEAR, AND USUALLY OF THE ORDER OF C ABOUT 1.324.... c THE FUNCTION F IS NEVER EVALUATED AT TWO POINTS CLOSER TOGETHER C THAN EPS*ABS(FMIN) + (TOLI/3), WHERE EPS IS APPROXIMATELY THE SQUARE C ROOT OF THE RELATIVE MACHINE PRECISION. IF F IS A UNIMODAL C FUNCTION AND THE COMPUTED VALUES OF F ARE ALWAYS UNIMODAL WHEN C SEPARATED BY AT LEAST EPS*ABS(X) + (TOLI/3), THEN FMIN APPROXIMATES C THE ABCISSA OF THE GLOBAL MINIMUM OF F ON THE INTERVAL [AX,BX] WITH C AN ERROR LESS THAN 3*EPS*ABS(FMIN) + TOLI. IF F IS NOT UNIMODAL, C THEN FMIN MAY APPROXIMATE A LOCAL, BUT PERHAPS NON-GLOBAL, MINIMUM TO C THE SAME ACCURACY. (Comments from F., M. & M.) c ------------------------------------------------------------------ c--S replaces "?": ?FMIN c Both versions use IERM1 c ------------------------------------------------------------------ external R1MACH real X,XORF,TOL, R1MACH real A,B,C,D,E,EPS,XM,P,Q,R,TOL1,TOL2,U,V,W real FU,FV,FW,FX,XI, TOLI, C3 real ASAVE, BSAVE real HALF, TWO, THREE, FIVE, ZERO, PT95 integer MODE, NEXT, IC, NP parameter(HALF = 0.5E0, TWO = 2.0E0, THREE = 3.0E0) parameter(FIVE = 5.0E0, ZERO = 0.0E0, PT95 = 0.95E0) save data EPS / ZERO /, NP / 0 /, IC / 0 / c ------------------------------------------------------------------ IF (MODE .LT. 0) THEN NP = -1 - MODE RETURN END IF IF (NP .GT. 0) THEN NP = NP - 1 print*, ' In SFMIN: X= ',X,', XORF= ',XORF,', IC= ',IC END IF if(MODE .eq. 1) then go to (301,302), NEXT endif if( MODE .ne. 0) then c Error: MODE > 1 on entry. call IERM1('SFMIN',MODE,0, * 'Input value of MODE exceeds 1.','MODE',MODE,'.') MODE = 4 return endif c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C is the squared inverse of the "Golden Ratio", 1.618... C C = 0.381966 c C = HALF*(THREE - sqrt(FIVE)) C C EPS IS THE SQUARE ROOT OF THE RELATIVE MACHINE PRECISION. C if (EPS .eq. ZERO) EPS = sqrt(R1MACH(4)) C C INITIALIZATION C A = X B = XORF if( A .gt. B) then XI = A A = B B = XI endif c Now we have A .le. B ASAVE = A BSAVE = B TOLI = TOL if(TOLI .le. ZERO) TOLI = ZERO C3 = TOLI/THREE + EPS ** 4 V = A + C*(B - A) W = V XI = V IC = 0 E = ZERO c Return to calling prog for computation of FX = f(XI) X = XI MODE = 1 NEXT = 1 return c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 301 continue NEXT = 2 FX = XORF FV = FX FW = FX C MAIN LOOP STARTS HERE 20 XM = HALF*(A + B) TOL1 = EPS*abs(XI) + C3 TOL2 = TWO*TOL1 C CHECK STOPPING CRITERION C if (abs(XI - XM) .le. (TOL2 - HALF*(B - A))) go to 90 C C IS GOLDEN-SECTION NECESSARY ? C if (abs(E) .le. TOL1) go to 40 C C FIT PARABOLA C R = (XI - W)*(FX - FV) Q = (XI - V)*(FX - FW) P = (XI - V)*Q - (XI - W)*R Q = TWO*(Q - R) if (Q .gt. ZERO) P = -P Q = abs(Q) R = E E = D C C IS PARABOLA ACCEPTABLE C if (abs(P) .ge. abs(HALF*Q*R)) go to 40 if (P .le. Q*(A - XI)) go to 40 if (P .ge. Q*(B - XI)) go to 40 C C A PARABOLIC INTERPOLATION STEP C D = P/Q U = XI + D C C F MUST NOT BE EVALUATED TOO CLOSE TO A OR B C if ((U - A) .lt. TOL2 .or. (B - U) .lt. TOL2) * D = sign(TOL1, XM - XI) go to 50 C C A GOLDEN-SECTION STEP C 40 IC = IC + 1 if (IC .gt. 3) then if (A .eq. ASAVE) then if (IC .eq. 4) then E = A + (EPS * abs(A) + C3) - XI else if (B .ne. W) go to 45 E = A - XI end if else if (B .eq. BSAVE) then if (IC .eq. 4) then E = B - (EPS * abs(B) + C3) - XI else if (A .ne. W) go to 45 E = B - XI end if else IC = -99 go to 45 end if D = E go to 50 end if 45 if (XI .ge. XM) then E = A - XI else E = B - XI endif D = C*E C C F MUST NOT BE EVALUATED TOO CLOSE TO XI C 50 if (abs(D) .ge. TOL1) then U = XI + D else if (B - A .gt. THREE * TOL1) TOL1 = TOL2 U = min(B, max(A, XI + sign(PT95*TOL1, D))) endif c c Return to calling prog for computation of FU = f(U) c Returning with MODE = 1 and NEXT = 2 X = U return c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 302 continue FU = XORF C C UPDATE A, B, V, W, AND XI C if (FU .le. FX) then if (FU .eq. FX) then if (B - min(XI, U) .gt. max(XI, U) - A) then B = max(XI, U) else A = min(XI, U) end if else if (U .ge. XI) then A = XI else B = XI endif end if V = W FV = FW W = XI FW = FX XI = U FX = FU go to 20 endif c if (U .lt. XI) then A = U else B = U endif if (FU .le. FW .or. W .eq. XI) then V = W FV = FW W = U FW = FU elseif (FU .le. FV .or. V .eq. XI .or. V .eq. W) then V = U FV = FU endif go to 20 c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 90 continue X = XI XORF = FX MODE = 2 if ((B - A .gt. THREE*TOLI) .and. (TOLI .ne. ZERO)) MODE = 3 return end