| [0b990d] | 1 | c Below are the contents of the original file that was used to generate
|
|---|
| 2 | c mcsearch.h and mcsearch.cc. Only the mcsrch and mcstep routines are
|
|---|
| 3 | c used. This file is not compiled or otherwise used.
|
|---|
| 4 | C ----------------------------------------------------------------------
|
|---|
| 5 | C This file contains the LBFGS algorithm and supporting routines
|
|---|
| 6 | C
|
|---|
| 7 | C ****************
|
|---|
| 8 | C LBFGS SUBROUTINE
|
|---|
| 9 | C ****************
|
|---|
| 10 | C
|
|---|
| 11 | SUBROUTINE LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG)
|
|---|
| 12 | C
|
|---|
| 13 | INTEGER N,M,IPRINT(2),IFLAG
|
|---|
| 14 | DOUBLE PRECISION X(N),G(N),DIAG(N),W(N*(2*M+1)+2*M)
|
|---|
| 15 | DOUBLE PRECISION F,EPS,XTOL
|
|---|
| 16 | INTEGER DIAGCO
|
|---|
| 17 | C
|
|---|
| 18 | C LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION
|
|---|
| 19 | C JORGE NOCEDAL
|
|---|
| 20 | C *** July 1990 ***
|
|---|
| 21 | C
|
|---|
| 22 | C
|
|---|
| 23 | C This subroutine solves the unconstrained minimization problem
|
|---|
| 24 | C
|
|---|
| 25 | C min F(x), x= (x1,x2,...,xN),
|
|---|
| 26 | C
|
|---|
| 27 | C using the limited memory BFGS method. The routine is especially
|
|---|
| 28 | C effective on problems involving a large number of variables. In
|
|---|
| 29 | C a typical iteration of this method an approximation Hk to the
|
|---|
| 30 | C inverse of the Hessian is obtained by applying M BFGS updates to
|
|---|
| 31 | C a diagonal matrix Hk0, using information from the previous M steps.
|
|---|
| 32 | C The user specifies the number M, which determines the amount of
|
|---|
| 33 | C storage required by the routine. The user may also provide the
|
|---|
| 34 | C diagonal matrices Hk0 if not satisfied with the default choice.
|
|---|
| 35 | C The algorithm is described in "On the limited memory BFGS method
|
|---|
| 36 | C for large scale optimization", by D. Liu and J. Nocedal,
|
|---|
| 37 | C Mathematical Programming B 45 (1989) 503-528.
|
|---|
| 38 | C
|
|---|
| 39 | C The user is required to calculate the function value F and its
|
|---|
| 40 | C gradient G. In order to allow the user complete control over
|
|---|
| 41 | C these computations, reverse communication is used. The routine
|
|---|
| 42 | C must be called repeatedly under the control of the parameter
|
|---|
| 43 | C IFLAG.
|
|---|
| 44 | C
|
|---|
| 45 | C The steplength is determined at each iteration by means of the
|
|---|
| 46 | C line search routine MCVSRCH, which is a slight modification of
|
|---|
| 47 | C the routine CSRCH written by More' and Thuente.
|
|---|
| 48 | C
|
|---|
| 49 | C The calling statement is
|
|---|
| 50 | C
|
|---|
| 51 | C CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG)
|
|---|
| 52 | C
|
|---|
| 53 | C where
|
|---|
| 54 | C
|
|---|
| 55 | C N is an INTEGER variable that must be set by the user to the
|
|---|
| 56 | C number of variables. It is not altered by the routine.
|
|---|
| 57 | C Restriction: N>0.
|
|---|
| 58 | C
|
|---|
| 59 | C M is an INTEGER variable that must be set by the user to
|
|---|
| 60 | C the number of corrections used in the BFGS update. It
|
|---|
| 61 | C is not altered by the routine. Values of M less than 3 are
|
|---|
| 62 | C not recommended; large values of M will result in excessive
|
|---|
| 63 | C computing time. 3<= M <=7 is recommended. Restriction: M>0.
|
|---|
| 64 | C
|
|---|
| 65 | C X is a DOUBLE PRECISION array of length N. On initial entry
|
|---|
| 66 | C it must be set by the user to the values of the initial
|
|---|
| 67 | C estimate of the solution vector. On exit with IFLAG=0, it
|
|---|
| 68 | C contains the values of the variables at the best point
|
|---|
| 69 | C found (usually a solution).
|
|---|
| 70 | C
|
|---|
| 71 | C F is a DOUBLE PRECISION variable. Before initial entry and on
|
|---|
| 72 | C a re-entry with IFLAG=1, it must be set by the user to
|
|---|
| 73 | C contain the value of the function F at the point X.
|
|---|
| 74 | C
|
|---|
| 75 | C G is a DOUBLE PRECISION array of length N. Before initial
|
|---|
| 76 | C entry and on a re-entry with IFLAG=1, it must be set by
|
|---|
| 77 | C the user to contain the components of the gradient G at
|
|---|
| 78 | C the point X.
|
|---|
| 79 | C
|
|---|
| 80 | C DIAGCO is a LOGICAL variable that must be set to .TRUE. if the
|
|---|
| 81 | C user wishes to provide the diagonal matrix Hk0 at each
|
|---|
| 82 | C iteration. Otherwise it should be set to .FALSE., in which
|
|---|
| 83 | C case LBFGS will use a default value described below. If
|
|---|
| 84 | C DIAGCO is set to .TRUE. the routine will return at each
|
|---|
| 85 | C iteration of the algorithm with IFLAG=2, and the diagonal
|
|---|
| 86 | C matrix Hk0 must be provided in the array DIAG.
|
|---|
| 87 | C
|
|---|
| 88 | C
|
|---|
| 89 | C DIAG is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE.,
|
|---|
| 90 | C then on initial entry or on re-entry with IFLAG=2, DIAG
|
|---|
| 91 | C it must be set by the user to contain the values of the
|
|---|
| 92 | C diagonal matrix Hk0. Restriction: all elements of DIAG
|
|---|
| 93 | C must be positive.
|
|---|
| 94 | C
|
|---|
| 95 | C IPRINT is an INTEGER array of length two which must be set by the
|
|---|
| 96 | C user.
|
|---|
| 97 | C
|
|---|
| 98 | C IPRINT(1) specifies the frequency of the output:
|
|---|
| 99 | C IPRINT(1) < 0 : no output is generated,
|
|---|
| 100 | C IPRINT(1) = 0 : output only at first and last iteration,
|
|---|
| 101 | C IPRINT(1) > 0 : output every IPRINT(1) iterations.
|
|---|
| 102 | C
|
|---|
| 103 | C IPRINT(2) specifies the type of output generated:
|
|---|
| 104 | C IPRINT(2) = 0 : iteration count, number of function
|
|---|
| 105 | C evaluations, function value, norm of the
|
|---|
| 106 | C gradient, and steplength,
|
|---|
| 107 | C IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of
|
|---|
| 108 | C variables and gradient vector at the
|
|---|
| 109 | C initial point,
|
|---|
| 110 | C IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of
|
|---|
| 111 | C variables,
|
|---|
| 112 | C IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector.
|
|---|
| 113 | C
|
|---|
| 114 | C
|
|---|
| 115 | C EPS is a positive DOUBLE PRECISION variable that must be set by
|
|---|
| 116 | C the user, and determines the accuracy with which the solution
|
|---|
| 117 | C is to be found. The subroutine terminates when
|
|---|
| 118 | C
|
|---|
| 119 | C ||G|| < EPS max(1,||X||),
|
|---|
| 120 | C
|
|---|
| 121 | C where ||.|| denotes the Euclidean norm.
|
|---|
| 122 | C
|
|---|
| 123 | C XTOL is a positive DOUBLE PRECISION variable that must be set by
|
|---|
| 124 | C the user to an estimate of the machine precision (e.g.
|
|---|
| 125 | C 10**(-16) on a SUN station 3/60). The line search routine will
|
|---|
| 126 | C terminate if the relative width of the interval of uncertainty
|
|---|
| 127 | C is less than XTOL.
|
|---|
| 128 | C
|
|---|
| 129 | C W is a DOUBLE PRECISION array of length N(2M+1)+2M used as
|
|---|
| 130 | C workspace for LBFGS. This array must not be altered by the
|
|---|
| 131 | C user.
|
|---|
| 132 | C
|
|---|
| 133 | C IFLAG is an INTEGER variable that must be set to 0 on initial entry
|
|---|
| 134 | C to the subroutine. A return with IFLAG<0 indicates an error,
|
|---|
| 135 | C and IFLAG=0 indicates that the routine has terminated without
|
|---|
| 136 | C detecting errors. On a return with IFLAG=1, the user must
|
|---|
| 137 | C evaluate the function F and gradient G. On a return with
|
|---|
| 138 | C IFLAG=2, the user must provide the diagonal matrix Hk0.
|
|---|
| 139 | C
|
|---|
| 140 | C The following negative values of IFLAG, detecting an error,
|
|---|
| 141 | C are possible:
|
|---|
| 142 | C
|
|---|
| 143 | C IFLAG=-1 The line search routine MCSRCH failed. The
|
|---|
| 144 | C parameter INFO provides more detailed information
|
|---|
| 145 | C (see also the documentation of MCSRCH):
|
|---|
| 146 | C
|
|---|
| 147 | C INFO = 0 IMPROPER INPUT PARAMETERS.
|
|---|
| 148 | C
|
|---|
| 149 | C INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF
|
|---|
| 150 | C UNCERTAINTY IS AT MOST XTOL.
|
|---|
| 151 | C
|
|---|
| 152 | C INFO = 3 MORE THAN 20 FUNCTION EVALUATIONS WERE
|
|---|
| 153 | C REQUIRED AT THE PRESENT ITERATION.
|
|---|
| 154 | C
|
|---|
| 155 | C INFO = 4 THE STEP IS TOO SMALL.
|
|---|
| 156 | C
|
|---|
| 157 | C INFO = 5 THE STEP IS TOO LARGE.
|
|---|
| 158 | C
|
|---|
| 159 | C INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.
|
|---|
| 160 | C THERE MAY NOT BE A STEP WHICH SATISFIES
|
|---|
| 161 | C THE SUFFICIENT DECREASE AND CURVATURE
|
|---|
| 162 | C CONDITIONS. TOLERANCES MAY BE TOO SMALL.
|
|---|
| 163 | C
|
|---|
| 164 | C
|
|---|
| 165 | C IFLAG=-2 The i-th diagonal element of the diagonal inverse
|
|---|
| 166 | C Hessian approximation, given in DIAG, is not
|
|---|
| 167 | C positive.
|
|---|
| 168 | C
|
|---|
| 169 | C IFLAG=-3 Improper input parameters for LBFGS (N or M are
|
|---|
| 170 | C not positive).
|
|---|
| 171 | C
|
|---|
| 172 | C
|
|---|
| 173 | C
|
|---|
| 174 | C ON THE DRIVER:
|
|---|
| 175 | C
|
|---|
| 176 | C The program that calls LBFGS must contain the declaration:
|
|---|
| 177 | C
|
|---|
| 178 | C EXTERNAL LB2
|
|---|
| 179 | C
|
|---|
| 180 | C LB2 is a BLOCK DATA that defines the default values of several
|
|---|
| 181 | C parameters described in the COMMON section.
|
|---|
| 182 | C
|
|---|
| 183 | C
|
|---|
| 184 | C
|
|---|
| 185 | C COMMON:
|
|---|
| 186 | C
|
|---|
| 187 | C The subroutine contains one common area, which the user may wish to
|
|---|
| 188 | C reference:
|
|---|
| 189 | C
|
|---|
| 190 | COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
|
|---|
| 191 | C
|
|---|
| 192 | C MP is an INTEGER variable with default value 6. It is used as the
|
|---|
| 193 | C unit number for the printing of the monitoring information
|
|---|
| 194 | C controlled by IPRINT.
|
|---|
| 195 | C
|
|---|
| 196 | C LP is an INTEGER variable with default value 6. It is used as the
|
|---|
| 197 | C unit number for the printing of error messages. This printing
|
|---|
| 198 | C may be suppressed by setting LP to a non-positive value.
|
|---|
| 199 | C
|
|---|
| 200 | C GTOL is a DOUBLE PRECISION variable with default value 0.9, which
|
|---|
| 201 | C controls the accuracy of the line search routine MCSRCH. If the
|
|---|
| 202 | C function and gradient evaluations are inexpensive with respect
|
|---|
| 203 | C to the cost of the iteration (which is sometimes the case when
|
|---|
| 204 | C solving very large problems) it may be advantageous to set GTOL
|
|---|
| 205 | C to a small value. A typical small value is 0.1. Restriction:
|
|---|
| 206 | C GTOL should be greater than 1.D-04.
|
|---|
| 207 | C
|
|---|
| 208 | C STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which
|
|---|
| 209 | C specify lower and uper bounds for the step in the line search.
|
|---|
| 210 | C Their default values are 1.D-20 and 1.D+20, respectively. These
|
|---|
| 211 | C values need not be modified unless the exponents are too large
|
|---|
| 212 | C for the machine being used, or unless the problem is extremely
|
|---|
| 213 | C badly scaled (in which case the exponents should be increased).
|
|---|
| 214 | C
|
|---|
| 215 | C
|
|---|
| 216 | C MACHINE DEPENDENCIES
|
|---|
| 217 | C
|
|---|
| 218 | C The only variables that are machine-dependent are XTOL,
|
|---|
| 219 | C STPMIN and STPMAX.
|
|---|
| 220 | C
|
|---|
| 221 | C
|
|---|
| 222 | C GENERAL INFORMATION
|
|---|
| 223 | C
|
|---|
| 224 | C Other routines called directly: DAXPY, DDOT, LB1, MCSRCH
|
|---|
| 225 | C
|
|---|
| 226 | C Input/Output : No input; diagnostic messages on unit MP and
|
|---|
| 227 | C error messages on unit LP.
|
|---|
| 228 | C
|
|---|
| 229 | C
|
|---|
| 230 | C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|---|
| 231 | C
|
|---|
| 232 | DOUBLE PRECISION GTOL,ONE,ZERO,GNORM,DDOT,STP1,FTOL,STPMIN,
|
|---|
| 233 | . STPMAX,STP,YS,YY,SQ,YR,BETA,XNORM
|
|---|
| 234 | INTEGER MP,LP,ITER,NFUN,POINT,ISPT,IYPT,MAXFEV,INFO,
|
|---|
| 235 | . BOUND,NPT,CP,I,NFEV,INMC,IYCN,ISCN
|
|---|
| 236 | LOGICAL FINISH
|
|---|
| 237 | C
|
|---|
| 238 | SAVE
|
|---|
| 239 | DATA ONE,ZERO/1.0D+0,0.0D+0/
|
|---|
| 240 | C
|
|---|
| 241 | C INITIALIZE
|
|---|
| 242 | C ----------
|
|---|
| 243 | C
|
|---|
| 244 | IF(IFLAG.EQ.0) GO TO 10
|
|---|
| 245 | GO TO (172,100) IFLAG
|
|---|
| 246 | 10 ITER= 0
|
|---|
| 247 | IF(N.LE.0.OR.M.LE.0) GO TO 196
|
|---|
| 248 | IF(GTOL.LE.1.D-04) THEN
|
|---|
| 249 | IF(LP.GT.0) WRITE(LP,245)
|
|---|
| 250 | GTOL=9.D-01
|
|---|
| 251 | ENDIF
|
|---|
| 252 | NFUN= 1
|
|---|
| 253 | POINT= 0
|
|---|
| 254 | FINISH= .FALSE.
|
|---|
| 255 | IF(DIAGCO.NE.0) THEN
|
|---|
| 256 | DO 30 I=1,N
|
|---|
| 257 | 30 IF (DIAG(I).LE.ZERO) GO TO 195
|
|---|
| 258 | ELSE
|
|---|
| 259 | DO 40 I=1,N
|
|---|
| 260 | 40 DIAG(I)= 1.0D0
|
|---|
| 261 | ENDIF
|
|---|
| 262 | C
|
|---|
| 263 | C THE WORK VECTOR W IS DIVIDED AS FOLLOWS:
|
|---|
| 264 | C ---------------------------------------
|
|---|
| 265 | C THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT AND
|
|---|
| 266 | C OTHER TEMPORARY INFORMATION.
|
|---|
| 267 | C LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO.
|
|---|
| 268 | C LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USED
|
|---|
| 269 | C IN THE FORMULA THAT COMPUTES H*G.
|
|---|
| 270 | C LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCH
|
|---|
| 271 | C STEPS.
|
|---|
| 272 | C LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST M
|
|---|
| 273 | C GRADIENT DIFFERENCES.
|
|---|
| 274 | C
|
|---|
| 275 | C THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN A
|
|---|
| 276 | C CIRCULAR ORDER CONTROLLED BY THE PARAMETER POINT.
|
|---|
| 277 | C
|
|---|
| 278 | ISPT= N+2*M
|
|---|
| 279 | IYPT= ISPT+N*M
|
|---|
| 280 | DO 50 I=1,N
|
|---|
| 281 | 50 W(ISPT+I)= -G(I)*DIAG(I)
|
|---|
| 282 | GNORM= DSQRT(DDOT(N,G,1,G,1))
|
|---|
| 283 | STP1= ONE/GNORM
|
|---|
| 284 | C
|
|---|
| 285 | C PARAMETERS FOR LINE SEARCH ROUTINE
|
|---|
| 286 | C
|
|---|
| 287 | FTOL= 1.0D-4
|
|---|
| 288 | MAXFEV= 20
|
|---|
| 289 | C
|
|---|
| 290 | IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,
|
|---|
| 291 | * GNORM,N,M,X,F,G,STP,FINISH)
|
|---|
| 292 | C
|
|---|
| 293 | C --------------------
|
|---|
| 294 | C MAIN ITERATION LOOP
|
|---|
| 295 | C --------------------
|
|---|
| 296 | C
|
|---|
| 297 | 80 ITER= ITER+1
|
|---|
| 298 | INFO=0
|
|---|
| 299 | BOUND=ITER-1
|
|---|
| 300 | IF(ITER.EQ.1) GO TO 165
|
|---|
| 301 | IF (ITER .GT. M)BOUND=M
|
|---|
| 302 | C
|
|---|
| 303 | YS= DDOT(N,W(IYPT+NPT+1),1,W(ISPT+NPT+1),1)
|
|---|
| 304 | IF(DIAGCO.EQ.0) THEN
|
|---|
| 305 | YY= DDOT(N,W(IYPT+NPT+1),1,W(IYPT+NPT+1),1)
|
|---|
| 306 | DO 90 I=1,N
|
|---|
| 307 | 90 DIAG(I)= YS/YY
|
|---|
| 308 | ELSE
|
|---|
| 309 | IFLAG=2
|
|---|
| 310 | RETURN
|
|---|
| 311 | ENDIF
|
|---|
| 312 | 100 CONTINUE
|
|---|
| 313 | IF(DIAGCO.NE.0) THEN
|
|---|
| 314 | DO 110 I=1,N
|
|---|
| 315 | 110 IF (DIAG(I).LE.ZERO) GO TO 195
|
|---|
| 316 | ENDIF
|
|---|
| 317 | C
|
|---|
| 318 | C COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980,
|
|---|
| 319 | C "Updating quasi-Newton matrices with limited storage",
|
|---|
| 320 | C Mathematics of Computation, Vol.24, No.151, pp. 773-782.
|
|---|
| 321 | C ---------------------------------------------------------
|
|---|
| 322 | C
|
|---|
| 323 | CP= POINT
|
|---|
| 324 | IF (POINT.EQ.0) CP=M
|
|---|
| 325 | W(N+CP)= ONE/YS
|
|---|
| 326 | DO 112 I=1,N
|
|---|
| 327 | 112 W(I)= -G(I)
|
|---|
| 328 | CP= POINT
|
|---|
| 329 | DO 125 I= 1,BOUND
|
|---|
| 330 | CP=CP-1
|
|---|
| 331 | IF (CP.EQ. -1)CP=M-1
|
|---|
| 332 | SQ= DDOT(N,W(ISPT+CP*N+1),1,W,1)
|
|---|
| 333 | INMC=N+M+CP+1
|
|---|
| 334 | IYCN=IYPT+CP*N
|
|---|
| 335 | W(INMC)= W(N+CP+1)*SQ
|
|---|
| 336 | CALL DAXPY(N,-W(INMC),W(IYCN+1),1,W,1)
|
|---|
| 337 | 125 CONTINUE
|
|---|
| 338 | C
|
|---|
| 339 | DO 130 I=1,N
|
|---|
| 340 | 130 W(I)=DIAG(I)*W(I)
|
|---|
| 341 | C
|
|---|
| 342 | DO 145 I=1,BOUND
|
|---|
| 343 | YR= DDOT(N,W(IYPT+CP*N+1),1,W,1)
|
|---|
| 344 | BETA= W(N+CP+1)*YR
|
|---|
| 345 | INMC=N+M+CP+1
|
|---|
| 346 | BETA= W(INMC)-BETA
|
|---|
| 347 | ISCN=ISPT+CP*N
|
|---|
| 348 | CALL DAXPY(N,BETA,W(ISCN+1),1,W,1)
|
|---|
| 349 | CP=CP+1
|
|---|
| 350 | IF (CP.EQ.M)CP=0
|
|---|
| 351 | 145 CONTINUE
|
|---|
| 352 | C
|
|---|
| 353 | C STORE THE NEW SEARCH DIRECTION
|
|---|
| 354 | C ------------------------------
|
|---|
| 355 | C
|
|---|
| 356 | DO 160 I=1,N
|
|---|
| 357 | 160 W(ISPT+POINT*N+I)= W(I)
|
|---|
| 358 | C
|
|---|
| 359 | C OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION
|
|---|
| 360 | C BY USING THE LINE SEARCH ROUTINE MCSRCH
|
|---|
| 361 | C ----------------------------------------------------
|
|---|
| 362 | 165 NFEV=0
|
|---|
| 363 | STP=ONE
|
|---|
| 364 | IF (ITER.EQ.1) STP=STP1
|
|---|
| 365 | DO 170 I=1,N
|
|---|
| 366 | 170 W(I)=G(I)
|
|---|
| 367 | 172 CONTINUE
|
|---|
| 368 | CALL MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,
|
|---|
| 369 | * XTOL,MAXFEV,INFO,NFEV,DIAG)
|
|---|
| 370 | IF (INFO .EQ. -1) THEN
|
|---|
| 371 | IFLAG=1
|
|---|
| 372 | RETURN
|
|---|
| 373 | ENDIF
|
|---|
| 374 | IF (INFO .NE. 1) GO TO 190
|
|---|
| 375 | NFUN= NFUN + NFEV
|
|---|
| 376 | C
|
|---|
| 377 | C COMPUTE THE NEW STEP AND GRADIENT CHANGE
|
|---|
| 378 | C -----------------------------------------
|
|---|
| 379 | C
|
|---|
| 380 | NPT=POINT*N
|
|---|
| 381 | DO 175 I=1,N
|
|---|
| 382 | W(ISPT+NPT+I)= STP*W(ISPT+NPT+I)
|
|---|
| 383 | 175 W(IYPT+NPT+I)= G(I)-W(I)
|
|---|
| 384 | POINT=POINT+1
|
|---|
| 385 | IF (POINT.EQ.M)POINT=0
|
|---|
| 386 | C
|
|---|
| 387 | C TERMINATION TEST
|
|---|
| 388 | C ----------------
|
|---|
| 389 | C
|
|---|
| 390 | GNORM= DSQRT(DDOT(N,G,1,G,1))
|
|---|
| 391 | XNORM= DSQRT(DDOT(N,X,1,X,1))
|
|---|
| 392 | XNORM= DMAX1(1.0D0,XNORM)
|
|---|
| 393 | IF (GNORM/XNORM .LE. EPS) FINISH=.TRUE.
|
|---|
| 394 | C
|
|---|
| 395 | IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,
|
|---|
| 396 | * GNORM,N,M,X,F,G,STP,FINISH)
|
|---|
| 397 | IF (FINISH) THEN
|
|---|
| 398 | IFLAG=0
|
|---|
| 399 | RETURN
|
|---|
| 400 | ENDIF
|
|---|
| 401 | GO TO 80
|
|---|
| 402 | C
|
|---|
| 403 | C ------------------------------------------------------------
|
|---|
| 404 | C END OF MAIN ITERATION LOOP. ERROR EXITS.
|
|---|
| 405 | C ------------------------------------------------------------
|
|---|
| 406 | C
|
|---|
| 407 | 190 IFLAG=-1
|
|---|
| 408 | IF(LP.GT.0) WRITE(LP,200) INFO
|
|---|
| 409 | RETURN
|
|---|
| 410 | 195 IFLAG=-2
|
|---|
| 411 | IF(LP.GT.0) WRITE(LP,235) I
|
|---|
| 412 | RETURN
|
|---|
| 413 | 196 IFLAG= -3
|
|---|
| 414 | IF(LP.GT.0) WRITE(LP,240)
|
|---|
| 415 | C
|
|---|
| 416 | C FORMATS
|
|---|
| 417 | C -------
|
|---|
| 418 | C
|
|---|
| 419 | 200 FORMAT(/' IFLAG= -1 ',/' LINE SEARCH FAILED. SEE'
|
|---|
| 420 | . ' DOCUMENTATION OF ROUTINE MCSRCH',/' ERROR RETURN'
|
|---|
| 421 | . ' OF LINE SEARCH: INFO= ',I2,/
|
|---|
| 422 | . ' POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT',/,
|
|---|
| 423 | . ' OR INCORRECT TOLERANCES')
|
|---|
| 424 | 235 FORMAT(/' IFLAG= -2',/' THE',I5,'-TH DIAGONAL ELEMENT OF THE',/,
|
|---|
| 425 | . ' INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE')
|
|---|
| 426 | 240 FORMAT(/' IFLAG= -3',/' IMPROPER INPUT PARAMETERS (N OR M',
|
|---|
| 427 | . ' ARE NOT POSITIVE)')
|
|---|
| 428 | 245 FORMAT(/' GTOL IS LESS THAN OR EQUAL TO 1.D-04',
|
|---|
| 429 | . / ' IT HAS BEEN RESET TO 9.D-01')
|
|---|
| 430 | RETURN
|
|---|
| 431 | END
|
|---|
| 432 | C
|
|---|
| 433 | C LAST LINE OF SUBROUTINE LBFGS
|
|---|
| 434 | C
|
|---|
| 435 | C
|
|---|
| 436 | SUBROUTINE LB1(IPRINT,ITER,NFUN,
|
|---|
| 437 | * GNORM,N,M,X,F,G,STP,FINISH)
|
|---|
| 438 | C
|
|---|
| 439 | C -------------------------------------------------------------
|
|---|
| 440 | C THIS ROUTINE PRINTS MONITORING INFORMATION. THE FREQUENCY AND
|
|---|
| 441 | C AMOUNT OF OUTPUT ARE CONTROLLED BY IPRINT.
|
|---|
| 442 | C -------------------------------------------------------------
|
|---|
| 443 | C
|
|---|
| 444 | INTEGER IPRINT(2),ITER,NFUN,LP,MP,N,M
|
|---|
| 445 | DOUBLE PRECISION X(N),G(N),F,GNORM,STP,GTOL,STPMIN,STPMAX
|
|---|
| 446 | LOGICAL FINISH
|
|---|
| 447 | COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
|
|---|
| 448 | C
|
|---|
| 449 | IF (ITER.EQ.0)THEN
|
|---|
| 450 | WRITE(MP,10)
|
|---|
| 451 | WRITE(MP,20) N,M
|
|---|
| 452 | WRITE(MP,30)F,GNORM
|
|---|
| 453 | IF (IPRINT(2).GE.1)THEN
|
|---|
| 454 | WRITE(MP,40)
|
|---|
| 455 | WRITE(MP,50) (X(I),I=1,N)
|
|---|
| 456 | WRITE(MP,60)
|
|---|
| 457 | WRITE(MP,50) (G(I),I=1,N)
|
|---|
| 458 | ENDIF
|
|---|
| 459 | WRITE(MP,10)
|
|---|
| 460 | WRITE(MP,70)
|
|---|
| 461 | ELSE
|
|---|
| 462 | IF ((IPRINT(1).EQ.0).AND.(ITER.NE.1.AND..NOT.FINISH))RETURN
|
|---|
| 463 | IF (IPRINT(1).NE.0)THEN
|
|---|
| 464 | IF(MOD(ITER-1,IPRINT(1)).EQ.0.OR.FINISH)THEN
|
|---|
| 465 | IF(IPRINT(2).GT.1.AND.ITER.GT.1) WRITE(MP,70)
|
|---|
| 466 | WRITE(MP,80)ITER,NFUN,F,GNORM,STP
|
|---|
| 467 | ELSE
|
|---|
| 468 | RETURN
|
|---|
| 469 | ENDIF
|
|---|
| 470 | ELSE
|
|---|
| 471 | IF( IPRINT(2).GT.1.AND.FINISH) WRITE(MP,70)
|
|---|
| 472 | WRITE(MP,80)ITER,NFUN,F,GNORM,STP
|
|---|
| 473 | ENDIF
|
|---|
| 474 | IF (IPRINT(2).EQ.2.OR.IPRINT(2).EQ.3)THEN
|
|---|
| 475 | IF (FINISH)THEN
|
|---|
| 476 | WRITE(MP,90)
|
|---|
| 477 | ELSE
|
|---|
| 478 | WRITE(MP,40)
|
|---|
| 479 | ENDIF
|
|---|
| 480 | WRITE(MP,50)(X(I),I=1,N)
|
|---|
| 481 | IF (IPRINT(2).EQ.3)THEN
|
|---|
| 482 | WRITE(MP,60)
|
|---|
| 483 | WRITE(MP,50)(G(I),I=1,N)
|
|---|
| 484 | ENDIF
|
|---|
| 485 | ENDIF
|
|---|
| 486 | IF (FINISH) WRITE(MP,100)
|
|---|
| 487 | ENDIF
|
|---|
| 488 | C
|
|---|
| 489 | 10 FORMAT('*************************************************')
|
|---|
| 490 | 20 FORMAT(' N=',I5,' NUMBER OF CORRECTIONS=',I2,
|
|---|
| 491 | . /, ' INITIAL VALUES')
|
|---|
| 492 | 30 FORMAT(' F= ',1PD10.3,' GNORM= ',1PD10.3)
|
|---|
| 493 | 40 FORMAT(' VECTOR X= ')
|
|---|
| 494 | 50 FORMAT(6(2X,1PD10.3))
|
|---|
| 495 | 60 FORMAT(' GRADIENT VECTOR G= ')
|
|---|
| 496 | 70 FORMAT(/' I NFN',4X,'FUNC',8X,'GNORM',7X,'STEPLENGTH'/)
|
|---|
| 497 | 80 FORMAT(2(I4,1X),3X,3(1PD10.3,2X))
|
|---|
| 498 | 90 FORMAT(' FINAL POINT X= ')
|
|---|
| 499 | 100 FORMAT(/' THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.',
|
|---|
| 500 | . /' IFLAG = 0')
|
|---|
| 501 | C
|
|---|
| 502 | RETURN
|
|---|
| 503 | END
|
|---|
| 504 | C ******
|
|---|
| 505 | C
|
|---|
| 506 | C
|
|---|
| 507 | C ----------------------------------------------------------
|
|---|
| 508 | C DATA
|
|---|
| 509 | C ----------------------------------------------------------
|
|---|
| 510 | C
|
|---|
| 511 | BLOCK DATA LB2
|
|---|
| 512 | INTEGER LP,MP
|
|---|
| 513 | DOUBLE PRECISION GTOL,STPMIN,STPMAX
|
|---|
| 514 | COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
|
|---|
| 515 | DATA MP,LP,GTOL,STPMIN,STPMAX/6,6,9.0D-01,1.0D-20,1.0D+20/
|
|---|
| 516 | END
|
|---|
| 517 | C
|
|---|
| 518 | C
|
|---|
| 519 | C ----------------------------------------------------------
|
|---|
| 520 | C
|
|---|
| 521 | subroutine daxpy(n,da,dx,incx,dy,incy)
|
|---|
| 522 | c
|
|---|
| 523 | c constant times a vector plus a vector.
|
|---|
| 524 | c uses unrolled loops for increments equal to one.
|
|---|
| 525 | c jack dongarra, linpack, 3/11/78.
|
|---|
| 526 | c
|
|---|
| 527 | double precision dx(1),dy(1),da
|
|---|
| 528 | integer i,incx,incy,ix,iy,m,mp1,n
|
|---|
| 529 | c
|
|---|
| 530 | if(n.le.0)return
|
|---|
| 531 | if (da .eq. 0.0d0) return
|
|---|
| 532 | if(incx.eq.1.and.incy.eq.1)go to 20
|
|---|
| 533 | c
|
|---|
| 534 | c code for unequal increments or equal increments
|
|---|
| 535 | c not equal to 1
|
|---|
| 536 | c
|
|---|
| 537 | ix = 1
|
|---|
| 538 | iy = 1
|
|---|
| 539 | if(incx.lt.0)ix = (-n+1)*incx + 1
|
|---|
| 540 | if(incy.lt.0)iy = (-n+1)*incy + 1
|
|---|
| 541 | do 10 i = 1,n
|
|---|
| 542 | dy(iy) = dy(iy) + da*dx(ix)
|
|---|
| 543 | ix = ix + incx
|
|---|
| 544 | iy = iy + incy
|
|---|
| 545 | 10 continue
|
|---|
| 546 | return
|
|---|
| 547 | c
|
|---|
| 548 | c code for both increments equal to 1
|
|---|
| 549 | c
|
|---|
| 550 | c
|
|---|
| 551 | c clean-up loop
|
|---|
| 552 | c
|
|---|
| 553 | 20 m = mod(n,4)
|
|---|
| 554 | if( m .eq. 0 ) go to 40
|
|---|
| 555 | do 30 i = 1,m
|
|---|
| 556 | dy(i) = dy(i) + da*dx(i)
|
|---|
| 557 | 30 continue
|
|---|
| 558 | if( n .lt. 4 ) return
|
|---|
| 559 | 40 mp1 = m + 1
|
|---|
| 560 | do 50 i = mp1,n,4
|
|---|
| 561 | dy(i) = dy(i) + da*dx(i)
|
|---|
| 562 | dy(i + 1) = dy(i + 1) + da*dx(i + 1)
|
|---|
| 563 | dy(i + 2) = dy(i + 2) + da*dx(i + 2)
|
|---|
| 564 | dy(i + 3) = dy(i + 3) + da*dx(i + 3)
|
|---|
| 565 | 50 continue
|
|---|
| 566 | return
|
|---|
| 567 | end
|
|---|
| 568 | C
|
|---|
| 569 | C
|
|---|
| 570 | C ----------------------------------------------------------
|
|---|
| 571 | C
|
|---|
| 572 | double precision function ddot(n,dx,incx,dy,incy)
|
|---|
| 573 | c
|
|---|
| 574 | c forms the dot product of two vectors.
|
|---|
| 575 | c uses unrolled loops for increments equal to one.
|
|---|
| 576 | c jack dongarra, linpack, 3/11/78.
|
|---|
| 577 | c
|
|---|
| 578 | double precision dx(1),dy(1),dtemp
|
|---|
| 579 | integer i,incx,incy,ix,iy,m,mp1,n
|
|---|
| 580 | c
|
|---|
| 581 | ddot = 0.0d0
|
|---|
| 582 | dtemp = 0.0d0
|
|---|
| 583 | if(n.le.0)return
|
|---|
| 584 | if(incx.eq.1.and.incy.eq.1)go to 20
|
|---|
| 585 | c
|
|---|
| 586 | c code for unequal increments or equal increments
|
|---|
| 587 | c not equal to 1
|
|---|
| 588 | c
|
|---|
| 589 | ix = 1
|
|---|
| 590 | iy = 1
|
|---|
| 591 | if(incx.lt.0)ix = (-n+1)*incx + 1
|
|---|
| 592 | if(incy.lt.0)iy = (-n+1)*incy + 1
|
|---|
| 593 | do 10 i = 1,n
|
|---|
| 594 | dtemp = dtemp + dx(ix)*dy(iy)
|
|---|
| 595 | ix = ix + incx
|
|---|
| 596 | iy = iy + incy
|
|---|
| 597 | 10 continue
|
|---|
| 598 | ddot = dtemp
|
|---|
| 599 | return
|
|---|
| 600 | c
|
|---|
| 601 | c code for both increments equal to 1
|
|---|
| 602 | c
|
|---|
| 603 | c
|
|---|
| 604 | c clean-up loop
|
|---|
| 605 | c
|
|---|
| 606 | 20 m = mod(n,5)
|
|---|
| 607 | if( m .eq. 0 ) go to 40
|
|---|
| 608 | do 30 i = 1,m
|
|---|
| 609 | dtemp = dtemp + dx(i)*dy(i)
|
|---|
| 610 | 30 continue
|
|---|
| 611 | if( n .lt. 5 ) go to 60
|
|---|
| 612 | 40 mp1 = m + 1
|
|---|
| 613 | do 50 i = mp1,n,5
|
|---|
| 614 | dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
|
|---|
| 615 | * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
|
|---|
| 616 | 50 continue
|
|---|
| 617 | 60 ddot = dtemp
|
|---|
| 618 | return
|
|---|
| 619 | end
|
|---|
| 620 | C ------------------------------------------------------------------
|
|---|
| 621 | C
|
|---|
| 622 | C **************************
|
|---|
| 623 | C LINE SEARCH ROUTINE MCSRCH
|
|---|
| 624 | C **************************
|
|---|
| 625 | C
|
|---|
| 626 | SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL,MAXFEV,INFO,NFEV,WA)
|
|---|
| 627 | INTEGER N,MAXFEV,INFO,NFEV
|
|---|
| 628 | DOUBLE PRECISION F,STP,FTOL,GTOL,XTOL,STPMIN,STPMAX
|
|---|
| 629 | DOUBLE PRECISION X(N),G(N),S(N),WA(N)
|
|---|
| 630 | COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
|
|---|
| 631 | SAVE
|
|---|
| 632 | C
|
|---|
| 633 | C SUBROUTINE MCSRCH
|
|---|
| 634 | C
|
|---|
| 635 | C A slight modification of the subroutine CSRCH of More' and Thuente.
|
|---|
| 636 | C The changes are to allow reverse communication, and do not affect
|
|---|
| 637 | C the performance of the routine.
|
|---|
| 638 | C
|
|---|
| 639 | C THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES
|
|---|
| 640 | C A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION.
|
|---|
| 641 | C
|
|---|
| 642 | C AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF
|
|---|
| 643 | C UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF
|
|---|
| 644 | C UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A
|
|---|
| 645 | C MINIMIZER OF THE MODIFIED FUNCTION
|
|---|
| 646 | C
|
|---|
| 647 | C F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
|
|---|
| 648 | C
|
|---|
| 649 | C IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION
|
|---|
| 650 | C HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE,
|
|---|
| 651 | C THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT
|
|---|
| 652 | C CONTAINS A MINIMIZER OF F(X+STP*S).
|
|---|
| 653 | C
|
|---|
| 654 | C THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES
|
|---|
| 655 | C THE SUFFICIENT DECREASE CONDITION
|
|---|
| 656 | C
|
|---|
| 657 | C F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S),
|
|---|
| 658 | C
|
|---|
| 659 | C AND THE CURVATURE CONDITION
|
|---|
| 660 | C
|
|---|
| 661 | C ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S).
|
|---|
| 662 | C
|
|---|
| 663 | C IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION
|
|---|
| 664 | C IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES
|
|---|
| 665 | C BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH
|
|---|
| 666 | C CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING
|
|---|
| 667 | C ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY
|
|---|
| 668 | C SATISFIES THE SUFFICIENT DECREASE CONDITION.
|
|---|
| 669 | C
|
|---|
| 670 | C THE SUBROUTINE STATEMENT IS
|
|---|
| 671 | C
|
|---|
| 672 | C SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL, MAXFEV,INFO,NFEV,WA)
|
|---|
| 673 | C WHERE
|
|---|
| 674 | C
|
|---|
| 675 | C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
|
|---|
| 676 | C OF VARIABLES.
|
|---|
| 677 | C
|
|---|
| 678 | C X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
|
|---|
| 679 | C BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS
|
|---|
| 680 | C X + STP*S.
|
|---|
| 681 | C
|
|---|
| 682 | C F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F
|
|---|
| 683 | C AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S.
|
|---|
| 684 | C
|
|---|
| 685 | C G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
|
|---|
| 686 | C GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT
|
|---|
| 687 | C OF F AT X + STP*S.
|
|---|
| 688 | C
|
|---|
| 689 | C S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE
|
|---|
| 690 | C SEARCH DIRECTION.
|
|---|
| 691 | C
|
|---|
| 692 | C STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN
|
|---|
| 693 | C INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT
|
|---|
| 694 | C STP CONTAINS THE FINAL ESTIMATE.
|
|---|
| 695 | C
|
|---|
| 696 | C FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. (In this reverse
|
|---|
| 697 | C communication implementation GTOL is defined in a COMMON
|
|---|
| 698 | C statement.) TERMINATION OCCURS WHEN THE SUFFICIENT DECREASE
|
|---|
| 699 | C CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
|
|---|
| 700 | C SATISFIED.
|
|---|
| 701 | C
|
|---|
| 702 | C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
|
|---|
| 703 | C WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
|
|---|
| 704 | C IS AT MOST XTOL.
|
|---|
| 705 | C
|
|---|
| 706 | C STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH
|
|---|
| 707 | C SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. (In this reverse
|
|---|
| 708 | C communication implementatin they are defined in a COMMON
|
|---|
| 709 | C statement).
|
|---|
| 710 | C
|
|---|
| 711 | C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
|
|---|
| 712 | C OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST
|
|---|
| 713 | C MAXFEV BY THE END OF AN ITERATION.
|
|---|
| 714 | C
|
|---|
| 715 | C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
|
|---|
| 716 | C
|
|---|
| 717 | C INFO = 0 IMPROPER INPUT PARAMETERS.
|
|---|
| 718 | C
|
|---|
| 719 | C INFO =-1 A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT.
|
|---|
| 720 | C
|
|---|
| 721 | C INFO = 1 THE SUFFICIENT DECREASE CONDITION AND THE
|
|---|
| 722 | C DIRECTIONAL DERIVATIVE CONDITION HOLD.
|
|---|
| 723 | C
|
|---|
| 724 | C INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
|
|---|
| 725 | C IS AT MOST XTOL.
|
|---|
| 726 | C
|
|---|
| 727 | C INFO = 3 NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
|
|---|
| 728 | C
|
|---|
| 729 | C INFO = 4 THE STEP IS AT THE LOWER BOUND STPMIN.
|
|---|
| 730 | C
|
|---|
| 731 | C INFO = 5 THE STEP IS AT THE UPPER BOUND STPMAX.
|
|---|
| 732 | C
|
|---|
| 733 | C INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.
|
|---|
| 734 | C THERE MAY NOT BE A STEP WHICH SATISFIES THE
|
|---|
| 735 | C SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
|
|---|
| 736 | C TOLERANCES MAY BE TOO SMALL.
|
|---|
| 737 | C
|
|---|
| 738 | C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
|
|---|
| 739 | C CALLS TO FCN.
|
|---|
| 740 | C
|
|---|
| 741 | C WA IS A WORK ARRAY OF LENGTH N.
|
|---|
| 742 | C
|
|---|
| 743 | C SUBPROGRAMS CALLED
|
|---|
| 744 | C
|
|---|
| 745 | C MCSTEP
|
|---|
| 746 | C
|
|---|
| 747 | C FORTRAN-SUPPLIED...ABS,MAX,MIN
|
|---|
| 748 | C
|
|---|
| 749 | C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
|
|---|
| 750 | C JORGE J. MORE', DAVID J. THUENTE
|
|---|
| 751 | C
|
|---|
| 752 | C **********
|
|---|
| 753 | INTEGER INFOC,J
|
|---|
| 754 | LOGICAL BRACKT,STAGE1
|
|---|
| 755 | DOUBLE PRECISION DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,
|
|---|
| 756 | * FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY,
|
|---|
| 757 | * STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,ZERO
|
|---|
| 758 | DATA P5,P66,XTRAPF,ZERO /0.5D0,0.66D0,4.0D0,0.0D0/
|
|---|
| 759 | IF(INFO.EQ.-1) GO TO 45
|
|---|
| 760 | INFOC = 1
|
|---|
| 761 | C
|
|---|
| 762 | C CHECK THE INPUT PARAMETERS FOR ERRORS.
|
|---|
| 763 | C
|
|---|
| 764 | IF (N .LE. 0 .OR. STP .LE. ZERO .OR. FTOL .LT. ZERO .OR.
|
|---|
| 765 | * GTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. STPMIN .LT. ZERO
|
|---|
| 766 | * .OR. STPMAX .LT. STPMIN .OR. MAXFEV .LE. 0) RETURN
|
|---|
| 767 | C
|
|---|
| 768 | C COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
|
|---|
| 769 | C AND CHECK THAT S IS A DESCENT DIRECTION.
|
|---|
| 770 | C
|
|---|
| 771 | DGINIT = ZERO
|
|---|
| 772 | DO 10 J = 1, N
|
|---|
| 773 | DGINIT = DGINIT + G(J)*S(J)
|
|---|
| 774 | 10 CONTINUE
|
|---|
| 775 | IF (DGINIT .GE. ZERO) then
|
|---|
| 776 | write(LP,15)
|
|---|
| 777 | 15 FORMAT(/' THE SEARCH DIRECTION IS NOT A DESCENT DIRECTION')
|
|---|
| 778 | RETURN
|
|---|
| 779 | ENDIF
|
|---|
| 780 | C
|
|---|
| 781 | C INITIALIZE LOCAL VARIABLES.
|
|---|
| 782 | C
|
|---|
| 783 | BRACKT = .FALSE.
|
|---|
| 784 | STAGE1 = .TRUE.
|
|---|
| 785 | NFEV = 0
|
|---|
| 786 | FINIT = F
|
|---|
| 787 | DGTEST = FTOL*DGINIT
|
|---|
| 788 | WIDTH = STPMAX - STPMIN
|
|---|
| 789 | WIDTH1 = WIDTH/P5
|
|---|
| 790 | DO 20 J = 1, N
|
|---|
| 791 | WA(J) = X(J)
|
|---|
| 792 | 20 CONTINUE
|
|---|
| 793 | C
|
|---|
| 794 | C THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
|
|---|
| 795 | C FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
|
|---|
| 796 | C THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
|
|---|
| 797 | C FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
|
|---|
| 798 | C THE INTERVAL OF UNCERTAINTY.
|
|---|
| 799 | C THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
|
|---|
| 800 | C FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
|
|---|
| 801 | C
|
|---|
| 802 | STX = ZERO
|
|---|
| 803 | FX = FINIT
|
|---|
| 804 | DGX = DGINIT
|
|---|
| 805 | STY = ZERO
|
|---|
| 806 | FY = FINIT
|
|---|
| 807 | DGY = DGINIT
|
|---|
| 808 | C
|
|---|
| 809 | C START OF ITERATION.
|
|---|
| 810 | C
|
|---|
| 811 | 30 CONTINUE
|
|---|
| 812 | C
|
|---|
| 813 | C SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
|
|---|
| 814 | C TO THE PRESENT INTERVAL OF UNCERTAINTY.
|
|---|
| 815 | C
|
|---|
| 816 | IF (BRACKT) THEN
|
|---|
| 817 | STMIN = MIN(STX,STY)
|
|---|
| 818 | STMAX = MAX(STX,STY)
|
|---|
| 819 | ELSE
|
|---|
| 820 | STMIN = STX
|
|---|
| 821 | STMAX = STP + XTRAPF*(STP - STX)
|
|---|
| 822 | END IF
|
|---|
| 823 | C
|
|---|
| 824 | C FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
|
|---|
| 825 | C
|
|---|
| 826 | STP = MAX(STP,STPMIN)
|
|---|
| 827 | STP = MIN(STP,STPMAX)
|
|---|
| 828 | C
|
|---|
| 829 | C IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
|
|---|
| 830 | C STP BE THE LOWEST POINT OBTAINED SO FAR.
|
|---|
| 831 | C
|
|---|
| 832 | IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
|
|---|
| 833 | * .OR. NFEV .GE. MAXFEV-1 .OR. INFOC .EQ. 0
|
|---|
| 834 | * .OR. (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX)) STP = STX
|
|---|
| 835 | C
|
|---|
| 836 | C EVALUATE THE FUNCTION AND GRADIENT AT STP
|
|---|
| 837 | C AND COMPUTE THE DIRECTIONAL DERIVATIVE.
|
|---|
| 838 | C We return to main program to obtain F and G.
|
|---|
| 839 | C
|
|---|
| 840 | DO 40 J = 1, N
|
|---|
| 841 | X(J) = WA(J) + STP*S(J)
|
|---|
| 842 | 40 CONTINUE
|
|---|
| 843 | INFO=-1
|
|---|
| 844 | RETURN
|
|---|
| 845 | C
|
|---|
| 846 | 45 INFO=0
|
|---|
| 847 | NFEV = NFEV + 1
|
|---|
| 848 | DG = ZERO
|
|---|
| 849 | DO 50 J = 1, N
|
|---|
| 850 | DG = DG + G(J)*S(J)
|
|---|
| 851 | 50 CONTINUE
|
|---|
| 852 | FTEST1 = FINIT + STP*DGTEST
|
|---|
| 853 | C
|
|---|
| 854 | C TEST FOR CONVERGENCE.
|
|---|
| 855 | C
|
|---|
| 856 | IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
|
|---|
| 857 | * .OR. INFOC .EQ. 0) INFO = 6
|
|---|
| 858 | IF (STP .EQ. STPMAX .AND.
|
|---|
| 859 | * F .LE. FTEST1 .AND. DG .LE. DGTEST) INFO = 5
|
|---|
| 860 | IF (STP .EQ. STPMIN .AND.
|
|---|
| 861 | * (F .GT. FTEST1 .OR. DG .GE. DGTEST)) INFO = 4
|
|---|
| 862 | IF (NFEV .GE. MAXFEV) INFO = 3
|
|---|
| 863 | IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2
|
|---|
| 864 | IF (F .LE. FTEST1 .AND. ABS(DG) .LE. GTOL*(-DGINIT)) INFO = 1
|
|---|
| 865 | C
|
|---|
| 866 | C CHECK FOR TERMINATION.
|
|---|
| 867 | C
|
|---|
| 868 | IF (INFO .NE. 0) RETURN
|
|---|
| 869 | C
|
|---|
| 870 | C IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
|
|---|
| 871 | C FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
|
|---|
| 872 | C
|
|---|
| 873 | IF (STAGE1 .AND. F .LE. FTEST1 .AND.
|
|---|
| 874 | * DG .GE. MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE.
|
|---|
| 875 | C
|
|---|
| 876 | C A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
|
|---|
| 877 | C WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
|
|---|
| 878 | C FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
|
|---|
| 879 | C DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
|
|---|
| 880 | C OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
|
|---|
| 881 | C
|
|---|
| 882 | IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN
|
|---|
| 883 | C
|
|---|
| 884 | C DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
|
|---|
| 885 | C
|
|---|
| 886 | FM = F - STP*DGTEST
|
|---|
| 887 | FXM = FX - STX*DGTEST
|
|---|
| 888 | FYM = FY - STY*DGTEST
|
|---|
| 889 | DGM = DG - DGTEST
|
|---|
| 890 | DGXM = DGX - DGTEST
|
|---|
| 891 | DGYM = DGY - DGTEST
|
|---|
| 892 | C
|
|---|
| 893 | C CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
|
|---|
| 894 | C AND TO COMPUTE THE NEW STEP.
|
|---|
| 895 | C
|
|---|
| 896 | CALL MCSTEP(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,
|
|---|
| 897 | * BRACKT,STMIN,STMAX,INFOC)
|
|---|
| 898 | C
|
|---|
| 899 | C RESET THE FUNCTION AND GRADIENT VALUES FOR F.
|
|---|
| 900 | C
|
|---|
| 901 | FX = FXM + STX*DGTEST
|
|---|
| 902 | FY = FYM + STY*DGTEST
|
|---|
| 903 | DGX = DGXM + DGTEST
|
|---|
| 904 | DGY = DGYM + DGTEST
|
|---|
| 905 | ELSE
|
|---|
| 906 | C
|
|---|
| 907 | C CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
|
|---|
| 908 | C AND TO COMPUTE THE NEW STEP.
|
|---|
| 909 | C
|
|---|
| 910 | CALL MCSTEP(STX,FX,DGX,STY,FY,DGY,STP,F,DG,
|
|---|
| 911 | * BRACKT,STMIN,STMAX,INFOC)
|
|---|
| 912 | END IF
|
|---|
| 913 | C
|
|---|
| 914 | C FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
|
|---|
| 915 | C INTERVAL OF UNCERTAINTY.
|
|---|
| 916 | C
|
|---|
| 917 | IF (BRACKT) THEN
|
|---|
| 918 | IF (ABS(STY-STX) .GE. P66*WIDTH1)
|
|---|
| 919 | * STP = STX + P5*(STY - STX)
|
|---|
| 920 | WIDTH1 = WIDTH
|
|---|
| 921 | WIDTH = ABS(STY-STX)
|
|---|
| 922 | END IF
|
|---|
| 923 | C
|
|---|
| 924 | C END OF ITERATION.
|
|---|
| 925 | C
|
|---|
| 926 | GO TO 30
|
|---|
| 927 | C
|
|---|
| 928 | C LAST LINE OF SUBROUTINE MCSRCH.
|
|---|
| 929 | C
|
|---|
| 930 | END
|
|---|
| 931 | SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
|
|---|
| 932 | * STPMIN,STPMAX,INFO)
|
|---|
| 933 | INTEGER INFO
|
|---|
| 934 | DOUBLE PRECISION STX,FX,DX,STY,FY,DY,STP,FP,DP,STPMIN,STPMAX
|
|---|
| 935 | LOGICAL BRACKT,BOUND
|
|---|
| 936 | C
|
|---|
| 937 | C SUBROUTINE MCSTEP
|
|---|
| 938 | C
|
|---|
| 939 | C THE PURPOSE OF MCSTEP IS TO COMPUTE A SAFEGUARDED STEP FOR
|
|---|
| 940 | C A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR
|
|---|
| 941 | C A MINIMIZER OF THE FUNCTION.
|
|---|
| 942 | C
|
|---|
| 943 | C THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION
|
|---|
| 944 | C VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS
|
|---|
| 945 | C ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE
|
|---|
| 946 | C DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A
|
|---|
| 947 | C MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY
|
|---|
| 948 | C WITH ENDPOINTS STX AND STY.
|
|---|
| 949 | C
|
|---|
| 950 | C THE SUBROUTINE STATEMENT IS
|
|---|
| 951 | C
|
|---|
| 952 | C SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
|
|---|
| 953 | C STPMIN,STPMAX,INFO)
|
|---|
| 954 | C
|
|---|
| 955 | C WHERE
|
|---|
| 956 | C
|
|---|
| 957 | C STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP,
|
|---|
| 958 | C THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED
|
|---|
| 959 | C SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION
|
|---|
| 960 | C OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE
|
|---|
| 961 | C SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY.
|
|---|
| 962 | C
|
|---|
| 963 | C STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP,
|
|---|
| 964 | C THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF
|
|---|
| 965 | C THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE
|
|---|
| 966 | C UPDATED APPROPRIATELY.
|
|---|
| 967 | C
|
|---|
| 968 | C STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP,
|
|---|
| 969 | C THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP.
|
|---|
| 970 | C IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE
|
|---|
| 971 | C BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP.
|
|---|
| 972 | C
|
|---|
| 973 | C BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER
|
|---|
| 974 | C HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED
|
|---|
| 975 | C THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER
|
|---|
| 976 | C IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE.
|
|---|
| 977 | C
|
|---|
| 978 | C STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER
|
|---|
| 979 | C AND UPPER BOUNDS FOR THE STEP.
|
|---|
| 980 | C
|
|---|
| 981 | C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
|
|---|
| 982 | C IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED
|
|---|
| 983 | C ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE
|
|---|
| 984 | C INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS.
|
|---|
| 985 | C
|
|---|
| 986 | C SUBPROGRAMS CALLED
|
|---|
| 987 | C
|
|---|
| 988 | C FORTRAN-SUPPLIED ... ABS,MAX,MIN,SQRT
|
|---|
| 989 | C
|
|---|
| 990 | C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
|
|---|
| 991 | C JORGE J. MORE', DAVID J. THUENTE
|
|---|
| 992 | C
|
|---|
| 993 | DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA
|
|---|
| 994 | INFO = 0
|
|---|
| 995 | C
|
|---|
| 996 | C CHECK THE INPUT PARAMETERS FOR ERRORS.
|
|---|
| 997 | C
|
|---|
| 998 | IF ((BRACKT .AND. (STP .LE. MIN(STX,STY) .OR.
|
|---|
| 999 | * STP .GE. MAX(STX,STY))) .OR.
|
|---|
| 1000 | * DX*(STP-STX) .GE. 0.0 .OR. STPMAX .LT. STPMIN) RETURN
|
|---|
| 1001 | C
|
|---|
| 1002 | C DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
|
|---|
| 1003 | C
|
|---|
| 1004 | SGND = DP*(DX/ABS(DX))
|
|---|
| 1005 | C
|
|---|
| 1006 | C FIRST CASE. A HIGHER FUNCTION VALUE.
|
|---|
| 1007 | C THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
|
|---|
| 1008 | C TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
|
|---|
| 1009 | C ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
|
|---|
| 1010 | C
|
|---|
| 1011 | IF (FP .GT. FX) THEN
|
|---|
| 1012 | INFO = 1
|
|---|
| 1013 | BOUND = .TRUE.
|
|---|
| 1014 | THETA = 3*(FX - FP)/(STP - STX) + DX + DP
|
|---|
| 1015 | S = MAX(ABS(THETA),ABS(DX),ABS(DP))
|
|---|
| 1016 | GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S))
|
|---|
| 1017 | IF (STP .LT. STX) GAMMA = -GAMMA
|
|---|
| 1018 | P = (GAMMA - DX) + THETA
|
|---|
| 1019 | Q = ((GAMMA - DX) + GAMMA) + DP
|
|---|
| 1020 | R = P/Q
|
|---|
| 1021 | STPC = STX + R*(STP - STX)
|
|---|
| 1022 | STPQ = STX + ((DX/((FX-FP)/(STP-STX)+DX))/2)*(STP - STX)
|
|---|
| 1023 | IF (ABS(STPC-STX) .LT. ABS(STPQ-STX)) THEN
|
|---|
| 1024 | STPF = STPC
|
|---|
| 1025 | ELSE
|
|---|
| 1026 | STPF = STPC + (STPQ - STPC)/2
|
|---|
| 1027 | END IF
|
|---|
| 1028 | BRACKT = .TRUE.
|
|---|
| 1029 | C
|
|---|
| 1030 | C SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
|
|---|
| 1031 | C OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
|
|---|
| 1032 | C STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
|
|---|
| 1033 | C THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
|
|---|
| 1034 | C
|
|---|
| 1035 | ELSE IF (SGND .LT. 0.0) THEN
|
|---|
| 1036 | INFO = 2
|
|---|
| 1037 | BOUND = .FALSE.
|
|---|
| 1038 | THETA = 3*(FX - FP)/(STP - STX) + DX + DP
|
|---|
| 1039 | S = MAX(ABS(THETA),ABS(DX),ABS(DP))
|
|---|
| 1040 | GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S))
|
|---|
| 1041 | IF (STP .GT. STX) GAMMA = -GAMMA
|
|---|
| 1042 | P = (GAMMA - DP) + THETA
|
|---|
| 1043 | Q = ((GAMMA - DP) + GAMMA) + DX
|
|---|
| 1044 | R = P/Q
|
|---|
| 1045 | STPC = STP + R*(STX - STP)
|
|---|
| 1046 | STPQ = STP + (DP/(DP-DX))*(STX - STP)
|
|---|
| 1047 | IF (ABS(STPC-STP) .GT. ABS(STPQ-STP)) THEN
|
|---|
| 1048 | STPF = STPC
|
|---|
| 1049 | ELSE
|
|---|
| 1050 | STPF = STPQ
|
|---|
| 1051 | END IF
|
|---|
| 1052 | BRACKT = .TRUE.
|
|---|
| 1053 | C
|
|---|
| 1054 | C THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
|
|---|
| 1055 | C SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
|
|---|
| 1056 | C THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
|
|---|
| 1057 | C IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
|
|---|
| 1058 | C IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
|
|---|
| 1059 | C EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
|
|---|
| 1060 | C COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
|
|---|
| 1061 | C CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
|
|---|
| 1062 | C
|
|---|
| 1063 | ELSE IF (ABS(DP) .LT. ABS(DX)) THEN
|
|---|
| 1064 | INFO = 3
|
|---|
| 1065 | BOUND = .TRUE.
|
|---|
| 1066 | THETA = 3*(FX - FP)/(STP - STX) + DX + DP
|
|---|
| 1067 | S = MAX(ABS(THETA),ABS(DX),ABS(DP))
|
|---|
| 1068 | C
|
|---|
| 1069 | C THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
|
|---|
| 1070 | C TO INFINITY IN THE DIRECTION OF THE STEP.
|
|---|
| 1071 | C
|
|---|
| 1072 | GAMMA = S*SQRT(MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S)))
|
|---|
| 1073 | IF (STP .GT. STX) GAMMA = -GAMMA
|
|---|
| 1074 | P = (GAMMA - DP) + THETA
|
|---|
| 1075 | Q = (GAMMA + (DX - DP)) + GAMMA
|
|---|
| 1076 | R = P/Q
|
|---|
| 1077 | IF (R .LT. 0.0 .AND. GAMMA .NE. 0.0) THEN
|
|---|
| 1078 | STPC = STP + R*(STX - STP)
|
|---|
| 1079 | ELSE IF (STP .GT. STX) THEN
|
|---|
| 1080 | STPC = STPMAX
|
|---|
| 1081 | ELSE
|
|---|
| 1082 | STPC = STPMIN
|
|---|
| 1083 | END IF
|
|---|
| 1084 | STPQ = STP + (DP/(DP-DX))*(STX - STP)
|
|---|
| 1085 | IF (BRACKT) THEN
|
|---|
| 1086 | IF (ABS(STP-STPC) .LT. ABS(STP-STPQ)) THEN
|
|---|
| 1087 | STPF = STPC
|
|---|
| 1088 | ELSE
|
|---|
| 1089 | STPF = STPQ
|
|---|
| 1090 | END IF
|
|---|
| 1091 | ELSE
|
|---|
| 1092 | IF (ABS(STP-STPC) .GT. ABS(STP-STPQ)) THEN
|
|---|
| 1093 | STPF = STPC
|
|---|
| 1094 | ELSE
|
|---|
| 1095 | STPF = STPQ
|
|---|
| 1096 | END IF
|
|---|
| 1097 | END IF
|
|---|
| 1098 | C
|
|---|
| 1099 | C FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
|
|---|
| 1100 | C SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
|
|---|
| 1101 | C NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
|
|---|
| 1102 | C IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
|
|---|
| 1103 | C
|
|---|
| 1104 | ELSE
|
|---|
| 1105 | INFO = 4
|
|---|
| 1106 | BOUND = .FALSE.
|
|---|
| 1107 | IF (BRACKT) THEN
|
|---|
| 1108 | THETA = 3*(FP - FY)/(STY - STP) + DY + DP
|
|---|
| 1109 | S = MAX(ABS(THETA),ABS(DY),ABS(DP))
|
|---|
| 1110 | GAMMA = S*SQRT((THETA/S)**2 - (DY/S)*(DP/S))
|
|---|
| 1111 | IF (STP .GT. STY) GAMMA = -GAMMA
|
|---|
| 1112 | P = (GAMMA - DP) + THETA
|
|---|
| 1113 | Q = ((GAMMA - DP) + GAMMA) + DY
|
|---|
| 1114 | R = P/Q
|
|---|
| 1115 | STPC = STP + R*(STY - STP)
|
|---|
| 1116 | STPF = STPC
|
|---|
| 1117 | ELSE IF (STP .GT. STX) THEN
|
|---|
| 1118 | STPF = STPMAX
|
|---|
| 1119 | ELSE
|
|---|
| 1120 | STPF = STPMIN
|
|---|
| 1121 | END IF
|
|---|
| 1122 | END IF
|
|---|
| 1123 | C
|
|---|
| 1124 | C UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT
|
|---|
| 1125 | C DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE.
|
|---|
| 1126 | C
|
|---|
| 1127 | IF (FP .GT. FX) THEN
|
|---|
| 1128 | STY = STP
|
|---|
| 1129 | FY = FP
|
|---|
| 1130 | DY = DP
|
|---|
| 1131 | ELSE
|
|---|
| 1132 | IF (SGND .LT. 0.0) THEN
|
|---|
| 1133 | STY = STX
|
|---|
| 1134 | FY = FX
|
|---|
| 1135 | DY = DX
|
|---|
| 1136 | END IF
|
|---|
| 1137 | STX = STP
|
|---|
| 1138 | FX = FP
|
|---|
| 1139 | DX = DP
|
|---|
| 1140 | END IF
|
|---|
| 1141 | C
|
|---|
| 1142 | C COMPUTE THE NEW STEP AND SAFEGUARD IT.
|
|---|
| 1143 | C
|
|---|
| 1144 | STPF = MIN(STPMAX,STPF)
|
|---|
| 1145 | STPF = MAX(STPMIN,STPF)
|
|---|
| 1146 | STP = STPF
|
|---|
| 1147 | IF (BRACKT .AND. BOUND) THEN
|
|---|
| 1148 | IF (STY .GT. STX) THEN
|
|---|
| 1149 | STP = MIN(STX+0.66*(STY-STX),STP)
|
|---|
| 1150 | ELSE
|
|---|
| 1151 | STP = MAX(STX+0.66*(STY-STX),STP)
|
|---|
| 1152 | END IF
|
|---|
| 1153 | END IF
|
|---|
| 1154 | RETURN
|
|---|
| 1155 | C
|
|---|
| 1156 | C LAST LINE OF SUBROUTINE MCSTEP.
|
|---|
| 1157 | C
|
|---|
| 1158 | END
|
|---|
| 1159 |
|
|---|
