# To unbundle, sh this file echo Readme 1>&2 sed 's/.//' >Readme <<'EOF Readme' - - This directory contains the - limited memory code for solving - large unconstrained optimization - problems. - - sdrive.f - is a simple driver - lbfgs.f - contains all the supporting - routines - makefile - a makefile for Unix users - - This directory also contains the paper - liu-nocedal.ps which describes - the limited memory algorithm. - EOF Readme echo makefile 1>&2 sed 's/.//' >makefile <<'EOF makefile' -fname = sdrive.o lbfgs.o - -main : $(fname) - f77 $(fname) $(extras) -o sdrive - -.f.o :; f77 -c $*.f - EOF makefile echo sdrive.f 1>&2 sed 's/.//' >sdrive.f <<'EOF sdrive.f' -C -C *********************** -C SIMPLE DRIVER FOR LBFGS -C *********************** -C -C Example of driver for LBFGS routine, using a -C simple test problem. The solution point is at -C X=(1,...,1) and the optimal function value of 0. -C -C JORGE NOCEDAL -C *** July 1990 *** -C - PROGRAM SDRIVE - PARAMETER(NDIM=2000,MSAVE=7,NWORK=NDIM*(2*MSAVE +1)+2*MSAVE) - DOUBLE PRECISION X(NDIM),G(NDIM),DIAG(NDIM),W(NWORK) - DOUBLE PRECISION F,EPS,XTOL,GTOL,T1,T2,STPMIN,STPMAX - INTEGER IPRINT(2),IFLAG,ICALL,N,M,MP,LP,J - LOGICAL DIAGCO -C -C The driver for LBFGS must always declare LB2 as EXTERNAL -C - EXTERNAL LB2 - COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX -C - N=100 - M=5 - IPRINT(1)= 1 - IPRINT(2)= 0 -C -C We do not wish to provide the diagonal matrices Hk0, and -C therefore set DIAGCO to FALSE. -C - DIAGCO= .FALSE. - EPS= 1.0D-5 - XTOL= 1.0D-16 - ICALL=0 - IFLAG=0 - DO 10 J=1,N,2 - X(J)=-1.2D0 - X(J+1)=1.D0 - 10 CONTINUE -C - 20 CONTINUE - F= 0.D0 - DO 30 J=1,N,2 - T1= 1.D0-X(J) - T2= 1.D1*(X(J+1)-X(J)**2) - G(J+1)= 2.D1*T2 - G(J)= -2.D0*(X(J)*G(J+1)+T1) - F= F+T1**2+T2**2 - 30 CONTINUE - CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG) - IF(IFLAG.LE.0) GO TO 50 - ICALL=ICALL + 1 -C We allow at most 2000 evaluations of F and G - IF(ICALL.GT.2000) GO TO 50 - GO TO 20 - 50 CONTINUE - END -C -C ** LAST LINE OF SIMPLE DRIVER (SDRIVE) ** - EOF sdrive.f echo lbfgs.f 1>&2 sed 's/.//' >lbfgs.f <<'EOF lbfgs.f' -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 - LOGICAL 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) 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(.NOT.DIAGCO) 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) 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 - EOF lbfgs.f