subroutine DFMIN(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 DFMIN Krogh Removed if (..) C3 = statement doing nil. c>> 2000-12-01 DFMIN Krogh Removed unused parameter ONE. c>> 1996-03-30 DFMIN Krogh Added external statement. C>> 1994-10-20 DFMIN Krogh Changes to use M77CON C>> 1994-04-20 DFMIN CLL Edited to make DP & SP files similar. C>> 1987-12-09 DFMIN 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(D1MACH(4)) + (2/3) * C3 c where C3 = TOLI/3.D0 + D1MACH(4)**2 where TOLI = max(TOL,0.D0). 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: D1MACH, 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--D replaces "?": ?FMIN c Both versions use IERM1 c ------------------------------------------------------------------ external D1MACH double precision X,XORF,TOL, D1MACH double precision A,B,C,D,E,EPS,XM,P,Q,R,TOL1,TOL2,U,V,W double precision FU,FV,FW,FX,XI, TOLI, C3 double precision ASAVE, BSAVE double precision HALF, TWO, THREE, FIVE, ZERO, PT95 integer MODE, NEXT, IC, NP parameter(HALF = 0.5D0, TWO = 2.0D0, THREE = 3.0D0) parameter(FIVE = 5.0D0, ZERO = 0.0D0, PT95 = 0.95D0) 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 DFMIN: 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('DFMIN',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(D1MACH(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