dcstep.f

Go to the documentation of this file.

c> \file dcstep.f
 
c> \brief This subroutine computes a safeguarded step for a search
c>        procedure and updates an interval that contains a step that
c>        satisfies a sufficient decrease and a curvature condition.
c>
c> This subroutine computes a safeguarded step for a search
c> procedure and updates an interval that contains a step that
c> satisfies a sufficient decrease and a curvature condition.
c> 
c> The parameter stx contains the step with the least function
c> value. If brackt is set to .true. then a minimizer has
c> been bracketed in an interval with endpoints stx and sty.
c> The parameter stp contains the current step. 
c> The subroutine assumes that if brackt is set to .true. then
c> 
c>       min(stx,sty) < stp < max(stx,sty),
c> 
c> and that the derivative at stx is negative in the direction 
c> of the step.
c> 
c> @param stx On entry stx is the best step obtained so far and is an
c>               endpoint of the interval that contains the minimizer.<br/>
c>            On exit stx is the updated best step.
c> 
c> @param fx On entry fx is the function at stx.<br/>
c>           On exit fx is the function at stx.
c> 
c> @param dx On entry dx is the derivative of the function at 
c>              stx. The derivative must be negative in the direction of 
c>              the step, that is, dx and stp - stx must have opposite 
c>              signs.<br/>
c>           On exit dx is the derivative of the function at stx.
c> 
c> @param sty On entry sty is the second endpoint of the interval that 
c>               contains the minimizer.<br/>
c>            On exit sty is the updated endpoint of the interval that 
c>               contains the minimizer.
c> 
c> @param fy On entry fy is the function at sty.<br/>
c>           On exit fy is the function at sty.
c> 
c> @param dy On entry dy is the derivative of the function at sty.<br/>
c>           On exit dy is the derivative of the function at the exit sty.
c> 
c> @param stp On entry stp is the current step. If brackt is set to .true.
c>               then on input stp must be between stx and sty.<br/>
c>            On exit stp is a new trial step.
c> 
c> @param fp On entry fp is the function at stp.<br/>
c>           On exit fp is unchanged.
c> 
c> @param dp On entry dp is the the derivative of the function at stp.<br/>
c>           On exit dp is unchanged.
c> 
c> @param brackt On entry brackt specifies if a minimizer has been bracketed.
c>                  Initially brackt must be set to .false.<br/>
c>               On exit brackt specifies if a minimizer has been bracketed.
c>                  When a minimizer is bracketed brackt is set to .true.
c> 
c> @param stpmin On entry stpmin is a lower bound for the step.<br/>
c>               On exit stpmin is unchanged.
c> 
c> @param stpmax On entry stpmax is an upper bound for the step.<br/>
c>               On exit stpmax is unchanged.
      subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,
     +                  stpmin,stpmax)
      logical brackt
      double precision stx,fx,dx,sty,fy,dy,stp,fp,dp,stpmin,stpmax
c
c     MINPACK-1 Project. June 1983
c     Argonne National Laboratory. 
c     Jorge J. More' and David J. Thuente.
c
c     MINPACK-2 Project. October 1993.
c     Argonne National Laboratory and University of Minnesota. 
c     Brett M. Averick and Jorge J. More'.
c
c     **********
      double precision zero,p66,two,three
      parameter(zero=0.0d0,p66=0.66d0,two=2.0d0,three=3.0d0)
   
      double precision gamma,p,q,r,s,sgnd,stpc,stpf,stpq,theta
 
      sgnd = dp*(dx/abs(dx))
 
c     First case: A higher function value. The minimum is bracketed. 
c     If the cubic step is closer to stx than the quadratic step, the 
c     cubic step is taken, otherwise the average of the cubic and 
c     quadratic steps is taken.
 
      if (fp .gt. fx) then
         theta = three*(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))/two)*
     +                                                       (stp - stx)
         if (abs(stpc-stx) .lt. abs(stpq-stx)) then
            stpf = stpc
         else
            stpf = stpc + (stpq - stpc)/two
         endif
         brackt = .true.
 
c     Second case: A lower function value and derivatives of opposite 
c     sign. The minimum is bracketed. If the cubic step is farther from 
c     stp than the secant step, the cubic step is taken, otherwise the 
c     secant step is taken.
 
      else if (sgnd .lt. zero) then
         theta = three*(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
         endif
         brackt = .true.
 
c     Third case: A lower function value, derivatives of the same sign,
c     and the magnitude of the derivative decreases.
 
      else if (abs(dp) .lt. abs(dx)) then
 
c        The cubic step is computed only 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 the 
c        secant step.
 
         theta = three*(fx - fp)/(stp - stx) + dx + dp
         s = max(abs(theta),abs(dx),abs(dp))
 
c        The case gamma = 0 only arises if the cubic does not tend
c        to infinity in the direction of the step.
 
         gamma = s*sqrt(max(zero,(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. zero .and. gamma .ne. zero) then
            stpc = stp + r*(stx - stp)
         else if (stp .gt. stx) then
            stpc = stpmax
         else
            stpc = stpmin
         endif
         stpq = stp + (dp/(dp - dx))*(stx - stp)
 
         if (brackt) then
 
c           A minimizer has been bracketed. If the cubic step is 
c           closer to stp than the secant step, the cubic step is 
c           taken, otherwise the secant step is taken.
 
            if (abs(stpc-stp) .lt. abs(stpq-stp)) then
               stpf = stpc
            else
               stpf = stpq
            endif
            if (stp .gt. stx) then
               stpf = min(stp+p66*(sty-stp),stpf)
            else
               stpf = max(stp+p66*(sty-stp),stpf)
            endif
         else
 
c           A minimizer has not been bracketed. If the cubic step is 
c           farther from stp than the secant step, the cubic step is 
c           taken, otherwise the secant step is taken.
 
            if (abs(stpc-stp) .gt. abs(stpq-stp)) then
               stpf = stpc
            else
               stpf = stpq
            endif
            stpf = min(stpmax,stpf)
            stpf = max(stpmin,stpf)
         endif
 
c     Fourth case: A lower function value, derivatives of the same sign, 
c     and the magnitude of the derivative does not decrease. If the 
c     minimum is not bracketed, the step is either stpmin or stpmax, 
c     otherwise the cubic step is taken.
 
      else
         if (brackt) then
            theta = three*(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
         endif
      endif
 
c     Update the interval which contains a minimizer.
 
      if (fp .gt. fx) then
         sty = stp
         fy = fp
         dy = dp
      else
         if (sgnd .lt. zero) then
            sty = stx
            fy = fx
            dy = dx
         endif
         stx = stp
         fx = fp
         dx = dp
      endif
 
c     Compute the new step.
 
      stp = stpf
 
      return
      end