check_jac_zero Subroutine

   subroutine check_jac_zero &
      (fcn, &
       n, m, np, q, &
       beta, xplusd, ifixb, ifixx, ldifx, &
       nrow, epsmac, j, lq, iswrtb, &
       tol, d, fd, typj, pvpstp, stp0, pv, &
       diffj, msg, istop, nfev, &
       wrk1, wrk2, wrk6)
   !! Recheck the derivatives in the case where the finite difference derivative disagrees with
   !! the analytic derivative and the analytic derivative is zero.
      ! Adapted from STARPAC subroutine DCKZRO.

      use odrpack_kinds, only: zero, three

      procedure(fcn_t) :: fcn
         !! User supplied subroutine for evaluating the model.
      integer, intent(in) :: n
         !! Number of observations.
      integer, intent(in) :: m
         !! Number of columns of data in the explanatory variable.
      integer, intent(in) :: np
         !! Number of function parameters.
      integer, intent(in) :: q
         !! Number of responses per observation.
      real(wp), intent(inout) :: beta(np)
         !! Function parameters.
      real(wp), intent(inout) :: xplusd(n, m)
         !! Values of `x + delta`.
      integer, intent(in) :: ifixb(np)
         !! Values designating whether the elements of `beta` are fixed at their input values or not.
      integer, intent(in) :: ifixx(ldifx, m)
         !! Values designating whether the elements of `x` are fixed at their input values or not.
      integer, intent(in) :: ldifx
         !! Leading dimension of array `ifixx`.
      integer, intent(in) :: nrow
         !! Row number of the explanatory variable array at which the derivative is to be checked.
      real(wp), intent(in) :: epsmac
         !! Value of machine precision.
      integer, intent(in) :: j
         !! Index of the partial derivative being examined.
      integer, intent(in) :: lq
         !! Response currently being examined.
      logical, intent(in) :: iswrtb
         !! Variable designating whether the derivatives wrt `beta` (`iswrtb = .true.`)
         !! or `delta` (`iswrtb = .false.`) are being checked.
      real(wp), intent(in) :: tol
         !! Agreement tolerance.
      real(wp), intent(in) :: d
         !! Derivative with respect to the `j`-th unknown parameter.
      real(wp), intent(in) :: fd
         !! Forward difference derivative wrt the `j`-th parameter.
      real(wp), intent(in) :: typj
         !! Typical size of the `j`-th unknown `beta` or `delta`.
      real(wp), intent(in) :: pvpstp
         !! Predicted value for row `nrow` of the model using the current parameter estimates
         !! for all but the `j`-th parameter value, which is `beta(j) + stp0`.
      real(wp), intent(in) :: stp0
         !! Initial step size for the finite difference derivative.
      real(wp), intent(in) :: pv
         !! Predicted value from the model for row `nrow`.
      real(wp), intent(out) :: diffj
         !! Relative differences between the user supplied and finite difference derivatives
         !! for the derivative being checked.
      integer, intent(out) :: msg(q, j)
         !! Error checking results.
      integer, intent(out) :: istop
         !! Variable designating whether there are problems computing the function at the
         !! current `beta` and `delta`.
      integer, intent(inout) :: nfev
         !! Number of function evaluations.
      real(wp), intent(out) :: wrk1(n, m, q)
         !! A work array of `(n, m, q)` elements.
      real(wp), intent(out) :: wrk2(n, q)
         !! A work array of `(n, q)` elements.
      real(wp), intent(out) :: wrk6(n, np, q)
         !! A work array of `(n, np, q)` elements.

      ! Local scalars
      real(wp) :: cd, pvmstp

      ! Variable Definitions (alphabetically)
      !  CD:      The central difference derivative wrt the Jth parameter.
      !  PVMSTP:  The predicted value for row NROW of the model using the current parameter
      !           estimates for all but the Jth parameter value, which is BETA(J) - STP0.

      ! Recalculate numerical derivative using central difference and step size of 2*STP0
      if (iswrtb) then
         ! Perform computations for derivatives wrt BETA
         call fpvb(fcn, &
                   n, m, np, q, &
                   beta, xplusd, ifixb, ifixx, ldifx, &
                   nrow, j, lq, -stp0, &
                   istop, nfev, pvmstp, &
                   wrk1, wrk2, wrk6)
      else
         ! Perform computations for derivatives wrt DELTA
         call fpvd(fcn, &
                   n, m, np, q, &
                   beta, xplusd, ifixb, ifixx, ldifx, &
                   nrow, j, lq, -stp0, &
                   istop, nfev, pvmstp, &
                   wrk1, wrk2, wrk6)
      end if

      if (istop /= 0) then
         return
      end if

      cd = (pvpstp - pvmstp)/(2*stp0)
      diffj = min(abs(cd - d), abs(fd - d))

      ! Check for agreement
      if (diffj <= tol*abs(d)) then
         ! Finite difference and analytic derivatives now agree
         if (d == zero) then
            msg(lq, j) = 1
         else
            msg(lq, j) = 0
         end if
      elseif (diffj*typj <= abs(pv*epsmac**(1/three))) then
         ! Derivatives are both close to zero
         msg(lq, j) = 2
      else
         ! Derivatives are not both close to zero
         msg(lq, j) = 3
      end if

   end subroutine check_jac_zero