Logo Search packages:      
Sourcecode: objcryst-fox version File versions  Download package

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
      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 methodC      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      whereC C     N       is an INTEGER variable that must be set by the user to theC             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 toC             the number of corrections used in the BFGS update. ItC             is not altered by the routine. Values of M less than 3 areC             not recommended; large values of M will result in excessiveC             computing time. 3<= M <=7 is recommended. Restriction: M>0.C C     X       is a DOUBLE PRECISION array of length N. On initial entryC             it must be set by the user to the values of the initialC             estimate of the solution vector. On exit with IFLAG=0, itC             contains the values of the variables at the best pointC             found (usually a solution).C C     F       is a DOUBLE PRECISION variable. Before initial entry and onC             a re-entry with IFLAG=1, it must be set by the user toC             contain the value of the function F at the point X.C C     G       is a DOUBLE PRECISION array of length N. Before initialC             entry and on a re-entry with IFLAG=1, it must be set byC             the user to contain the components of the gradient G atC             the point X.C C     DIAGCO  is a LOGICAL variable that must be set to .TRUE. if theC             user  wishes to provide the diagonal matrix Hk0 at eachC             iteration. Otherwise it should be set to .FALSE., in whichC             case  LBFGS will use a default value described below. IfC             DIAGCO is set to .TRUE. the routine will return at eachC             iteration of the algorithm with IFLAG=2, and the diagonalC              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, DIAGC             it must be set by the user to contain the values of the C             diagonal matrix Hk0.  Restriction: all elements of DIAGC             must be positive.C C     IPRINT  is an INTEGER array of length two which must be set by theC             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 theC                                gradient, and steplength,C                IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector ofC                                variables and  gradient vector at theC                                initial point,C                IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector ofC                                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 byC             the user, and determines the accuracy with which the solutionC             is to be found. The subroutine terminates whenCC                         ||G|| < EPS max(1,||X||),CC             where ||.|| denotes the Euclidean norm.C C     XTOL    is a  positive DOUBLE PRECISION variable that must be set byC             the user to an estimate of the machine precision (e.g.C             10**(-16) on a SUN station 3/60). The line search routine willC             terminate if the relative width of the interval of uncertaintyC             is less than XTOL.C C     W       is a DOUBLE PRECISION array of length N(2M+1)+2M used asC             workspace for LBFGS. This array must not be altered by theC             user.C C     IFLAG   is an INTEGER variable that must be set to 0 on initial entryC             to the subroutine. A return with IFLAG<0 indicates an error,C             and IFLAG=0 indicates that the routine has terminated withoutC             detecting errors. On a return with IFLAG=1, the user mustC             evaluate the function F and gradient G. On a return withC             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. TheC                        parameter INFO provides more detailed informationC                        (see also the documentation of MCSRCH):CC                       INFO = 0  IMPROPER INPUT PARAMETERS.CC                       INFO = 2  RELATIVE WIDTH OF THE INTERVAL OFC                                 UNCERTAINTY IS AT MOST XTOL.CC                       INFO = 3  MORE THAN 20 FUNCTION EVALUATIONS WEREC                                 REQUIRED AT THE PRESENT ITERATION.CC                       INFO = 4  THE STEP IS TOO SMALL.CC                       INFO = 5  THE STEP IS TOO LARGE.CC                       INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS. C                                 THERE MAY NOT BE A STEP WHICH SATISFIESC                                 THE SUFFICIENT DECREASE AND CURVATUREC                                 CONDITIONS. TOLERANCES MAY BE TOO SMALL.CC C              IFLAG=-2  The i-th diagonal element of the diagonal inverseC                        Hessian approximation, given in DIAG, is notC                        positive.C           C              IFLAG=-3  Improper input parameters for LBFGS (N or M areC                        not positive).C CCC    ON THE DRIVER:CC    The program that calls LBFGS must contain the declaration:CC                       EXTERNAL LB2CC    LB2 is a BLOCK DATA that defines the default values of severalC    parameters described in the COMMON section. CC C C    COMMON:C C     The subroutine contains one common area, which the user may wish toC    reference:C          COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAXC C    MP  is an INTEGER variable with default value 6. It is used as theC        unit number for the printing of the monitoring informationC        controlled by IPRINT.C C    LP  is an INTEGER variable with default value 6. It is used as theC        unit number for the printing of error messages. This printingC        may be suppressed by setting LP to a non-positive value.C C    GTOL is a DOUBLE PRECISION variable with default value 0.9, whichC        controls the accuracy of the line search routine MCSRCH. If theC        function and gradient evaluations are inexpensive with respectC        to the cost of the iteration (which is sometimes the case whenC        solving very large problems) it may be advantageous to set GTOLC        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 whichC        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. TheseC        values need not be modified unless the exponents are too largeC        for the machine being used, or unless the problem is extremelyC        badly scaled (in which case the exponents should be increased).C CC  MACHINE DEPENDENCIESCC        The only variables that are machine-dependent are XTOL,C        STPMIN and STPMAX.C CC  GENERAL INFORMATIONC C    Other routines called directly:  DAXPY, DDOT, LB1, MCSRCHC C    Input/Output  :  No input; diagnostic messages on unit MP andC                     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 FINISHC      SAVE      DATA ONE,ZERO/1.0D+0,0.0D+0/CC     INITIALIZEC     ----------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      ENDIFCC     THE WORK VECTOR W IS DIVIDED AS FOLLOWS:C     ---------------------------------------C     THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT ANDC         OTHER TEMPORARY INFORMATION.C     LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO.C     LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USEDC         IN THE FORMULA THAT COMPUTES H*G.C     LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCHC         STEPS.C     LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST MC         GRADIENT DIFFERENCES.CC     THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN AC     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/GNORMCC     PARAMETERS FOR LINE SEARCH ROUTINEC           FTOL= 1.0D-4      MAXFEV= 20C      IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,     *                     GNORM,N,M,X,F,G,STP,FINISH)CC    --------------------C     MAIN ITERATION LOOPC    --------------------C 80   ITER= ITER+1      INFO=0      BOUND=ITER-1      IF(ITER.EQ.1) GO TO 165      IF (ITER .GT. M)BOUND=MC         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      ENDIFCC     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  CONTINUEC      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  CONTINUECC     STORE THE NEW SEARCH DIRECTIONC     ------------------------------C       DO 160 I=1,N 160   W(ISPT+POINT*N+I)= W(I)CC     OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION C     BY USING THE LINE SEARCH ROUTINE MCSRCHC     ---------------------------------------------------- 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 + NFEVCC     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=0CC     TERMINATION TESTC     ----------------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 80CC     ------------------------------------------------------------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)CC     FORMATSC     -------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 affectC     the performance of the routine. CC     THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIESC     A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION.CC     AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OFC     UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OFC     UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS AC     MINIMIZER OF THE MODIFIED FUNCTIONCC          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),CC     AND THE CURVATURE CONDITIONCC           ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)
























































































'S).CC     IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTIONC     IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIESC     BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTHC     CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDINGC     ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLYC     SATISFIES THE SUFFICIENT DECREASE CONDITION.CC     THE SUBROUTINE STATEMENT ISCC        SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL, MAXFEV,INFO,NFEV,WA)C     WHERECC       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBERC         OF VARIABLES.CC       X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THEC         BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINSC         X + STP*S.CC       F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF FC         AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S.CC       G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THEC         GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENTC         OF F AT X + STP*S.CC       S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THEC         SEARCH DIRECTION.CC       STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS ANC         INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUTC         STP CONTAINS THE FINAL ESTIMATE.CC       FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. (In this reverseC         communication implementation GTOL is defined in a COMMONC         statement.) TERMINATION OCCURS WHEN THE SUFFICIENT DECREASEC         CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION AREC         SATISFIED.CC       XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURSC         WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTYC         IS AT MOST XTOL.CC       STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICHC         SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. (In this reverseC         communication implementatin they are defined in a COMMONC         statement).CC       MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATIONC         OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEASTC         MAXFEV BY THE END OF AN ITERATION.CC       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:CC         INFO = 0  IMPROPER INPUT PARAMETERS.CC         INFO =-1  A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT.CC         INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THEC                   DIRECTIONAL DERIVATIVE CONDITION HOLD.CC         INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTYC                   IS AT MOST XTOL.CC         INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.CC         INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN.CC         INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX.CC         INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.C                   THERE MAY NOT BE A STEP WHICH SATISFIES THEC                   SUFFICIENT DECREASE AND CURVATURE CONDITIONS.C                   TOLERANCES MAY BE TOO SMALL.CC       NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OFC         CALLS TO FCN.CC       WA IS A WORK ARRAY OF LENGTH N.CC     SUBPROGRAMS CALLEDCC       MCSTEPCC       FORTRAN-SUPPLIED...ABS,MAX,MINCC     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983C     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









































































































































































Generated by  Doxygen 1.6.0   Back to index