stepnf Subroutine

   subroutine stepnf(state, callbacks, sspar)
   !! This subroutine takes one step along the zero curve of the homotopy map using a
   !! predictor-corrector algorithm. The predictor uses a Hermite cubic interpolant, and
   !! the corrector returns to the zero curve along the flow normal to the Davidenko flow.
   !! [[stepnf]] also estimates a step size `h` for the next step along the zero curve.
   !! Normally, [[stepnf]] is used indirectly through [[fixpnf]], and should be called
   !! directly only if it is necessary to modify the stepping algorithm's parameters.

      use hompack_core, only: qofs
      implicit none

      type(nf_state), intent(inout) :: state
         !! State variables for [[fixpnf]].
      type(hompack_f_callbacks), intent(in) :: callbacks
         !! User-supplied function and Jacobian evaluation subroutines.
      real(dp), intent(in) :: sspar(:)
         !! Step-size control parameters:
         !! `(lideal, rideal, dideal, hmin, hmax, bmin, bmax, p)`.

      real(dp), parameter :: twou = 2*eps64, fouru = 4*eps64
      real(dp) :: dcalc, hfail, ht, lcalc, rcalc, rholen, temp
      integer :: itnum, j, judy, np1
      logical :: fail

      ! Newton iterations allowed before reducing the step size 'h'
      integer, parameter:: max_newton = 4

      np1 = state%n + 1
      state%crash = .true.

      associate (ws => state%workspace)

         ! The arclength 's' must be nonnegative
         if (state%s < zero) return

         ! If step size is too small, determine an acceptable one
         if (state%h < fouru*(one + state%s)) then
            state%h = fouru*(one + state%s)
            return
         end if

         ! If error tolerances are too small, increase them to acceptable values
         temp = norm2(state%y) + one
         if (0.5_dp*(state%relerr*temp + state%abserr) < twou*temp) then
            if (state%relerr /= zero) then
               state%relerr = fouru*(one + fouru)
               state%abserr = max(state%abserr, zero)
            else
               state%abserr = fouru*temp
            end if
            return
         end if

         ! STARTUP SECTION (FIRST STEP ALONG ZERO CURVE)

         state%crash = .false.
         if (state%start) then
            fail = .false.
            state%start = .false.

            ! Determine suitable initial step size
            state%h = min(state%h, 0.1_dp, sqrt(sqrt(state%relerr*temp + state%abserr)))

            ! Use linear predictor along tangent direction to start Newton iteration
            state%ypold(1) = one
            state%ypold(2:np1) = zero
            call tangnf(callbacks, state%s, state%n, state%y, state%yp, state%ypold, state%a, &
                        state%nfe, state%iflag, &
                        ws%qr, ws%alpha, ws%tz, ws%pivot)

            if (state%iflag > 0) return

            do

               ws%w = state%y + state%h*state%yp
               ws%z0 = ws%w

               do judy = 1, max_newton

                  ! Calculate the Newton step 'tz' at the current point 'w'
                  rholen = -one
                  call tangnf(callbacks, rholen, &
                              state%n, ws%w, ws%wp, state%ypold, state%a, &
                              state%nfe, state%iflag, &
                              ws%qr, ws%alpha, ws%tz, ws%pivot)
                  if (state%iflag > 0) return

                  ! Take Newton step and check convergence
                  ws%w = ws%w + ws%tz
                  itnum = judy

                  ! Compute quantities used for optimal step size estimation
                  if (judy == 1) then
                     lcalc = norm2(ws%tz)
                     rcalc = rholen
                     ws%z1 = ws%w
                  else if (judy == 2) then
                     lcalc = norm2(ws%tz)/lcalc
                     rcalc = rholen/rcalc
                  end if

                  ! Go to mop-up section after convergence
                  if (norm2(ws%tz) <= state%relerr*norm2(ws%w) + state%abserr) then
                     go to 600
                  end if

               end do

               ! No convergence in 'max_newton' iterations. Reduce h and try again.
               if (state%h <= fouru*(one + state%s)) then
                  state%iflag = 6
                  return
               end if
               state%h = state%h/2

            end do

         end if

         ! PREDICTOR SECTION

         fail = .false.

         do

            ! Compute point predicted by Hermite interpolant. Use step size 'h' computed on
            ! last call to 'stepnf'.
            do j = 1, np1
               ws%w(j) = qofs(state%yold(j), state%ypold(j), state%y(j), state%yp(j), &
                              state%hold, state%hold + state%h)
            end do

            ws%z0 = ws%w

            ! CORRECTOR SECTION
            do judy = 1, max_newton

               ! Calculate the Newton step 'tz' at the current point 'w'
               rholen = -one
               call tangnf(callbacks, rholen, &
                           state%n, ws%w, ws%wp, state%yp, state%a, state%nfe, &
                           state%iflag, ws%qr, ws%alpha, ws%tz, ws%pivot)
               if (state%iflag > 0) return

               ! Take Newton step and check convergence
               ws%w = ws%w + ws%tz
               itnum = judy

               ! Compute quantities used for optimal step size estimation.
               if (judy == 1) then
                  lcalc = norm2(ws%tz)
                  rcalc = rholen
                  ws%z1 = ws%w
               else if (judy == 2) then
                  lcalc = norm2(ws%tz)/lcalc
                  rcalc = rholen/rcalc
               end if

               ! Go to mop-up section after convergence
               if (norm2(ws%tz) <= state%relerr*norm2(ws%w) + state%abserr) then
                  go to 600
               end if

            end do

            ! No convergence in 'max_newton' iterations. Record failure at calculated 'h'
            ! Save this step size, reduce 'h' and try again.

            fail = .true.
            hfail = state%h

            if (state%h <= fouru*(one + state%s)) then
               state%iflag = 6
               return
            end if

            state%h = state%h/2

         end do

         ! MOP-UP SECTION

         ! 'yold'  and  'y'  always contain the last two points found on the zero curve of
         ! the homotopy map. 'ypold' and 'yp' contain the tangent vectors to the zero curve
         ! at  'yold'  and  'y' , respectively.
600      state%ypold = state%yp
         state%yold = state%y
         state%y = ws%w
         state%yp = ws%wp
         ws%w = state%y - state%yold

         ! Update arc length
         state%hold = norm2(ws%w)
         state%s = state%s + state%hold

         ! OPTIMAL STEP SIZE ESTIMATION SECTION

         ! Calculate the distance factor 'dcalc'
         ws%tz = ws%z0 - state%y
         ws%w = ws%z1 - state%y
         dcalc = norm2(ws%tz)
         if (dcalc /= zero) dcalc = norm2(ws%w)/dcalc

         ! The optimal step size hbar is defined by
         !
         !   ht = hold * [min(lideal/lcalc, rideal/rcalc, dideal/dcalc)]**(1/p)
         !
         !   hbar = min [ max(ht, bmin*hold, hmin), bmax*hold, hmax ]

         ! If convergence had occurred after 1 iteration, set the contraction factor 'lcalc'
         ! to zero.
         if (itnum == 1) lcalc = zero

         ! Formula for optimal step size
         if (lcalc + rcalc + dcalc == zero) then
            ht = sspar(7)*state%hold
         else
            ht = (one/max(lcalc/sspar(1), rcalc/sspar(2), &
                          dcalc/sspar(3)))**(one/sspar(8))*state%hold
         end if

         ! 'ht' contains the estimated optimal step size. Put it within reasonable bounds
         state%h = min( &
                   max(ht, sspar(6)*state%hold, sspar(4)), &
                   sspar(7)*state%hold, &
                   sspar(5) &
                   )

         if (itnum == 1) then
            ! If convergence had occurred after 1 iteration, don't decrease 'h'.
            state%h = max(state%h, state%hold)
         else if (itnum == max_newton) then
            ! If convergence required the maximum 'max_newton' iterations, don't increase 'h'.
            state%h = min(state%h, state%hold)
         end if

         ! If convergence did not occur in 'max_newton' iterations for a particular
         ! 'h = hfail', don't choose the new step size larger than 'hfail'.
         if (fail) state%h = min(state%h, hfail)

      end associate

   end subroutine stepnf