check_jac_curv Subroutine

   subroutine check_jac_curv &
      (fcn, &
       n, m, np, q, &
       beta, xplusd, ifixb, ifixx, ldifx, &
       eta, tol, nrow, epsmac, j, lq, hc, iswrtb, &
       fd, typj, pvpstp, stp0, &
       pv, d, &
       diffj, msg, istop, nfev, &
       wrk1, wrk2, wrk6)
   !! Check whether high curvature could be the cause of the disagreement between the numerical
   !! and analytic derviatives.
      ! Adapted from STARPAC subroutine DCKCRV.

      use odrpack_kinds, only: one, two, ten

      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`.
      real(wp), intent(in) :: eta
         !! Relative noise in the model.
      real(wp), intent(in) :: tol
         !! Agreement tolerance.
      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.
      real(wp), intent(in) :: hc
         !! Relative step size for central finite differences.
      logical, intent(in) :: iswrtb
         !! Variable designating whether the derivatives wrt `beta` (`iswrtb = .true.`) or
         !! `delta` (`iswrtb = .false.`) are being checked.
      real(wp), intent(out) :: 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(out) :: pvpstp
         !! Predicted value for row `nrow` of the model based on 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 of the model for row `nrow`.
      real(wp), intent(in) :: d
         !! Derivative with respect to the `j`-th unknown parameter.
      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), parameter :: p01 = 0.01_wp
      real(wp) :: curve, pvmcrv, pvpcrv, stp, stpcrv

      ! Variable Definitions (alphabetically)
      !  CURVE:   A measure of the curvature in the model.
      !  PVMCRV:  The predicted value for row  NROW of the model based on the current parameter
      !           estimates for all but the Jth parameter value, which is BETA(J)-STPCRV.
      !  PVPCRV:  The predicted value for row NROW of the model based on the current parameter
      !           estimates for all but the Jth parameter value, which is BETA(J)+STPCRV.
      !  STP:     A step size for the finite difference derivative.
      !  STPCRV:  The step size selected to check for curvature in the model.

      if (iswrtb) then

         ! Perform central difference computations for derivatives wrt BETA
         stpcrv = (hc*typj*sign(one, beta(j)) + beta(j)) - beta(j)
         call fpvb(fcn, &
                   n, m, np, q, &
                   beta, xplusd, ifixb, ifixx, ldifx, &
                   nrow, j, lq, stpcrv, &
                   istop, nfev, pvpcrv, &
                   wrk1, wrk2, wrk6)
         if (istop /= 0) then
            return
         end if
         call fpvb(fcn, &
                   n, m, np, q, &
                   beta, xplusd, ifixb, ifixx, ldifx, &
                   nrow, j, lq, -stpcrv, &
                   istop, nfev, pvmcrv, &
                   wrk1, wrk2, wrk6)
         if (istop /= 0) then
            return
         end if
      else

         ! Perform central difference computations for derivatives wrt DELTA
         stpcrv = (hc*typj*sign(one, xplusd(nrow, j)) + xplusd(nrow, j)) - xplusd(nrow, j)
         call fpvd(fcn, &
                   n, m, np, q, &
                   beta, xplusd, ifixb, ifixx, ldifx, &
                   nrow, j, lq, stpcrv, &
                   istop, nfev, pvpcrv, &
                   wrk1, wrk2, wrk6)
         if (istop /= 0) then
            return
         end if
         call fpvd(fcn, &
                   n, m, np, q, &
                   beta, xplusd, ifixb, ifixx, ldifx, &
                   nrow, j, lq, -stpcrv, &
                   istop, nfev, pvmcrv, &
                   wrk1, wrk2, wrk6)
         if (istop /= 0) then
            return
         end if
      end if

      ! Estimate curvature by second derivative of model
      curve = abs((pvpcrv - pv) + (pvmcrv - pv))/(stpcrv*stpcrv)
      curve = curve + eta*(abs(pvpcrv) + abs(pvmcrv) + two*abs(pv))/(stpcrv**2)

      ! Check if finite precision arithmetic could be the culprit.
      call check_jac_fp(fcn, &
                        n, m, np, q, &
                        beta, xplusd, ifixb, ifixx, ldifx, &
                        eta, tol, nrow, j, lq, iswrtb, &
                        fd, typj, pvpstp, stp0, curve, pv, d, &
                        diffj, msg, istop, nfev, &
                        wrk1, wrk2, wrk6)
      if (istop /= 0) then
         return
      end if
      if (msg(lq, j) == 0) then
         return
      end if

      ! Check if high curvature could be the problem.
      stp = two*max(tol*abs(d)/curve, epsmac)
      if (stp < abs(ten*stp0)) then
         stp = min(stp, p01*abs(stp0))
      end if

      if (iswrtb) then
         ! Perform computations for derivatives wrt BETA
         stp = (stp*sign(one, beta(j)) + beta(j)) - beta(j)
         call fpvb(fcn, &
                   n, m, np, q, &
                   beta, xplusd, ifixb, ifixx, ldifx, &
                   nrow, j, lq, stp, &
                   istop, nfev, pvpstp, &
                   wrk1, wrk2, wrk6)
         if (istop /= 0) then
            return
         end if
      else

         ! Perform computations for derivatives wrt DELTA
         stp = (stp*sign(one, xplusd(nrow, j)) + xplusd(nrow, j)) - xplusd(nrow, j)
         call fpvd(fcn, &
                   n, m, np, q, &
                   beta, xplusd, ifixb, ifixx, ldifx, &
                   nrow, j, lq, stp, &
                   istop, nfev, pvpstp, &
                   wrk1, wrk2, wrk6)
         if (istop /= 0) then
            return
         end if
      end if

      ! Compute the new numerical derivative
      fd = (pvpstp - pv)/stp
      diffj = min(diffj, abs(fd - d)/abs(d))

      ! Check whether the new numerical derivative is ok
      if (abs(fd - d) <= tol*abs(d)) then
         msg(lq, j) = 0

         ! Check if finite precision may be the culprit (fudge factor = 2)
      elseif (abs(stp*(fd - d)) < two*eta*(abs(pv) + abs(pvpstp)) + &
              curve*(epsmac*typj)**2) then
         msg(lq, j) = 5
      end if

   end subroutine check_jac_curv