dcstep Subroutine

      subroutine dcstep(Stx,Fx,Dx,Sty,Fy,Dy,Stp,Fp,Dp,Brackt,Stpmin, &
                        Stpmax)
      implicit none

      logical,intent(inout) :: Brackt  !! On entry `brackt` specifies if a minimizer has been bracketed.
                                       !! Initially `brackt` must be set to .false.
                                       !! On exit `brackt` specifies if a minimizer has been bracketed.
                                       !! When a minimizer is bracketed `brackt` is set to .true.
      real(wp),intent(inout) :: Stx !! On entry `stx` is the best step obtained so far and is an
                                    !! endpoint of the interval that contains the minimizer.
                                    !! On exit `stx is the updated best step.
      real(wp),intent(inout) :: Fx !! On entry `fx` is the function at `stx`.
                                   !! On exit `fx` is the function at `stx`.
      real(wp),intent(inout) :: Dx !! On entry `dx` is the derivative of the function at
                                   !! `stx`. The derivative must be negative in the direction of
                                   !! the step, that is, `dx` and `stp - stx` must have opposite
                                   !! signs.
                                   !! On exit `dx` is the derivative of the function at `stx`.
      real(wp),intent(inout) :: Sty !! On entry `sty` is the second endpoint of the interval that contains the minimizer.
                                    !! On exit `sty` is the updated endpoint of the interval that contains the minimizer.
      real(wp),intent(inout) :: Fy !! On entry `fy` is the function at `sty`.
                                   !! On exit `fy` is the function at `sty`.
      real(wp),intent(inout) :: Dy !! On entry `dy` is the derivative of the function at `sty`.
                                   !! On exit `dy` is the derivative of the function at the exit `sty`.
      real(wp),intent(inout) :: Stp !! On entry `stp` is the current step. If `brackt` is set to .true.
                                    !! then on input `stp` must be between `stx` and `sty`.
                                    !! On exit `stp` is a new trial step.
      real(wp),intent(in) :: Fp !! the function at `stp`.
      real(wp),intent(in) :: Dp !! the derivative of the function at `stp`.
      real(wp),intent(in) :: Stpmin !! a lower bound for the step.
      real(wp),intent(in) :: Stpmax !! an upper bound for the step.

      real(wp),parameter :: p66 = 0.66_wp

      real(wp) :: gamma , p , q , r , s , sgnd , stpc , stpf , &
                  stpq , theta

      sgnd = Dp*(Dx/abs(Dx))

      ! First case: A higher function value. The minimum is bracketed.
      ! If the cubic step is closer to stx than the quadratic step, the
      ! cubic step is taken, otherwise the average of the cubic and
      ! quadratic steps is taken.

      if ( Fp>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<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)<abs(stpq-Stx) ) then
            stpf = stpc
         else
            stpf = stpc + (stpq-stpc)/two
         endif
         Brackt = .true.

         ! Second case: A lower function value and derivatives of opposite
         ! sign. The minimum is bracketed. If the cubic step is farther from
         ! stp than the secant step, the cubic step is taken, otherwise the
         ! secant step is taken.

      else if ( sgnd<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>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)>abs(stpq-Stp) ) then
            stpf = stpc
         else
            stpf = stpq
         endif
         Brackt = .true.

         ! Third case: A lower function value, derivatives of the same sign,
         ! and the magnitude of the derivative decreases.

      else if ( abs(Dp)<abs(Dx) ) then

         ! The cubic step is computed only if the cubic tends to infinity
         ! in the direction of the step or if the minimum of the cubic
         ! is beyond stp. Otherwise the cubic step is defined to be the
         ! secant step.

         theta = three*(Fx-Fp)/(Stp-Stx) + Dx + Dp
         s = max(abs(theta),abs(Dx),abs(Dp))

         ! The case gamma = 0 only arises if the cubic does not tend
         ! to infinity in the direction of the step.

         gamma = s*sqrt(max(zero,(theta/s)**2-(Dx/s)*(Dp/s)))
         if ( Stp>Stx ) gamma = -gamma
         p = (gamma-Dp) + theta
         q = (gamma+(Dx-Dp)) + gamma
         r = p/q
         if ( r<zero .and. gamma/=zero ) then
            stpc = Stp + r*(Stx-Stp)
         else if ( Stp>Stx ) then
            stpc = Stpmax
         else
            stpc = Stpmin
         endif
         stpq = Stp + (Dp/(Dp-Dx))*(Stx-Stp)

         if ( Brackt ) then

            ! A minimizer has been bracketed. If the cubic step is
            ! closer to stp than the secant step, the cubic step is
            ! taken, otherwise the secant step is taken.

            if ( abs(stpc-Stp)<abs(stpq-Stp) ) then
               stpf = stpc
            else
               stpf = stpq
            endif
            if ( Stp>Stx ) then
               stpf = min(Stp+p66*(Sty-Stp),stpf)
            else
               stpf = max(Stp+p66*(Sty-Stp),stpf)
            endif
         else

            ! A minimizer has not been bracketed. If the cubic step is
            ! farther from stp than the secant step, the cubic step is
            ! taken, otherwise the secant step is taken.

            if ( abs(stpc-Stp)>abs(stpq-Stp) ) then
               stpf = stpc
            else
               stpf = stpq
            endif
            stpf = min(Stpmax,stpf)
            stpf = max(Stpmin,stpf)
         endif

         ! Fourth case: A lower function value, derivatives of the same sign,
         ! and the magnitude of the derivative does not decrease. If the
         ! minimum is not bracketed, the step is either stpmin or stpmax,
         ! 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>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>Stx ) then
            stpf = Stpmax
         else
            stpf = Stpmin
         endif
      endif

      ! Update the interval which contains a minimizer.

      if ( Fp>Fx ) then
         Sty = Stp
         Fy = Fp
         Dy = Dp
      else
         if ( sgnd<zero ) then
            Sty = Stx
            Fy = Fx
            Dy = Dx
         endif
         Stx = Stp
         Fx = Fp
         Dx = Dp
      endif

      ! Compute the new step.

      Stp = stpf

      end subroutine dcstep