fixpnf Subroutine

   subroutine fixpnf(state, callbacks, n, x0, iflag, config, a)
   !! This subroutine finds a fixed point or zero of the N-dimensional vector function
   !! \( F(x) \), or tracks a zero curve of a general homotopy map \( \rho(a,\lambda,x) \).
   !!
   !! For the fixed-point problem, \( F(x) \) is assumed to be a \(C^2\) map of some ball
   !! into itself. The equation \( x = F(x) \) is solved by following the zero curve of
   !! the homotopy map
   !!
   !! $$ \rho_a(\lambda, x) = \lambda (x - F(x)) + (1 - \lambda) (x - a) $$
   !!
   !! starting from \( \lambda=0, x=a \). The curve is parameterized by arc length \(s\),
   !! and is followed by solving the ordinary differential equation \( d \rho/ d s = 0 \)
   !! for \( y(s) = (\lambda(s), x(s)) \) using a Hermite cubic predictor and a corrector
   !! which returns to the zero curve along the flow normal to the Davidenko flow (which
   !! consists of the integral curves of \( d \rho/ d s \)).
   !!
   !! For the zero-finding problem, \( F(x) \) is assumed to be a \(C^2\) map such that
   !! for some \(r > 0\),  \(x F(x) \ge 0\) whenever \( ||x|| = r \). The equation
   !! \( F(x) = 0 \) is solved by following the zero curve of the homotopy map
   !!
   !! $$ \rho_a(\lambda, x) = \lambda F(x) + (1 - \lambda) (x - a) $$
   !!
   !! emanating from \( \lambda = 0 \), \(x = a \). Parameter \(a\)  must be an interior
   !! point of the above mentioned balls.
   !!
   !! For the curve tracking problem, \( \rho(a,\lambda,x) \) is assumed to be a \(C^2\)
   !! map from \( E^m \times [0,1) \times E^n \) into \( E^n \), which for almost all
   !! parameter vectors \( a \) in some nonempty open subset of \( E^m \) satisfies
   !!
   !!  $$\operatorname{rank} \left[ \frac{\partial \rho(a,\lambda,x)}{\partial \lambda}, \frac{\partial \rho(a,\lambda,x)}{\partial x} \right] = N$$
   !!
   !! for all points \( (\lambda,x) \) such that \( \rho(a,\lambda,x)=0 \). It is further
   !! assumed that
   !!
   !! $$\operatorname{rank} \left[ \frac{\partial \rho(a,0,x_0)}{\partial x} \right] = N $$
   !!
   !! With \(a\) fixed, the zero curve of \( \rho(a,\lambda,x) \) emanating from
   !! \( \lambda=0 \), \( x=x_0 \) is tracked until \( \lambda=1 \) by solving the
   !! ordinary differential equation \( d\rho(a,\lambda(s),x(s))/ds = 0 \) for
   !! \( y(s) = (\lambda(s), x(s)) \), where \(s\) is the arc length along the zero curve.
   !! Also the homotopy map \( \rho(a,\lambda,x) \) is assumed to be constructed such that
   !! \( d\lambda(0) / ds > 0 \).
   !!
   !! For the fixed point and zero finding problems, the user must supply a subroutine
   !! `f(x,v)` which evaluates \(F(x)\) at \(x\) and returns the vector \(F(x)\) in `v`,
   !! and a subroutine `fjac(x,v,k)` which returns in `v` the `k`th column of the Jacobian
   !! matrix of \(F(x)\) evaluated at \(x\).
   !!
   !! For the curve tracking problem, the user must supply a subroutine `rho(a,lambda,x,v)`
   !! which evaluates the homotopy map \(\rho\) at \( (a,\lambda,x) \) and returns the
   !! corresponding vector in `v`, and a subroutine `rhojac(a,lambda,x,v,k)` which
   !! returns in `v` the `k`th column of the \( N \times (N+1) \) Jacobian matrix
   !! \( [\partial \rho/\partial \lambda, \partial \rho/\partial x] \) evaluated
   !! at \( (a,\lambda,x) \).

      implicit none

      type(nf_state), intent(inout) :: state
         !! Preallocated state for [[fixpnf]]. Set to start conditions on the first call, and
         !! updated on subsequent calls when `iflag=2` or `iflag=3`.
      type(hompack_f_callbacks), intent(in) :: callbacks
         !! User-supplied function and Jacobian evaluation subroutines.
      integer, intent(in) :: n
         !! Problem dimension.
      real(dp), intent(in) :: x0(:)
         !! Initial point for the solver. For fixed-point and zero-finding problems, this is
         !! identical to the initial guess \(a\). For curve-tracking problems, this is the
         !! initial solution at \(\lambda=0\).
      integer, intent(inout) :: iflag
         !! Problem type and status flag.
         !!
         !! On input:
         !! * `0`  : solve `x = f(x)` (first call).
         !! * `-1` : solve `f(x) = 0` (first call).
         !! * `-2` : track a zero curve of `rho(a,lambda,x) = 0` (first call).
         !!
         !! On output:
         !! * `1` : normal return.
         !! * `2` : requested tolerances cannot be achieved; tolerances have been
         !!         increased and the routine may be called again.
         !! * `3` : iteration limit reached; call again to continue.
         !! * `4` : Jacobian matrix lost full rank.
         !! * `5` : tracking algorithm lost the zero curve and is not making progress; the
         !!         tolerances `arc?e` and `ans?e` were too lenient. Retry with smaller
         !!         tolerances.
         !! * `6` : normal flow Newton iteration failed to converge; tolerances `ans?e` may be
         !!         too stringent.
         !! * `7` : illegal input parameters.
         !! * `8` : memory allocation failure.
      type(nf_config), intent(inout) :: config
         !! Configuration parameters for the solver. Optional settings will be set to default
         !! values. Tolerances may be updated (increased) on output if `iflag=2`.
      real(dp), intent(in), optional :: a(:)
         !! Parameter vector `a` for curve-tracking problems. Ignored for fixed-point and
         !! zero-finding problems.

      integer :: iter, np1, dima

      np1 = n + 1

      ! INPUT CHECKS (intentionally light; extensive in the callers)

      ! Validate problem dimension
      if (n < 1 .or. size(x0) /= n .or. state%n /= n) then
         iflag = fixpnf_input_error
         return
      end if

      ! Validate 'iflag' and launch job
      if (iflag == 0 .or. iflag == -1) then
         dima = n
         goto 20
      else if (iflag == -2) then
         if (.not. present(a)) then
            iflag = fixpnf_input_error
            return
         end if
         dima = size(a)
         goto 20
      else if (iflag == 2) then
         goto 120
      else if (iflag == 3) then
         goto 90
      else
         iflag = fixpnf_input_error
         return
      end if

      ! FIRST CALL

      ! Set state variables to starting values
20    state%abserr = zero
      state%relerr = zero
      state%curtol = zero
      state%h = 0.1_dp
      state%hold = zero
      state%s = zero
      state%nfe = 0
      state%start = .true.
      state%crash = .false.

      state%iflag = iflag

      state%y = zero
      state%yold = zero
      state%yp = zero
      state%ypold = zero

      state%y(2:np1) = x0
      state%yp(1) = one
      state%ypold(1) = one

      ! Load 'a'
      if (state%iflag == 0 .or. state%iflag == -1) then
         state%a = x0
      else
         state%a = a
      end if

      ! Default arc tolerances
      if (config%arcre <= zero) config%arcre = sqrt(config%ansre)/2
      if (config%arcae <= zero) config%arcae = sqrt(config%ansae)/2

      ! Set optimal step size estimation parameters
      ! 'z[k]' denote the Newton iterates along the flow normal to the Davidenko flow
      ! and 'y' their limit

      associate (sspar => config%sspar)
         ! Ideal contraction factor: ||z[2] - z[1]|| / ||z[1] - z[0]||
         if (sspar(1) <= zero) sspar(1) = 0.5_dp
         ! Ideal residual factor: ||rho(A, z[1])|| / ||rho(a, z[0])||
         if (sspar(2) <= zero) sspar(2) = 0.01_dp
         ! Ideal distance factor: ||z[1] - y|| / ||z[0] - y||
         if (sspar(3) <= zero) sspar(3) = 0.5_dp
         ! Minimum step size 'hmin'
         if (sspar(4) <= zero) sspar(4) = (sqrt(real(n + 1, dp)) + 4.0_dp)*eps64
         ! Maximum step size 'hmax'
         if (sspar(5) <= zero) sspar(5) = one
         ! Minimum step size reduction factor 'bmin'
         if (sspar(6) <= zero) sspar(6) = 0.1_dp
         ! Maximum step size expansion factor 'bmax'
         if (sspar(7) <= zero) sspar(7) = 3.0_dp
         ! Assumed operating order 'p'
         if (sspar(8) <= zero) sspar(8) = 2.0_dp
      end associate

      ! COMMON PART FOR FIRST CALL AND RESTARTS

      ! Set default iteration limit
90    state%limit = config%max_steps

      ! Main loop
120   do iter = 1, state%limit

         ! Tracking algorithm lost the zero curve
         if (state%y(1) < zero) then
            iflag = 5
            return
         end if

         ! Set different error tolerance if the trajectory y(s) has any high curvature
         ! components
         state%curtol = config%cursw*state%hold
         state%relerr = config%arcre
         state%abserr = config%arcae
         if (any(abs(state%yp - state%ypold) > state%curtol)) then
            state%relerr = config%ansre
            state%abserr = config%ansae
         end if

         ! Take a step along the curve
         call stepnf(state, callbacks, config%sspar)

         ! Print latest point on curve if requested
         if (config%lunit /= 0) then
            write (config%lunit, 217) iter, state%nfe, state%s, state%y(1), state%y(2:np1)
217         format(/' STEP', i5, 3x, 'NFE =', i5, 3x, 'ARC LENGTH =', f9.4, 3x, &
                    'LAMBDA =', f7.4, 5x, 'X VECTOR:'/(1x, 6es12.4))
         end if

         ! Check if the step was successful
         if (state%iflag > 0) then
            iflag = state%iflag
            return
         end if

         if (state%crash) then
            ! Return code for error tolerance too small
            iflag = 2
            ! Change error tolerances
            if (config%arcre < state%relerr) config%arcre = state%relerr
            if (config%ansre < state%relerr) config%ansre = state%relerr
            if (config%arcae < state%abserr) config%arcae = state%abserr
            if (config%ansae < state%abserr) config%ansae = state%abserr
            ! Change limit on number of iterations
            state%limit = state%limit - iter
            return
         end if

         ! Use Hermite cubic interpolation and Newton iteration to get the answer at lambda=1
         if (state%y(1) >= one) then

            associate (ws => state%workspace)

               ! Save 'yold' for arc length calculation later
               ws%z0 = state%yold
               call rootnf(state, callbacks, config%ansre, config%ansae)
               iflag = 1

               ! Set error flag if 'rootnf' could not get the point on the zero curve at
               ! lambda=1
               if (state%iflag > 0) iflag = state%iflag

               ! Calculate final arc length
               ws%w = state%y - ws%z0
               state%s = state%s - state%hold + norm2(ws%w)
               return

            end associate

         end if

         ! For polynomial systems and the 'polsys1h' homotopy map, dlambda/ds>= 0
         ! necessarily; this condition is enforced here
         if (config%polsys) then
            if (state%yp(1) < zero) then
               ! Reverse tangent direction so D LAMBDA/DS = YP(1) > 0
               state%yp = -state%yp
               state%ypold = state%yp
               ! Force 'stepnf' to use the linear predictor for the next step only
               state%start = .true.
            end if
         end if

      end do

      ! Lambda has not reached 1 in 'limitd' steps
      iflag = 3

   end subroutine fixpnf