check_jac_value Subroutine

   subroutine check_jac_value &
      (fcn, &
       n, m, np, q, &
       beta, xplusd, ifixb, ifixx, ldifx, &
       eta, tol, nrow, epsmac, j, lq, typj, h0, hc0, &
       iswrtb, pv, d, &
       diffj, msg1, msg, istop, nfev, &
       wrk1, wrk2, wrk6, interval)
   !! Check user supplied analytic derivatives against numerical derivatives.
      ! Adapted from STARPAC subroutine DCKMN.

      use odrpack_kinds, only: zero, one, two, three, ten, hundred

      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) :: typj
         !! Typical size of the `j`-th unknown `beta` or `delta`.
      real(wp), intent(in) :: h0
         !! Initial step size for the finite difference derivative.
      real(wp), intent(in) :: hc0
         !! 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(in) :: pv
         !! Predicted value 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) :: msg1
         !! First set of error checking results.
      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.
      integer, intent(in) :: interval(np)
         !! Specifies which checks can be performed when checking derivatives based on the
         !! interval of the bound constraints.

      ! Local scalars
      real(wp), parameter :: p01 = 0.01_wp, p1 = 0.1_wp
      real(wp), parameter :: big = 1.0E19_wp, tol2 = 5.0E-2_wp
      real(wp) :: fd, h, hc, h1, hc1, pvpstp, stp0
      integer :: i

      ! Variable Definitions (alphabetically)
      !  BIG:      A big value, used to initialize DIFFJ.
      !  FD:       The forward difference derivative wrt the Jth parameter.
      !  H:        The relative step size for forward differences.
      !  H1:       The default relative step size for forward differences.
      !  HC:       The relative step size for central differences.
      !  HC1:      The default relative step size for central differences.
      !  PVPSTP:   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.
      !  STP0:     The initial step size for the finite difference derivative.
      !  TOL2:     A minimum agreement tolerance.

      ! Calculate the Jth partial derivative using forward difference quotients and decide if
      ! it agrees with user supplied values

      h1 = sqrt(eta)
      hc1 = eta**(1/three)

      msg(lq, j) = 7
      diffj = big

      do i = 1, 3

         if (i == 1) then
            ! Try initial relative step size
            h = h0
            hc = hc0
         elseif (i == 2) then
            ! Try larger relative step size
            h = max(ten*h1, min(hundred*h0, one))
            hc = max(ten*hc1, min(hundred*hc0, one))
         elseif (i == 3) then
            ! Try smaller relative step size
            h = min(p1*h1, max(p01*h, two*epsmac))
            hc = min(p1*hc1, max(p01*hc, two*epsmac))
         end if

         if (iswrtb) then
            ! Perform computations for derivatives wrt BETA
            stp0 = (h*typj*sign(one, beta(j)) + beta(j)) - beta(j)
            call fpvb(fcn, &
                      n, m, np, q, &
                      beta, xplusd, ifixb, ifixx, ldifx, &
                      nrow, j, lq, stp0, &
                      istop, nfev, pvpstp, &
                      wrk1, wrk2, wrk6)
         else
            ! Perform computations for derivatives wrt DELTA
            stp0 = (h*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, stp0, &
                      istop, nfev, pvpstp, &
                      wrk1, wrk2, wrk6)
         end if
         if (istop /= 0) then
            return
         end if

         fd = (pvpstp - pv)/stp0

         ! Check for agreement

         ! Numerical and analytic derivatives agree
         if (abs(fd - d) <= tol*abs(d)) then

            ! Set relative difference for derivative checking report
            if ((d == zero) .or. (fd == zero)) then
               diffj = abs(fd - d)
            else
               diffj = abs(fd - d)/abs(d)
            end if

            ! Set message flag
            if (d == zero) then
               ! JTH analytic and numerical derivatives are both zero.
               msg(lq, j) = 1

            else
               ! JTH analytic and numerical derivatives are both nonzero.
               msg(lq, j) = 0
            end if

         else
            ! Numerical and analytic derivatives disagree.  Check why
            if ((d == zero) .or. (fd == zero)) then
               if (interval(j) >= 10 .or. .not. iswrtb) then
                  call 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)
               else
                  msg(lq, j) = 8
               end if
            else
               if (interval(j) >= 100 .or. .not. iswrtb) then
                  call 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)
               else
                  msg(lq, j) = 8
               end if
            end if

            if (msg(lq, j) <= 2) then
               exit
            end if

         end if
      end do

      ! Set summary flag to indicate questionable results
      if ((msg(lq, j) >= 7) .and. (diffj <= tol2)) then
         msg(lq, j) = 6
      end if
      if ((msg(lq, j) >= 1) .and. (msg(lq, j) <= 6)) then
         msg1 = max(msg1, 1)
      elseif (msg(lq, j) >= 7) then
         msg1 = 2
      end if

   end subroutine check_jac_value