C ---------------------------------------------------------------------- C This file contains the LBFGS algorithm and supporting routines C C **************** C LBFGS SUBROUTINE C **************** C SUBROUTINE LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG) C INTEGER N,M,IPRINT(2),IFLAG DOUBLE PRECISION X(N),G(N),DIAG(N),W(N*(2*M+1)+2*M) DOUBLE PRECISION F,EPS,XTOL INTEGER DIAGCO C C LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION C JORGE NOCEDAL C *** July 1990 *** C C C This subroutine solves the unconstrained minimization problem C C min F(x), x= (x1,x2,...,xN), C C using the limited memory BFGS method. The routine is especially C effective on problems involving a large number of variables. In C a typical iteration of this method an approximation Hk to the C inverse of the Hessian is obtained by applying M BFGS updates to C a diagonal matrix Hk0, using information from the previous M steps. C The user specifies the number M, which determines the amount of C storage required by the routine. The user may also provide the C diagonal matrices Hk0 if not satisfied with the default choice. C The algorithm is described in "On the limited memory BFGS method C for large scale optimization", by D. Liu and J. Nocedal, C Mathematical Programming B 45 (1989) 503-528. C C The user is required to calculate the function value F and its C gradient G. In order to allow the user complete control over C these computations, reverse communication is used. The routine C must be called repeatedly under the control of the parameter C IFLAG. C C The steplength is determined at each iteration by means of the C line search routine MCVSRCH, which is a slight modification of C the routine CSRCH written by More' and Thuente. C C The calling statement is C C CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG) C C where C C N is an INTEGER variable that must be set by the user to the C number of variables. It is not altered by the routine. C Restriction: N>0. C C M is an INTEGER variable that must be set by the user to C the number of corrections used in the BFGS update. It C is not altered by the routine. Values of M less than 3 are C not recommended; large values of M will result in excessive C computing time. 3<= M <=7 is recommended. Restriction: M>0. C C X is a DOUBLE PRECISION array of length N. On initial entry C it must be set by the user to the values of the initial C estimate of the solution vector. On exit with IFLAG=0, it C contains the values of the variables at the best point C found (usually a solution). C C F is a DOUBLE PRECISION variable. Before initial entry and on C a re-entry with IFLAG=1, it must be set by the user to C contain the value of the function F at the point X. C C G is a DOUBLE PRECISION array of length N. Before initial C entry and on a re-entry with IFLAG=1, it must be set by C the user to contain the components of the gradient G at C the point X. C C DIAGCO is a LOGICAL variable that must be set to .TRUE. if the C user wishes to provide the diagonal matrix Hk0 at each C iteration. Otherwise it should be set to .FALSE., in which C case LBFGS will use a default value described below. If C DIAGCO is set to .TRUE. the routine will return at each C iteration of the algorithm with IFLAG=2, and the diagonal C matrix Hk0 must be provided in the array DIAG. C C C DIAG is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE., C then on initial entry or on re-entry with IFLAG=2, DIAG C it must be set by the user to contain the values of the C diagonal matrix Hk0. Restriction: all elements of DIAG C must be positive. C C IPRINT is an INTEGER array of length two which must be set by the C user. C C IPRINT(1) specifies the frequency of the output: C IPRINT(1) < 0 : no output is generated, C IPRINT(1) = 0 : output only at first and last iteration, C IPRINT(1) > 0 : output every IPRINT(1) iterations. C C IPRINT(2) specifies the type of output generated: C IPRINT(2) = 0 : iteration count, number of function C evaluations, function value, norm of the C gradient, and steplength, C IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of C variables and gradient vector at the C initial point, C IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of C variables, C IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector. C C C EPS is a positive DOUBLE PRECISION variable that must be set by C the user, and determines the accuracy with which the solution C is to be found. The subroutine terminates when C C ||G|| < EPS max(1,||X||), C C where ||.|| denotes the Euclidean norm. C C XTOL is a positive DOUBLE PRECISION variable that must be set by C the user to an estimate of the machine precision (e.g. C 10**(-16) on a SUN station 3/60). The line search routine will C terminate if the relative width of the interval of uncertainty C is less than XTOL. C C W is a DOUBLE PRECISION array of length N(2M+1)+2M used as C workspace for LBFGS. This array must not be altered by the C user. C C IFLAG is an INTEGER variable that must be set to 0 on initial entry C to the subroutine. A return with IFLAG<0 indicates an error, C and IFLAG=0 indicates that the routine has terminated without C detecting errors. On a return with IFLAG=1, the user must C evaluate the function F and gradient G. On a return with C IFLAG=2, the user must provide the diagonal matrix Hk0. C C The following negative values of IFLAG, detecting an error, C are possible: C C IFLAG=-1 The line search routine MCSRCH failed. The C parameter INFO provides more detailed information C (see also the documentation of MCSRCH): C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF C UNCERTAINTY IS AT MOST XTOL. C C INFO = 3 MORE THAN 20 FUNCTION EVALUATIONS WERE C REQUIRED AT THE PRESENT ITERATION. C C INFO = 4 THE STEP IS TOO SMALL. C C INFO = 5 THE STEP IS TOO LARGE. C C INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS. C THERE MAY NOT BE A STEP WHICH SATISFIES C THE SUFFICIENT DECREASE AND CURVATURE C CONDITIONS. TOLERANCES MAY BE TOO SMALL. C C C IFLAG=-2 The i-th diagonal element of the diagonal inverse C Hessian approximation, given in DIAG, is not C positive. C C IFLAG=-3 Improper input parameters for LBFGS (N or M are C not positive). C C C C ON THE DRIVER: C C The program that calls LBFGS must contain the declaration: C C EXTERNAL LB2 C C LB2 is a BLOCK DATA that defines the default values of several C parameters described in the COMMON section. C C C C COMMON: C C The subroutine contains one common area, which the user may wish to C reference: C COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX C C MP is an INTEGER variable with default value 6. It is used as the C unit number for the printing of the monitoring information C controlled by IPRINT. C C LP is an INTEGER variable with default value 6. It is used as the C unit number for the printing of error messages. This printing C may be suppressed by setting LP to a non-positive value. C C GTOL is a DOUBLE PRECISION variable with default value 0.9, which C controls the accuracy of the line search routine MCSRCH. If the C function and gradient evaluations are inexpensive with respect C to the cost of the iteration (which is sometimes the case when C solving very large problems) it may be advantageous to set GTOL C to a small value. A typical small value is 0.1. Restriction: C GTOL should be greater than 1.D-04. C C STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which C specify lower and uper bounds for the step in the line search. C Their default values are 1.D-20 and 1.D+20, respectively. These C values need not be modified unless the exponents are too large C for the machine being used, or unless the problem is extremely C badly scaled (in which case the exponents should be increased). C C C MACHINE DEPENDENCIES C C The only variables that are machine-dependent are XTOL, C STPMIN and STPMAX. C C C GENERAL INFORMATION C C Other routines called directly: DAXPY, DDOT, LB1, MCSRCH C C Input/Output : No input; diagnostic messages on unit MP and C error messages on unit LP. C C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C DOUBLE PRECISION GTOL,ONE,ZERO,GNORM,DDOT,STP1,FTOL,STPMIN, . STPMAX,STP,YS,YY,SQ,YR,BETA,XNORM INTEGER MP,LP,ITER,NFUN,POINT,ISPT,IYPT,MAXFEV,INFO, . BOUND,NPT,CP,I,NFEV,INMC,IYCN,ISCN LOGICAL FINISH C SAVE DATA ONE,ZERO/1.0D+0,0.0D+0/ C C INITIALIZE C ---------- C IF(IFLAG.EQ.0) GO TO 10 GO TO (172,100) IFLAG 10 ITER= 0 IF(N.LE.0.OR.M.LE.0) GO TO 196 IF(GTOL.LE.1.D-04) THEN IF(LP.GT.0) WRITE(LP,245) GTOL=9.D-01 ENDIF NFUN= 1 POINT= 0 FINISH= .FALSE. IF(DIAGCO.NE.0) THEN DO 30 I=1,N 30 IF (DIAG(I).LE.ZERO) GO TO 195 ELSE DO 40 I=1,N 40 DIAG(I)= 1.0D0 ENDIF C C THE WORK VECTOR W IS DIVIDED AS FOLLOWS: C --------------------------------------- C THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT AND C OTHER TEMPORARY INFORMATION. C LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO. C LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USED C IN THE FORMULA THAT COMPUTES H*G. C LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCH C STEPS. C LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST M C GRADIENT DIFFERENCES. C C THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN A C CIRCULAR ORDER CONTROLLED BY THE PARAMETER POINT. C ISPT= N+2*M IYPT= ISPT+N*M DO 50 I=1,N 50 W(ISPT+I)= -G(I)*DIAG(I) GNORM= DSQRT(DDOT(N,G,1,G,1)) STP1= ONE/GNORM C C PARAMETERS FOR LINE SEARCH ROUTINE C FTOL= 1.0D-4 MAXFEV= 20 C IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN, * GNORM,N,M,X,F,G,STP,FINISH) C C -------------------- C MAIN ITERATION LOOP C -------------------- C 80 ITER= ITER+1 INFO=0 BOUND=ITER-1 IF(ITER.EQ.1) GO TO 165 IF (ITER .GT. M)BOUND=M C YS= DDOT(N,W(IYPT+NPT+1),1,W(ISPT+NPT+1),1) IF(DIAGCO.EQ.0) THEN YY= DDOT(N,W(IYPT+NPT+1),1,W(IYPT+NPT+1),1) DO 90 I=1,N 90 DIAG(I)= YS/YY ELSE IFLAG=2 RETURN ENDIF 100 CONTINUE IF(DIAGCO.NE.0) THEN DO 110 I=1,N 110 IF (DIAG(I).LE.ZERO) GO TO 195 ENDIF C C COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, C "Updating quasi-Newton matrices with limited storage", C Mathematics of Computation, Vol.24, No.151, pp. 773-782. C --------------------------------------------------------- C CP= POINT IF (POINT.EQ.0) CP=M W(N+CP)= ONE/YS DO 112 I=1,N 112 W(I)= -G(I) CP= POINT DO 125 I= 1,BOUND CP=CP-1 IF (CP.EQ. -1)CP=M-1 SQ= DDOT(N,W(ISPT+CP*N+1),1,W,1) INMC=N+M+CP+1 IYCN=IYPT+CP*N W(INMC)= W(N+CP+1)*SQ CALL DAXPY(N,-W(INMC),W(IYCN+1),1,W,1) 125 CONTINUE C DO 130 I=1,N 130 W(I)=DIAG(I)*W(I) C DO 145 I=1,BOUND YR= DDOT(N,W(IYPT+CP*N+1),1,W,1) BETA= W(N+CP+1)*YR INMC=N+M+CP+1 BETA= W(INMC)-BETA ISCN=ISPT+CP*N CALL DAXPY(N,BETA,W(ISCN+1),1,W,1) CP=CP+1 IF (CP.EQ.M)CP=0 145 CONTINUE C C STORE THE NEW SEARCH DIRECTION C ------------------------------ C DO 160 I=1,N 160 W(ISPT+POINT*N+I)= W(I) C C OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION C BY USING THE LINE SEARCH ROUTINE MCSRCH C ---------------------------------------------------- 165 NFEV=0 STP=ONE IF (ITER.EQ.1) STP=STP1 DO 170 I=1,N 170 W(I)=G(I) 172 CONTINUE CALL MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL, * XTOL,MAXFEV,INFO,NFEV,DIAG) IF (INFO .EQ. -1) THEN IFLAG=1 RETURN ENDIF IF (INFO .NE. 1) GO TO 190 NFUN= NFUN + NFEV C C COMPUTE THE NEW STEP AND GRADIENT CHANGE C ----------------------------------------- C NPT=POINT*N DO 175 I=1,N W(ISPT+NPT+I)= STP*W(ISPT+NPT+I) 175 W(IYPT+NPT+I)= G(I)-W(I) POINT=POINT+1 IF (POINT.EQ.M)POINT=0 C C TERMINATION TEST C ---------------- C GNORM= DSQRT(DDOT(N,G,1,G,1)) XNORM= DSQRT(DDOT(N,X,1,X,1)) XNORM= DMAX1(1.0D0,XNORM) IF (GNORM/XNORM .LE. EPS) FINISH=.TRUE. C IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN, * GNORM,N,M,X,F,G,STP,FINISH) IF (FINISH) THEN IFLAG=0 RETURN ENDIF GO TO 80 C C ------------------------------------------------------------ C END OF MAIN ITERATION LOOP. ERROR EXITS. C ------------------------------------------------------------ C 190 IFLAG=-1 IF(LP.GT.0) WRITE(LP,200) INFO RETURN 195 IFLAG=-2 IF(LP.GT.0) WRITE(LP,235) I RETURN 196 IFLAG= -3 IF(LP.GT.0) WRITE(LP,240) C C FORMATS C ------- C 200 FORMAT(/' IFLAG= -1 ',/' LINE SEARCH FAILED. SEE' . ' DOCUMENTATION OF ROUTINE MCSRCH',/' ERROR RETURN' . ' OF LINE SEARCH: INFO= ',I2,/ . ' POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT',/, . ' OR INCORRECT TOLERANCES') 235 FORMAT(/' IFLAG= -2',/' THE',I5,'-TH DIAGONAL ELEMENT OF THE',/, . ' INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE') 240 FORMAT(/' IFLAG= -3',/' IMPROPER INPUT PARAMETERS (N OR M', . ' ARE NOT POSITIVE)') 245 FORMAT(/' GTOL IS LESS THAN OR EQUAL TO 1.D-04', . / ' IT HAS BEEN RESET TO 9.D-01') RETURN END C C LAST LINE OF SUBROUTINE LBFGS C C SUBROUTINE LB1(IPRINT,ITER,NFUN, * GNORM,N,M,X,F,G,STP,FINISH) C C ------------------------------------------------------------- C THIS ROUTINE PRINTS MONITORING INFORMATION. THE FREQUENCY AND C AMOUNT OF OUTPUT ARE CONTROLLED BY IPRINT. C ------------------------------------------------------------- C INTEGER IPRINT(2),ITER,NFUN,LP,MP,N,M DOUBLE PRECISION X(N),G(N),F,GNORM,STP,GTOL,STPMIN,STPMAX LOGICAL FINISH COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX C IF (ITER.EQ.0)THEN WRITE(MP,10) WRITE(MP,20) N,M WRITE(MP,30)F,GNORM IF (IPRINT(2).GE.1)THEN WRITE(MP,40) WRITE(MP,50) (X(I),I=1,N) WRITE(MP,60) WRITE(MP,50) (G(I),I=1,N) ENDIF WRITE(MP,10) WRITE(MP,70) ELSE IF ((IPRINT(1).EQ.0).AND.(ITER.NE.1.AND..NOT.FINISH))RETURN IF (IPRINT(1).NE.0)THEN IF(MOD(ITER-1,IPRINT(1)).EQ.0.OR.FINISH)THEN IF(IPRINT(2).GT.1.AND.ITER.GT.1) WRITE(MP,70) WRITE(MP,80)ITER,NFUN,F,GNORM,STP ELSE RETURN ENDIF ELSE IF( IPRINT(2).GT.1.AND.FINISH) WRITE(MP,70) WRITE(MP,80)ITER,NFUN,F,GNORM,STP ENDIF IF (IPRINT(2).EQ.2.OR.IPRINT(2).EQ.3)THEN IF (FINISH)THEN WRITE(MP,90) ELSE WRITE(MP,40) ENDIF WRITE(MP,50)(X(I),I=1,N) IF (IPRINT(2).EQ.3)THEN WRITE(MP,60) WRITE(MP,50)(G(I),I=1,N) ENDIF ENDIF IF (FINISH) WRITE(MP,100) ENDIF C 10 FORMAT('*************************************************') 20 FORMAT(' N=',I5,' NUMBER OF CORRECTIONS=',I2, . /, ' INITIAL VALUES') 30 FORMAT(' F= ',1PD10.3,' GNORM= ',1PD10.3) 40 FORMAT(' VECTOR X= ') 50 FORMAT(6(2X,1PD10.3)) 60 FORMAT(' GRADIENT VECTOR G= ') 70 FORMAT(/' I NFN',4X,'FUNC',8X,'GNORM',7X,'STEPLENGTH'/) 80 FORMAT(2(I4,1X),3X,3(1PD10.3,2X)) 90 FORMAT(' FINAL POINT X= ') 100 FORMAT(/' THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.', . /' IFLAG = 0') C RETURN END C ****** C C C ---------------------------------------------------------- C DATA C ---------------------------------------------------------- C BLOCK DATA LB2 INTEGER LP,MP DOUBLE PRECISION GTOL,STPMIN,STPMAX COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX DATA MP,LP,GTOL,STPMIN,STPMAX/6,6,9.0D-01,1.0D-20,1.0D+20/ END C C C ---------------------------------------------------------- C subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end C C C ---------------------------------------------------------- C double precision function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end C ------------------------------------------------------------------ C C ************************** C LINE SEARCH ROUTINE MCSRCH C ************************** C SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL,MAXFEV,INFO,NFEV,WA) INTEGER N,MAXFEV,INFO,NFEV DOUBLE PRECISION F,STP,FTOL,GTOL,XTOL,STPMIN,STPMAX DOUBLE PRECISION X(N),G(N),S(N),WA(N) COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX SAVE C C SUBROUTINE MCSRCH C C A slight modification of the subroutine CSRCH of More' and Thuente. C The changes are to allow reverse communication, and do not affect C the performance of the routine. C C THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES C A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION. C C AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF C UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF C UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A C MINIMIZER OF THE MODIFIED FUNCTION C C F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S). C C IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION C HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE, C THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT C CONTAINS A MINIMIZER OF F(X+STP*S). C C THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES C THE SUFFICIENT DECREASE CONDITION C C F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S), C C AND THE CURVATURE CONDITION C C ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S). C C IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION C IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES C BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH C CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING C ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY C SATISFIES THE SUFFICIENT DECREASE CONDITION. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL, MAXFEV,INFO,NFEV,WA) C WHERE C C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER C OF VARIABLES. C C X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE C BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS C X + STP*S. C C F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F C AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S. C C G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE C GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT C OF F AT X + STP*S. C C S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE C SEARCH DIRECTION. C C STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN C INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT C STP CONTAINS THE FINAL ESTIMATE. C C FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. (In this reverse C communication implementation GTOL is defined in a COMMON C statement.) TERMINATION OCCURS WHEN THE SUFFICIENT DECREASE C CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE C SATISFIED. C C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS C WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY C IS AT MOST XTOL. C C STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH C SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. (In this reverse C communication implementatin they are defined in a COMMON C statement). C C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION C OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST C MAXFEV BY THE END OF AN ITERATION. C C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS: C C INFO = 0 IMPROPER INPUT PARAMETERS. C C INFO =-1 A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT. C C INFO = 1 THE SUFFICIENT DECREASE CONDITION AND THE C DIRECTIONAL DERIVATIVE CONDITION HOLD. C C INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY C IS AT MOST XTOL. C C INFO = 3 NUMBER OF CALLS TO FCN HAS REACHED MAXFEV. C C INFO = 4 THE STEP IS AT THE LOWER BOUND STPMIN. C C INFO = 5 THE STEP IS AT THE UPPER BOUND STPMAX. C C INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS. C THERE MAY NOT BE A STEP WHICH SATISFIES THE C SUFFICIENT DECREASE AND CURVATURE CONDITIONS. C TOLERANCES MAY BE TOO SMALL. C C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF C CALLS TO FCN. C C WA IS A WORK ARRAY OF LENGTH N. C C SUBPROGRAMS CALLED C C MCSTEP C C FORTRAN-SUPPLIED...ABS,MAX,MIN C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983 C JORGE J. MORE', DAVID J. THUENTE C C ********** INTEGER INFOC,J LOGICAL BRACKT,STAGE1 DOUBLE PRECISION DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM, * FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY, * STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,ZERO DATA P5,P66,XTRAPF,ZERO /0.5D0,0.66D0,4.0D0,0.0D0/ IF(INFO.EQ.-1) GO TO 45 INFOC = 1 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF (N .LE. 0 .OR. STP .LE. ZERO .OR. FTOL .LT. ZERO .OR. * GTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. STPMIN .LT. ZERO * .OR. STPMAX .LT. STPMIN .OR. MAXFEV .LE. 0) RETURN C C COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION C AND CHECK THAT S IS A DESCENT DIRECTION. C DGINIT = ZERO DO 10 J = 1, N DGINIT = DGINIT + G(J)*S(J) 10 CONTINUE IF (DGINIT .GE. ZERO) then write(LP,15) 15 FORMAT(/' THE SEARCH DIRECTION IS NOT A DESCENT DIRECTION') RETURN ENDIF C C INITIALIZE LOCAL VARIABLES. C BRACKT = .FALSE. STAGE1 = .TRUE. NFEV = 0 FINIT = F DGTEST = FTOL*DGINIT WIDTH = STPMAX - STPMIN WIDTH1 = WIDTH/P5 DO 20 J = 1, N WA(J) = X(J) 20 CONTINUE C C THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP, C FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP. C THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP, C FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF C THE INTERVAL OF UNCERTAINTY. C THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP, C FUNCTION, AND DERIVATIVE AT THE CURRENT STEP. C STX = ZERO FX = FINIT DGX = DGINIT STY = ZERO FY = FINIT DGY = DGINIT C C START OF ITERATION. C 30 CONTINUE C C SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND C TO THE PRESENT INTERVAL OF UNCERTAINTY. C IF (BRACKT) THEN STMIN = MIN(STX,STY) STMAX = MAX(STX,STY) ELSE STMIN = STX STMAX = STP + XTRAPF*(STP - STX) END IF C C FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN. C STP = MAX(STP,STPMIN) STP = MIN(STP,STPMAX) C C IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET C STP BE THE LOWEST POINT OBTAINED SO FAR. C IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX)) * .OR. NFEV .GE. MAXFEV-1 .OR. INFOC .EQ. 0 * .OR. (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX)) STP = STX C C EVALUATE THE FUNCTION AND GRADIENT AT STP C AND COMPUTE THE DIRECTIONAL DERIVATIVE. C We return to main program to obtain F and G. C DO 40 J = 1, N X(J) = WA(J) + STP*S(J) 40 CONTINUE INFO=-1 RETURN C 45 INFO=0 NFEV = NFEV + 1 DG = ZERO DO 50 J = 1, N DG = DG + G(J)*S(J) 50 CONTINUE FTEST1 = FINIT + STP*DGTEST C C TEST FOR CONVERGENCE. C IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX)) * .OR. INFOC .EQ. 0) INFO = 6 IF (STP .EQ. STPMAX .AND. * F .LE. FTEST1 .AND. DG .LE. DGTEST) INFO = 5 IF (STP .EQ. STPMIN .AND. * (F .GT. FTEST1 .OR. DG .GE. DGTEST)) INFO = 4 IF (NFEV .GE. MAXFEV) INFO = 3 IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2 IF (F .LE. FTEST1 .AND. ABS(DG) .LE. GTOL*(-DGINIT)) INFO = 1 C C CHECK FOR TERMINATION. C IF (INFO .NE. 0) RETURN C C IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED C FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE. C IF (STAGE1 .AND. F .LE. FTEST1 .AND. * DG .GE. MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE. C C A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF C WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED C FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE C DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN C OBTAINED BUT THE DECREASE IS NOT SUFFICIENT. C IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN C C DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES. C FM = F - STP*DGTEST FXM = FX - STX*DGTEST FYM = FY - STY*DGTEST DGM = DG - DGTEST DGXM = DGX - DGTEST DGYM = DGY - DGTEST C C CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY C AND TO COMPUTE THE NEW STEP. C CALL MCSTEP(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM, * BRACKT,STMIN,STMAX,INFOC) C C RESET THE FUNCTION AND GRADIENT VALUES FOR F. C FX = FXM + STX*DGTEST FY = FYM + STY*DGTEST DGX = DGXM + DGTEST DGY = DGYM + DGTEST ELSE C C CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY C AND TO COMPUTE THE NEW STEP. C CALL MCSTEP(STX,FX,DGX,STY,FY,DGY,STP,F,DG, * BRACKT,STMIN,STMAX,INFOC) END IF C C FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE C INTERVAL OF UNCERTAINTY. C IF (BRACKT) THEN IF (ABS(STY-STX) .GE. P66*WIDTH1) * STP = STX + P5*(STY - STX) WIDTH1 = WIDTH WIDTH = ABS(STY-STX) END IF C C END OF ITERATION. C GO TO 30 C C LAST LINE OF SUBROUTINE MCSRCH. C END SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT, * STPMIN,STPMAX,INFO) INTEGER INFO DOUBLE PRECISION STX,FX,DX,STY,FY,DY,STP,FP,DP,STPMIN,STPMAX LOGICAL BRACKT,BOUND C C SUBROUTINE MCSTEP C C THE PURPOSE OF MCSTEP IS TO COMPUTE A SAFEGUARDED STEP FOR C A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR C A MINIMIZER OF THE FUNCTION. C C THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION C VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS C ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE C DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A C MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY C WITH ENDPOINTS STX AND STY. C C THE SUBROUTINE STATEMENT IS C C SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT, C STPMIN,STPMAX,INFO) C C WHERE C C STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP, C THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED C SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION C OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE C SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY. C C STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP, C THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF C THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE C UPDATED APPROPRIATELY. C C STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP, C THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP. C IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE C BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP. C C BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER C HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED C THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER C IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE. C C STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER C AND UPPER BOUNDS FOR THE STEP. C C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS: C IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED C ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE C INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS. C C SUBPROGRAMS CALLED C C FORTRAN-SUPPLIED ... ABS,MAX,MIN,SQRT C C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983 C JORGE J. MORE', DAVID J. THUENTE C DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA INFO = 0 C C CHECK THE INPUT PARAMETERS FOR ERRORS. C IF ((BRACKT .AND. (STP .LE. MIN(STX,STY) .OR. * STP .GE. MAX(STX,STY))) .OR. * DX*(STP-STX) .GE. 0.0 .OR. STPMAX .LT. STPMIN) RETURN C C DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN. C SGND = DP*(DX/ABS(DX)) C C FIRST CASE. A HIGHER FUNCTION VALUE. C THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER C TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN, C ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN. C IF (FP .GT. FX) THEN INFO = 1 BOUND = .TRUE. THETA = 3*(FX - FP)/(STP - STX) + DX + DP S = MAX(ABS(THETA),ABS(DX),ABS(DP)) GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S)) IF (STP .LT. STX) GAMMA = -GAMMA P = (GAMMA - DX) + THETA Q = ((GAMMA - DX) + GAMMA) + DP R = P/Q STPC = STX + R*(STP - STX) STPQ = STX + ((DX/((FX-FP)/(STP-STX)+DX))/2)*(STP - STX) IF (ABS(STPC-STX) .LT. ABS(STPQ-STX)) THEN STPF = STPC ELSE STPF = STPC + (STPQ - STPC)/2 END IF BRACKT = .TRUE. C C SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF C OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC C STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP, C THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN. C ELSE IF (SGND .LT. 0.0) THEN INFO = 2 BOUND = .FALSE. THETA = 3*(FX - FP)/(STP - STX) + DX + DP S = MAX(ABS(THETA),ABS(DX),ABS(DP)) GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S)) IF (STP .GT. STX) GAMMA = -GAMMA P = (GAMMA - DP) + THETA Q = ((GAMMA - DP) + GAMMA) + DX R = P/Q STPC = STP + R*(STX - STP) STPQ = STP + (DP/(DP-DX))*(STX - STP) IF (ABS(STPC-STP) .GT. ABS(STPQ-STP)) THEN STPF = STPC ELSE STPF = STPQ END IF BRACKT = .TRUE. C C THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE C SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES. C THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY C IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC C IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE C EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO C COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP C CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN. C ELSE IF (ABS(DP) .LT. ABS(DX)) THEN INFO = 3 BOUND = .TRUE. THETA = 3*(FX - FP)/(STP - STX) + DX + DP S = MAX(ABS(THETA),ABS(DX),ABS(DP)) C C THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND C TO INFINITY IN THE DIRECTION OF THE STEP. C GAMMA = S*SQRT(MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S))) IF (STP .GT. STX) GAMMA = -GAMMA P = (GAMMA - DP) + THETA Q = (GAMMA + (DX - DP)) + GAMMA R = P/Q IF (R .LT. 0.0 .AND. GAMMA .NE. 0.0) THEN STPC = STP + R*(STX - STP) ELSE IF (STP .GT. STX) THEN STPC = STPMAX ELSE STPC = STPMIN END IF STPQ = STP + (DP/(DP-DX))*(STX - STP) IF (BRACKT) THEN IF (ABS(STP-STPC) .LT. ABS(STP-STPQ)) THEN STPF = STPC ELSE STPF = STPQ END IF ELSE IF (ABS(STP-STPC) .GT. ABS(STP-STPQ)) THEN STPF = STPC ELSE STPF = STPQ END IF END IF C C FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE C SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES C NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP C IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN. C ELSE INFO = 4 BOUND = .FALSE. IF (BRACKT) THEN THETA = 3*(FP - FY)/(STY - STP) + DY + DP S = MAX(ABS(THETA),ABS(DY),ABS(DP)) GAMMA = S*SQRT((THETA/S)**2 - (DY/S)*(DP/S)) IF (STP .GT. STY) GAMMA = -GAMMA P = (GAMMA - DP) + THETA Q = ((GAMMA - DP) + GAMMA) + DY R = P/Q STPC = STP + R*(STY - STP) STPF = STPC ELSE IF (STP .GT. STX) THEN STPF = STPMAX ELSE STPF = STPMIN END IF END IF C C UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT C DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE. C IF (FP .GT. FX) THEN STY = STP FY = FP DY = DP ELSE IF (SGND .LT. 0.0) THEN STY = STX FY = FX DY = DX END IF STX = STP FX = FP DX = DP END IF C C COMPUTE THE NEW STEP AND SAFEGUARD IT. C STPF = MIN(STPMAX,STPF) STPF = MAX(STPMIN,STPF) STP = STPF IF (BRACKT .AND. BOUND) THEN IF (STY .GT. STX) THEN STP = MIN(STX+0.66*(STY-STX),STP) ELSE STP = MAX(STX+0.66*(STY-STX),STP) END IF END IF RETURN C C LAST LINE OF SUBROUTINE MCSTEP. C END