rootnf Subroutine

   subroutine rootnf(state, callbacks, relerr, abserr)
   !! This subroutine finds the point `ybar = (1, xbar)` on the zero curve of the homotopy
   !! map. It starts with two points `yold = (lambdaold, xold)` and `y = (lambda, x)` such
   !! that `lambdaold < 1 <= lambda` , and alternates between secant estimates of `ybar`
   !! and Newton iteration until convergence.

      use hompack_core, only: root, root_state, qofs
      implicit none

      type(nf_state), intent(inout) :: state
         !! State variables for [[fixpnf]].
      type(hompack_f_callbacks) :: callbacks
         !! User-supplied function and Jacobian evaluation subroutines.
      real(dp), intent(in) :: relerr
         !! Relative convergence tolerance.
         !! Iteration is considered converged when `|y(1) - 1| <= relerr + abserr` and the
         !! Newton correction satisfies `||z|| <= relerr*||x|| + abserr`.
      real(dp), intent(in) :: abserr
         !! Absolute convergence tolerance.
         !! Used together with `relerr` in the convergence criteria.

      real(dp) :: dels, qsout, aerr, rerr, sa, sb, sout
      integer :: judy, jw, lcode, limit, np1
      logical :: bracket
      type(root_state) :: state_root

      rerr = max(relerr, eps64)
      aerr = max(abserr, zero)
      np1 = state%n + 1

      ! The limit on the number of iterations allowed may be changed by changing the
      ! following parameter statement
      limit = 2*(int(abs(log10(aerr + rerr))) + 1)

      associate (ws => state%workspace)

         ws%tz = state%y - state%yold
         dels = norm2(ws%tz)

         ! Using two points and tangents on the homotopy zero curve, construct the Hermite
         ! cubic interpolant q(s). Then use 'root' to find the 's' corresponding to
         ! 'lambda=1'. The two points on the zero curve are always chosen to bracket
         ! 'lambda=1', with the bracketing interval always being [0, dels].
         sa = zero
         sb = dels
         lcode = 1 ! forces initialization of 'root'
         do
            call root(sout, qsout, sa, sb, rerr, aerr, lcode, state_root)
            if (lcode > 0) exit
            qsout = qofs(state%yold(1), state%ypold(1), state%y(1), state%yp(1), dels, sout) - one
         end do

         ! If lambda=1 were bracketed, root cannot fail
         if (lcode > 2) then
            state%iflag = 6
            return
         end if

         ! Calculate 'q(sa)' as the initial point for a Newton iteration
         do jw = 1, np1
            ws%w(jw) = qofs(state%yold(jw), state%ypold(jw), state%y(jw), state%yp(jw), dels, sa)
         end do

         ! Tangent information 'yp' is no longer needed. Hereafter, 'yp' represents the most
         ! recent point which is on the opposite side of the hyperplane 'lambda=1' from 'y'

         ! Prepare for main loop
         state%yp = state%yold

         ! Initialize bracket to indicate that the points 'y' and 'yold' bracket 'lambda=1',
         ! thus 'yold = yp'
         bracket = .true.

         ! Main loop
         do judy = 1, limit

            ! Calculate Newton step at current estimate 'w'
            call tangnf(callbacks, sa, &
                        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

            ! Next point = current point + Newton step
            ws%w = ws%w + ws%tz

            ! Check for convergence
            if ((abs(ws%w(1) - one) <= rerr + aerr) .and. &
                (norm2(ws%tz) <= rerr*norm2(ws%w(2:np1)) + aerr)) then
               state%y = ws%w
               return
            end if

            ! Prepare for next iteration
            if (abs(ws%w(1) - one) <= rerr + aerr) then
               state%ypold = ws%wp
               cycle
            end if

            ! Update 'y' and 'yold'
            state%yold = state%y
            state%y = ws%w

            ! Update 'yp' such that 'yp' is the most recent point opposite of 'lambda=1'
            ! from 'y'. Set bracket=.true. iff 'y' and 'yold' bracket 'lambda=1' so that
            ! yp = yold .
            if ((state%y(1) - one)*(state%yold(1) - one) > 0) then
               bracket = .false.
            else
               bracket = .true.
               state%yp = state%yold
            end if

            ! Compute dels=||y-yp||
            ws%tz = state%y - state%yp
            dels = norm2(ws%tz)

            ! Compute tz for the linear predictor w = y + tz, where tz = sa*(yold-y).
            sa = (one - state%y(1))/(state%yold(1) - state%y(1))
            ws%tz = sa*(state%yold - state%y)

            ! To insure stability, the linear prediction must be no farther from y than
            ! yp is. This is guaranteed if bracket=true. If linear prediction is too far
            ! away, use bracketing points to compute linear prediction.
            if (.not. bracket) then
               if (norm2(ws%tz) > dels) then
                  ! Compute tz = sa*(yp-y)
                  sa = (one - state%y(1))/(state%yp(1) - state%y(1))
                  ws%tz = sa*(state%yp - state%y)
               end if
            end if

            ! Compute estimate w = y + tz  and save old tangent vector.
            ws%w = ws%w + ws%tz
            state%ypold = ws%wp

         end do

         ! The alternating secant estimation and Newton iteration has not converged in
         ! 'limit' iterations
         state%iflag = 6

      end associate

   end subroutine rootnf