SUBROUTINE MINST(F,DXA,NC) * * Minimizer for onedimensional functions, especially for Newton * steps in multidimensional minimization * * Subroutine MINST searches for the minimum of a function fun(x) * of a single variable X in the region * XSTART .... XSTART+SDX*ST * for a given start value xstart and given values of SDX and ST. The * function fun(x) is evaluated up to 11 times. * * Increasing or decreasing steps are taken, until a function value * lower than fun(XSTART) is found. Then the minimum is searched for * near the minimum of a function, interpolating all previous * function values. The following convergence criteria are used: * * 1) the slope of the interpolating function is reduced to * 1/100 of the slope ate XSTART; * 2) the relative function variation of the last 3 function * evaluations if about 1.0E-6. * * Application for n-dimensional minimization: MINST may be used to * find the minimum along the direction of the Newton step. In this * case X (and XSTART) and ST are n-dimensional vectors, while DX * (and SDX) remains a scalar. For a Newton step ST the minimum * should be near XSTART+ST. * * Usage * ===== * initialization: * CALL MINSTI(XLIM1,XLIM2,IPRI) * XLIM1 = limit of step in negative direction; * XLIM2 = limit of step in positive direction, * assuming start at 0.0 * IPRI = print flag (0,1,2) * * The minimization starts at the start value XSTART (see below), * but the limits above correspond to the start value assumed as 0.0. * The sign of XLIM1 should be negative (or zero), the sign of * XLIM2 should be positove (or zero). * * XLIM2 -1.0 0.0 +1.0 XLIM2 * --------------|----------+---------+---------+--------|-------> X * forbidden | XSTART-ST XSTART XSTART-ST | forbidden * * * X = XSTART start value for X * ST= step size step ( * 10 F = FUN(X) evaluate function * CALL MINST(F,DX,NC) get factor DX * X = X + DX*ST modify X * IF(NC.LT.0) GOTO 10 test convergence *) * * F is the lowest function value and x the corresp. argument * * after final return: * CALL MINSTE(DFUN,DPAR,SECD) * DFUN = total function change * DPAR = total step * SECD = second derivative * * *) the value of NC has the following meaning: * NC = -1 additional function value required * = 0 conv. criterium satisfied * = 1 conv. criterium not satisfied * * SUBROUTINE CSPNK IS USED INTERNALLY. * PARAMETER (TAU=0.618033886,CASY=0.05) REAL X(11),Y(11),C(5,11) CHARACTER*4 TEXT COMMON/CMINST/DFUNST,DPARST,SECDST,IPRNST,XLIMN,XLIMP,NEVAST * save all ... SAVE DATA EPSMA/1.0E-6/,SLPMIN/0.01/ DATA IFG/1000/ * ... NC=-1 * initialization for NEVAST=0 IF(NEVAST.EQ.0) THEN XC=0.0 IC=0 END IF * insert new function value at correct place (INS) DO INS=1,NEVAST IF(XC.LT.X(INS)) THEN DO J=NEVAST,INS,-1 X(J+1)=X(J) Y(J+1)=Y(J) END DO GOTO 10 END IF END DO INS=NEVAST+1 10 NEVAST =NEVAST+1 * XC is current internal x-value X(INS)=XC Y(INS)=F * * find smallest function value * 20 IMIN=1 IBEG=1 DO I=2,NEVAST IF(X(I).EQ.0.0 ) IBEG=I IF(Y(I).LT.Y(IMIN)) IMIN=I IF(X(I)-X(I-1).LT.0.00001) THEN * remove double entry DO J=I+1,NEVAST X(J-1)=X(J) Y(J-1)=Y(J) END DO NEVAST=NEVAST-1 GOTO 20 END IF END DO * * check roundoff error region for last three points * IF(NEVAST.GE.3) THEN FQ=(F+FP+FPP)/3.0 RF=SQRT(0.5*((F-FQ)**2+(FP-FQ)**2+(FPP-FQ)**2))/(ABS(FQ)+1.0) AF=(ABS(F)+ABS(FP)+ABS(FPP))/3.0 IF(AF.NE.0.0) THEN AF=SQRT(0.5*((F-FQ)**2+(FP-FQ)**2+(FPP-FQ)**2))/AF END IF * test following statement RF=AF IF(RF.LE.EPSMA) GOTO 90 END IF * save last DX DXL=DX * XC is current value, determine now XN = next value =============== IF(NEVAST.EQ.1) THEN * second point (first point was at zero) ------------------------ XN=MIN(TAU,XLIMP) IF(XN.EQ.0.0) XN=MAX(-TAU,XLIMN) ELSE IF(NEVAST.EQ.2) THEN * third point IF(IMIN.EQ.INS) THEN * success - search in same direction with increased step DX=DXL/TAU IF(DX.GT.0.0) THEN * ... but to center if at limit IF(XC.EQ.XLIMP) DX=-DX/TAU ELSE IF(XC.EQ.XLIMN) DX=-DX/TAU END IF XN=XC+DX ELSE * failure - search in opposite direction IF(DX.GT.0.0) THEN XN=-1.0 IF(XLIMN.EQ.0.0) XN=0.5*XC ELSE IF(XC.EQ.XLIMN) THEN XN=0.5*XC ELSE * XN=MAX(XC+DXL/TAU,XLIMN) END IF END IF END IF ELSE IF(IMIN.EQ.NEVAST) THEN * monotonically decreasing function in positive direction IF(X(NEVAST).NE.XLIMP) THEN * ... not at upper limit DX=(X(NEVAST)-X(NEVAST-1))/TAU XN=MIN(X(NEVAST)+DX,XLIMP) ELSE * ... at upper limit - backward DX=(X(NEVAST-1)-X(NEVAST))*TAU*TAU*TAU XN=X(NEVAST)+DX END IF ELSE IF(IMIN.EQ.1) THEN * monotonically decreasing function in negative direction IF(X(1).NE.XLIMN) THEN * ... not at lower limit DX=(X(1)-X(2))/TAU XN=MAX(X(1)+DX,XLIMN) ELSE * ... at lower limit DX=(X(2)-X(1))*TAU*TAU*TAU XN=X(1)+DX END IF ELSE * inner minimum ------------------------------------------------- CALL CSPNK(X,Y,C,NEVAST) * slope ratio SLPR=ABS(C(2,IMIN)/C(2,IBEG)) IF(SLPR.GT.SLPMIN) IC=0 IF(IC.EQ.1) GOTO 90 * XT = approximate golden section value as alternative IF(X(IMIN+1)-X(IMIN).GT.X(IMIN)-X(IMIN-1)) THEN XT=X(IMIN)+(TAU-TAU*TAU)*(X(IMIN+1)-X(IMIN-1)) ELSE XT=X(IMIN)+(TAU-TAU*TAU)*(X(IMIN-1)-X(IMIN+1)) END IF * XN = minimum of interpolating function J=IMIN IF(C(2,J).GT.0.0) J=J-1 DIS=C(3,J)**2-3.0*C(2,J)*C(4,J) XDIS=SQRT(DIS)+ABS(C(3,J)) H=-C(2,J)/SIGN(XDIS,C(3,J)) XN=C(5,J)+H IF(SLPR.LT.0.5*SLPMIN.AND.NEVAST.GE.4) THEN * flag for final step IC=1 ELSE * add bias to step, to obtain good spacing BX=C(5,J+1)-C(5,J) XN=XN-0.25*BX*(H/BX-0.5)**3 END IF * step to XN DX=XN-XC IF(DX.EQ.0.0) GOTO 90 * calculate asymmetry of new point XN IF(XN.GT.X(IMIN)) THEN ASYM=AMIN1(C(5,IMIN+1)-XN,XN-C(5,IMIN)) + /(C(5,IMIN+1)-C(5,IMIN)) ASYL=AMIN1(XN-C(5,IMIN),C(5,IMIN)-C(5,IMIN-1)) + /(XN-C(5,IMIN-1)) ELSE ASYM=AMIN1(C(5,IMIN)-XN,XN-C(5,IMIN-1)) + /(C(5,IMIN)-C(5,IMIN-1)) ASYL=AMIN1(C(5,IMIN+1)-C(5,IMIN),C(5,IMIN)-XN) + /(C(5,IMIN+1)-XN) END IF IF(IC.EQ.0.AND.MIN(ASYM,ASYL).LT.CASY) THEN * minimum point replaced by golden section point XN=XT END IF END IF * next function evaluation ----------------------------------------- DX=XN-XC IF(NEVAST.GE.11) THEN NC=1 GOTO 94 END IF * save last point XL=XC * check limits and define next point XC ... IF(DX.NE.0.0) THEN XC=MIN(XC+DX,XLIMP) ELSE XC=MAX(XC+DX,XLIMN) END IF * ... and step to XC DX=XC-XL FPP=FP FP =F GOTO 100 * last return with preparation of final parameters ----------------- 90 NC=0 94 IF(IPRNST.NE.0) CALL CSPNK(X,Y,C,NEVAST) DX=X(IMIN)-XC F =Y(IMIN) * estimate second derivative SUM1=0.0 SUM2=0.0 DO I=1,NEVAST CC=(C(5,I)-C(5,IMIN))**2 SUM1=SUM1+CC*(C(1,1)-C(1,IMIN)) SUM2=SUM2+CC*CC END DO SECDST=0.0 IF(SUM2.NE.0.0) SECDST=2.0*SUM1/SUM2 * set common values DFUNST=C(1,IBEG)-C(1,IMIN) DPARST=C(5,IMIN) * printout --------------------------------------------------------- IF(IPRNST.NE.0) THEN WRITE(*,101) NC,NEVAST,C(5,IMIN),C(1,IMIN), + XLIMN,XLIMP,DFUNST IF(IPRNST.EQ.2) THEN WRITE(*,102) DO I=1,NEVAST TEXT=' ' IF(I.EQ.IMIN) TEXT='Min' WRITE(*,103) TEXT,C(5,I),C(1,I),C(2,I),C(3,I) END DO WRITE(*,103) END IF END IF 100 DXA=DX RETURN 101 FORMAT('0Minst',7X,'Return code',I3,' after',I3,' Function', 1 ' evaluations'/ 2 11X,' Minimum at ',G15.6,' Function F =',G15.6/ 3 11X,' Limits ',2G15.6,' DF =',G15.6/) 102 FORMAT(13X,'Table of function values and derivatives'/ 1 15X,5H X ,10X,5H Fun ,10X,5HFun' ,10X,5HFun''/) 103 FORMAT(7X,A4,1X,G13.6,3G15.6) END SUBROUTINE MINSTI(XLIM1,XLIM2,IPRI) * initialize MINST: limits and print flag COMMON/CMINST/DFUNST,DPARST,SECDST,IPRNST,XLIMN,XLIMP,NEVAST SAVE /CMINST/ * ... XLIMN=-ABS(XLIM1) XLIMP=+ABS(XLIM2) IF(XLIMP.EQ.0.0.AND.XLIMN.EQ.0.0) THEN XLIMP=+5.0 XLIMN=-5.0 END IF IPRNST=IPRI NEVAST=0 END SUBROUTINE MINSTE(DFUN,DPAR,SECD) * return final parameters COMMON/CMINST/DFUNST,DPARST,SECDST,IPRNST,XLIMN,XLIMP,NEVAST SAVE /CMINST/ * ... DFUN=DFUNST DPAR=DPARST SECD=SECDST END