tangnf Subroutine

   subroutine tangnf( &
      callbacks, rholen, n, y, yp, ypold, a, nfe, iflag, qr, alpha, tz, pivot)
   !! This subroutine builds the Jacobian matrix of the homotopy map, computes a QR
   !! decomposition of that matrix, and then calculates the (unit) tangent vector and the
   !! Newton step.

      use lapack_interfaces, only: dgeqpf, dormqr
      implicit none

      type(hompack_f_callbacks) :: callbacks
         !! User-supplied function and Jacobian evaluation subroutines.
      real(dp), intent(inout) :: rholen
         !! Controls computation of the homotopy residual norm.
         !!
         !! On input:
         !! * `rholen < 0` requests evaluation of the homotopy map norm.
         !! * `rholen >= 0` suppresses norm evaluation.
         !!
         !! On output, if `rholen < 0` on entry:
         !! `rholen = ||rho(a, lambda, x)||`. Otherwise the value is unchanged.
      integer, intent(in) :: n
         !! Problem dimension.
      real(dp), intent(in) :: y(:)
         !! Current point on the homotopy zero curve. `Shape: (n+1)`.
         !! Contains `(lambda, x)`.
      real(dp), intent(out) :: yp(:)
         !! Unit tangent vector to the zero curve at `y`. `Shape: (n+1)`.
         !! Represents `dy/ds`, where `s` is arc length along the zero curve.
      real(dp), intent(in) :: ypold(:)
         !! Unit tangent vector at the previous point on the zero curve. `Shape: (n+1)`.
         !! Used to determine a consistent orientation for the newly computed tangent
         !! vector.
      real(dp), intent(in) :: a(:)
         !! Parameter vector used in the homotopy map.
      integer, intent(inout) :: nfe
         !! Number of homotopy/Jacobian evaluations performed.
         !! Incremented by one on every successful call.
      integer, intent(inout) :: iflag
         !! Problem type and return status flag.
         !!
         !! On input:
         !! * `0`  : fixed-point problem, `F(x) = x`.
         !! * `-1` : zero-finding problem, `F(x) = 0`.
         !! * `-2` : general homotopy curve-tracking problem.
         !!
         !! On output:
         !! * unchanged (`0`, `-1`, or `-2`) on normal return.
         !! * `4` : Jacobian matrix lost full rank (`rank < n`); iteration not completed.
      real(dp), intent(inout) :: qr(:, :)
         !! Workspace containing the Jacobian matrix and its QR factorization.
         !! `Shape: (n, n+2)`.
      real(dp), intent(inout) :: alpha(:)
         !! Workspace array used during QR factorization and related linear algebra
         !! operations. `Shape: (3*n+3)`.
      real(dp), intent(out) :: tz(:)
         !! Newton correction vector. `Shape: (n+1)`.
         !! Equal to the negative pseudoinverse of the homotopy Jacobian applied to the
         !! homotopy residual.
      integer, intent(inout) :: pivot(:)
         !! Pivot indices produced by the QR factorization. `Shape: (n+1)`.

      real(dp) :: lambda, sigma, ypnorm
      integer :: i, j, k, kp1, np1, np2, info

      lambda = y(1)
      np1 = n + 1
      np2 = n + 2
      nfe = nfe + 1

      ! Compute the jacobian matrix, store it and homotopy map in QR
      if (iflag == -2) then

         ! QR = [ d rho(a,lambda,x)/d lambda , d rho(a,lambda,x)/dx , rho(a,lambda,x) ]
         do k = 1, np1
            call callbacks%rhojac(a, lambda, y(2:np1), qr(:, k), k, callbacks%data)
         end do

         call callbacks%rho(a, lambda, y(2:np1), qr(:, np2), callbacks%data)

      else

         call callbacks%f(y(2:np1), tz(1:n), callbacks%data)

         if (iflag == 0) then

            ! QR = [ a - f(x), i - lambdadf(x), x - a + lambda(a - f(x)) ]
            qr(:, 1) = a - tz(1:n)
            qr(:, np2) = y(2:np1) - a + lambda*qr(:, 1)
            do k = 1, n
               call callbacks%fjac(y(2:np1), tz(1:n), k, callbacks%data)
               kp1 = k + 1
               qr(:, kp1) = -lambda*tz(1:n)
               qr(k, kp1) = one + qr(k, kp1)
            end do

         else

            ! QR = [ f(x) - x + a, lambda*df(x) + (1 - lambda)i , x - a + lambda(f(x) - x + a) ]
            qr(:, 1) = tz(1:n) - y(2:np1) + a
            qr(:, np2) = y(2:np1) - a + lambda*qr(:, 1)
            do k = 1, n
               call callbacks%fjac(y(2:np1), tz(1:n), k, callbacks%data)
               kp1 = k + 1
               qr(:, kp1) = lambda*tz(1:n)
               qr(k, kp1) = one - lambda + qr(k, kp1)
            end do

         end if

      end if

      ! Compute the norm of the homotopy map if it was requested
      if (rholen < zero) rholen = norm2(qr(:, np2))

      ! Reduce the Jacobian matrix to upper triangular form
      pivot = 0
      call dgeqpf(n, np1, qr, n, pivot, yp, alpha, info)

      if (abs(qr(n, n)) <= abs(qr(1, 1))*eps64) then
         iflag = 4
         return
      end if

      call dormqr('L', 'T', n, 1, n, qr, n, yp, qr(:, np2), n, alpha, 3*n + 3, info)

      do i = 1, n
         alpha(i) = qr(i, i)
      end do

      ! Compute kernel of Jacobian, which specifies yp=dy/ds.
      tz(np1) = one
      do i = n, 1, -1
         j = i + 1
         tz(i) = -dot_product(qr(i, j:np1), tz(j:np1))/alpha(i)
      end do
      ypnorm = norm2(tz)
      yp(pivot) = tz/ypnorm
      if (dot_product(yp, ypold) < zero) yp = -yp

      ! 'yp' is the unit tangent vector in the correct direction
      !
      ! Compute the minimum norm solution of [d rho(y(s))] v = -rho(y(s)).
      ! 'v' is given by 'p - (p,q)q', where 'p' is any solution of [d rho] v = -rho
      ! and 'q' is a unit vector in the kernel of [d rho].
      alpha(np1) = one
      do i = n, 1, -1
         j = i + 1
         alpha(i) = -(dot_product(qr(i, j:np1), alpha(j:np1)) + qr(i, np2))/alpha(i)
      end do
      tz(pivot) = alpha(1:np1)

      ! 'tz' now contains a particular solution 'p', and 'yp' contains a vector 'q'
      ! in the kernel (the tangent)
      sigma = dot_product(tz, yp)

      ! 'tz' is the Newton step from the current point y(s) = (lambda(s), x(s)).
      tz = tz - sigma*yp

   end subroutine tangnf