share/lbfgs/mcsrch.f

   1 C    ------------------------------------------------------------------
   2 C
   3 C     **************************
   4 C     LINE SEARCH ROUTINE MCSRCH
   5 C     **************************
   6 C
   7       SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL,MAXFEV,INFO,NFEV,WA)
   8       INTEGER N,MAXFEV,INFO,NFEV
   9       DOUBLE PRECISION F,STP,FTOL,GTOL,XTOL,STPMIN,STPMAX
  10       DOUBLE PRECISION X(N),G(N),S(N),WA(N)
  11       COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
  12       SAVE
  13 C
  14 C                     SUBROUTINE MCSRCH
  15 C
  16 C     A slight modification of the subroutine CSRCH of More' and Thuente.
  17 C     The changes are to allow reverse communication, and do not affect
  18 C     the performance of the routine.
  19 C
  20 C     THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES
  21 C     A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION.
  22 C
  23 C     AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF
  24 C     UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF
  25 C     UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A
  26 C     MINIMIZER OF THE MODIFIED FUNCTION
  27 C
  28 C          F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
  29 C
  30 C     IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION
  31 C     HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE,
  32 C     THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT
  33 C     CONTAINS A MINIMIZER OF F(X+STP*S).
  34 C
  35 C     THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES
  36 C     THE SUFFICIENT DECREASE CONDITION
  37 C
  38 C           F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S),
  39 C
  40 C     AND THE CURVATURE CONDITION
  41 C
  42 C           ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S).
  43 C
  44 C     IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION
  45 C     IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES
  46 C     BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH
  47 C     CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING
  48 C     ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY
  49 C     SATISFIES THE SUFFICIENT DECREASE CONDITION.
  50 C
  51 C     THE SUBROUTINE STATEMENT IS
  52 C
  53 C        SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL, MAXFEV,INFO,NFEV,WA)
  54 C     WHERE
  55 C
  56 C       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
  57 C         OF VARIABLES.
  58 C
  59 C       X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
  60 C         BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS
  61 C         X + STP*S.
  62 C
  63 C       F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F
  64 C         AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S.
  65 C
  66 C       G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
  67 C         GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT
  68 C         OF F AT X + STP*S.
  69 C
  70 C       S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE
  71 C         SEARCH DIRECTION.
  72 C
  73 C       STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN
  74 C         INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT
  75 C         STP CONTAINS THE FINAL ESTIMATE.
  76 C
  77 C       FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. (In this reverse
  78 C         communication implementation GTOL is defined in a COMMON
  79 C         statement.) TERMINATION OCCURS WHEN THE SUFFICIENT DECREASE
  80 C         CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
  81 C         SATISFIED.
  82 C
  83 C       XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
  84 C         WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
  85 C         IS AT MOST XTOL.
  86 C
  87 C       STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH
  88 C         SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. (In this reverse
  89 C         communication implementatin they are defined in a COMMON
  90 C         statement).
  91 C
  92 C       MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
  93 C         OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST
  94 C         MAXFEV BY THE END OF AN ITERATION.
  95 C
  96 C       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
  97 C
  98 C         INFO = 0  IMPROPER INPUT PARAMETERS.
  99 C
 100 C         INFO =-1  A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT.
 101 C
 102 C         INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THE
 103 C                   DIRECTIONAL DERIVATIVE CONDITION HOLD.
 104 C
 105 C         INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
 106 C                   IS AT MOST XTOL.
 107 C
 108 C         INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
 109 C
 110 C         INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN.
 111 C
 112 C         INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX.
 113 C
 114 C         INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.
 115 C                   THERE MAY NOT BE A STEP WHICH SATISFIES THE
 116 C                   SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
 117 C                   TOLERANCES MAY BE TOO SMALL.
 118 C
 119 C       NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
 120 C         CALLS TO FCN.
 121 C
 122 C       WA IS A WORK ARRAY OF LENGTH N.
 123 C
 124 C     SUBPROGRAMS CALLED
 125 C
 126 C       MCSTEP
 127 C
 128 C       FORTRAN-SUPPLIED...ABS,MAX,MIN
 129 C
 130 C     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
 131 C     JORGE J. MORE', DAVID J. THUENTE
 132 C
 133 C     **********
 134       INTEGER INFOC,J
 135       LOGICAL BRACKT,STAGE1
 136       DOUBLE PRECISION DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,
 137      *       FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY,
 138      *       STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,ZERO
 139       DATA P5,P66,XTRAPF,ZERO /0.5D0,0.66D0,4.0D0,0.0D0/
 140       IF(INFO.EQ.-1) GO TO 45
 141       INFOC = 1
 142 C
 143 C     CHECK THE INPUT PARAMETERS FOR ERRORS.
 144 C
 145       IF (N .LE. 0 .OR. STP .LE. ZERO .OR. FTOL .LT. ZERO .OR.
 146      *    GTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. STPMIN .LT. ZERO
 147      *    .OR. STPMAX .LT. STPMIN .OR. MAXFEV .LE. 0) RETURN
 148 C
 149 C     COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
 150 C     AND CHECK THAT S IS A DESCENT DIRECTION.
 151 C
 152       DGINIT = ZERO
 153       DO 10 J = 1, N
 154          DGINIT = DGINIT + G(J)*S(J)
 155    10    CONTINUE
 156       IF (DGINIT .GE. ZERO) then
 157          write(LP,15)
 158    15    FORMAT(/'  THE SEARCH DIRECTION IS NOT A DESCENT DIRECTION')
 159          RETURN
 160          ENDIF
 161 C
 162 C     INITIALIZE LOCAL VARIABLES.
 163 C
 164       BRACKT = .FALSE.
 165       STAGE1 = .TRUE.
 166       NFEV = 0
 167       FINIT = F
 168       DGTEST = FTOL*DGINIT
 169       WIDTH = STPMAX - STPMIN
 170       WIDTH1 = WIDTH/P5
 171       DO 20 J = 1, N
 172          WA(J) = X(J)
 173    20    CONTINUE
 174 C
 175 C     THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
 176 C     FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
 177 C     THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
 178 C     FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
 179 C     THE INTERVAL OF UNCERTAINTY.
 180 C     THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
 181 C     FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
 182 C
 183       STX = ZERO
 184       FX = FINIT
 185       DGX = DGINIT
 186       STY = ZERO
 187       FY = FINIT
 188       DGY = DGINIT
 189 C
 190 C     START OF ITERATION.
 191 C
 192    30 CONTINUE
 193 C
 194 C        SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
 195 C        TO THE PRESENT INTERVAL OF UNCERTAINTY.
 196 C
 197          IF (BRACKT) THEN
 198             STMIN = MIN(STX,STY)
 199             STMAX = MAX(STX,STY)
 200          ELSE
 201             STMIN = STX
 202             STMAX = STP + XTRAPF*(STP - STX)
 203             END IF
 204 C
 205 C        FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
 206 C
 207          STP = MAX(STP,STPMIN)
 208          STP = MIN(STP,STPMAX)
 209 C
 210 C        IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
 211 C        STP BE THE LOWEST POINT OBTAINED SO FAR.
 212 C
 213          IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
 214      *      .OR. NFEV .GE. MAXFEV-1 .OR. INFOC .EQ. 0
 215      *      .OR. (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX)) STP = STX
 216 C
 217 C        EVALUATE THE FUNCTION AND GRADIENT AT STP
 218 C        AND COMPUTE THE DIRECTIONAL DERIVATIVE.
 219 C        We return to main program to obtain F and G.
 220 C
 221          DO 40 J = 1, N
 222             X(J) = WA(J) + STP*S(J)
 223    40       CONTINUE
 224          INFO=-1
 225          RETURN
 226 C
 227    45    INFO=0
 228          NFEV = NFEV + 1
 229          DG = ZERO
 230          DO 50 J = 1, N
 231             DG = DG + G(J)*S(J)
 232    50       CONTINUE
 233          FTEST1 = FINIT + STP*DGTEST
 234 C
 235 C        TEST FOR CONVERGENCE.
 236 C
 237          IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
 238      *      .OR. INFOC .EQ. 0) INFO = 6
 239          IF (STP .EQ. STPMAX .AND.
 240      *       F .LE. FTEST1 .AND. DG .LE. DGTEST) INFO = 5
 241          IF (STP .EQ. STPMIN .AND.
 242      *       (F .GT. FTEST1 .OR. DG .GE. DGTEST)) INFO = 4
 243          IF (NFEV .GE. MAXFEV) INFO = 3
 244          IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2
 245          IF (F .LE. FTEST1 .AND. DABS(DG) .LE. GTOL*(-DGINIT)) INFO = 1
 246 C
 247 C        CHECK FOR TERMINATION.
 248 C
 249          IF (INFO .NE. 0) RETURN
 250 C
 251 C        IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
 252 C        FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
 253 C
 254          IF (STAGE1 .AND. F .LE. FTEST1 .AND.
 255      *       DG .GE. MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE.
 256 C
 257 C        A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
 258 C        WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
 259 C        FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
 260 C        DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
 261 C        OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
 262 C
 263          IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN
 264 C
 265 C           DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
 266 C
 267             FM = F - STP*DGTEST
 268             FXM = FX - STX*DGTEST
 269             FYM = FY - STY*DGTEST
 270             DGM = DG - DGTEST
 271             DGXM = DGX - DGTEST
 272             DGYM = DGY - DGTEST
 273 C
 274 C           CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
 275 C           AND TO COMPUTE THE NEW STEP.
 276 C
 277             CALL MCSTEP(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,
 278      *                 BRACKT,STMIN,STMAX,INFOC)
 279 C
 280 C           RESET THE FUNCTION AND GRADIENT VALUES FOR F.
 281 C
 282             FX = FXM + STX*DGTEST
 283             FY = FYM + STY*DGTEST
 284             DGX = DGXM + DGTEST
 285             DGY = DGYM + DGTEST
 286          ELSE
 287 C
 288 C           CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
 289 C           AND TO COMPUTE THE NEW STEP.
 290 C
 291             CALL MCSTEP(STX,FX,DGX,STY,FY,DGY,STP,F,DG,
 292      *                 BRACKT,STMIN,STMAX,INFOC)
 293             END IF
 294 C
 295 C        FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
 296 C        INTERVAL OF UNCERTAINTY.
 297 C
 298          IF (BRACKT) THEN
 299             IF (DABS(STY-STX) .GE. P66*WIDTH1)
 300      *         STP = STX + P5*(STY - STX)
 301             WIDTH1 = WIDTH
 302             WIDTH = DABS(STY-STX)
 303             END IF
 304 C
 305 C        END OF ITERATION.
 306 C
 307          GO TO 30
 308 C
 309 C     LAST LINE OF SUBROUTINE MCSRCH.
 310 C
 311       END