check_jac Subroutine

   subroutine check_jac &
      (fcn, &
       n, m, np, q, &
       beta, betaj, xplusd, &
       ifixb, ifixx, ldifx, stpb, stpd, ldstpd, &
       ssf, tt, ldtt, &
       eta, neta, ntol, nrow, isodr, epsmac, &
       pv0i, fjacb, fjacd, &
       msgb, msgd, diff, istop, nfev, njev, &
       wrk1, wrk2, wrk6, &
       interval)
   !! Driver routine for the derivative checking process.
      ! Adapted from STARPAC subroutine DCKCNT.

      use odrpack_kinds, only: zero, one, p5 => half

      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) :: betaj(np)
         !! Function parameters offset such that steps don't cross bounds.
      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) :: stpb(np)
         !! Step size for finite difference derivatives wrt `beta`.
      real(wp), intent(in) :: stpd(ldstpd, m)
         !! Step size for finite difference derivatives wrt `delta`.
      integer, intent(in) :: ldstpd
         !! Leading dimension of array `stpd`.
      real(wp), intent(in) :: ssf(np)
         !! Scaling values used for `beta`.
      real(wp), intent(in) :: tt(ldtt, m)
         !! Scaling values used for `delta`.
      integer, intent(in) :: ldtt
         !! Leading dimension of array `tt`.
      real(wp), intent(in) :: eta
         !! Relative noise in the function results.
      integer, intent(in) :: neta
         !! Number of reliable digits in the model results.
      integer, intent(out) :: ntol
         !! Number of digits of agreement required between the numerical derivatives and the
         !! user supplied derivatives.
      integer, intent(in) :: nrow
         !! Row number of the explanatory variable array at which the derivative is checked.
      logical, intent(in) :: isodr
         !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`).
      real(wp), intent(in) :: epsmac
         !! Value of machine precision.
      real(wp), intent(in) :: pv0i(n, q)
         !! Predicted values using the user supplied parameter estimates.
      real(wp), intent(out) :: fjacb(n, np, q)
         !! Jacobian with respect to `beta`.
      real(wp), intent(out) :: fjacd(n, m, q)
         !! Jacobian with respect to `delta`.
      integer, intent(out) :: msgb(1 + q*np)
         !! Error checking results for the Jacobian wrt `beta`.
      integer, intent(out) :: msgd(1 + q*m)
         !! Error checking results for the Jacobian wrt `delta`.
      real(wp), intent(out) :: diff(q, np + m)
         !! Relative differences between the user supplied and finite difference derivatives
         !! for each derivative checked.
      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.
      integer, intent(inout) :: njev
         !! Number of Jacobian 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) :: diffj, h0, hc0, pv, tol, typj
      integer :: ideval, j, lq, msgb1, msgd1
      logical :: isfixd, iswrtb

      ! Local arrays
      real(wp) :: pv0(n, q)

      ! Variable Definitions (alphabetically)
      !  DIFFJ:    The relative differences between the user supplied and finite difference
      !            derivatives for the derivative being checked.
      !  H0:       The initial relative step size for forward differences.
      !  HC0:      The initial relative step size for central differences.
      !  IDEVAL:   The variable designating what computations are to be performed by user supplied
      !            subroutine FCN.
      !  ISFIXD:   The variable designating whether the parameter is fixed (TRUE) or not (FALSE).
      !  ISWRTB:   The variable designating whether the derivatives wrt BETA (TRUE) or DELTA
      !            (FALSE) are being checked.
      !  J:        An index variable.
      !  LQ:       The response currently being examined.
      !  MSGB1:    The error checking results for the Jacobian wrt BETA.
      !  MSGD1:    The error checking results for the Jacobian wrt DELTA.
      !  PV:       The scalar in which the predicted value from the model for row NROW is stored.
      !  PV0:      The predicted values using the current parameter estimates (possibly offset from
      !            the user supplied estimates to create distance between parameters and the bounds
      !            on the parameters).
      !  TOL:      The agreement tolerance.
      !  TYPJ:     The typical size of the J-th unknown BETA or DELTA.

      ! Set tolerance for checking derivatives
      tol = eta**(0.25E0_wp)
      ntol = int(max(one, p5 - log10(tol)))

      ! Compute, if necessary, PV0
      pv0 = pv0i
      if (any(beta /= betaj)) then
         istop = 0
         ideval = 001
         call fcn(betaj, xplusd, ifixb, ifixx, ideval, pv0, fjacb, fjacd, istop)
         if (istop /= 0) then
            return
         else
            njev = njev + 1
         end if
      end if

      ! Compute user-supplied derivative values
      istop = 0
      if (isodr) then
         ideval = 110
      else
         ideval = 010
      end if
      call fcn(betaj, xplusd, ifixb, ifixx, ideval, wrk2, fjacb, fjacd, istop)
      if (istop /= 0) then
         return
      else
         njev = njev + 1
      end if

      ! Check derivatives wrt BETA for each response of observation NROW
      msgb1 = 0
      msgd1 = 0

      do lq = 1, q

         ! Set predicted value of model at current parameter estimates
         pv = pv0(nrow, lq)

         iswrtb = .true.
         do j = 1, np

            if (ifixb(1) < 0) then
               isfixd = .false.
            elseif (ifixb(j) == 0) then
               isfixd = .true.
            else
               isfixd = .false.
            end if

            if (isfixd) then
               msgb(1 + lq + (j - 1)*q) = -1
            else
               if (beta(j) == zero) then
                  if (ssf(1) < zero) then
                     typj = one/abs(ssf(1))
                  else
                     typj = one/ssf(j)
                  end if
               else
                  typj = abs(beta(j))
               end if

               h0 = hstep(0, neta, 1, j, stpb, 1)
               hc0 = h0

               ! Check derivative wrt the J-th parameter at the NROW-th row
               if (interval(j) >= 1) then
                  call check_jac_value(fcn, &
                                       n, m, np, q, &
                                       betaj, xplusd, &
                                       ifixb, ifixx, ldifx, &
                                       eta, tol, nrow, epsmac, j, lq, typj, h0, hc0, &
                                       iswrtb, pv, fjacb(nrow, j, lq), &
                                       diffj, msgb1, msgb(2), istop, nfev, &
                                       wrk1, wrk2, wrk6, interval)
                  if (istop /= 0) then
                     msgb(1) = -1
                     return
                  else
                     diff(lq, j) = diffj
                  end if
               else
                  msgb(1 + j) = 9
               end if
            end if

         end do

         ! Check derivatives wrt X for each response of observation NROW
         if (isodr) then

            iswrtb = .false.
            do j = 1, m

               if (ifixx(1, 1) < 0) then
                  isfixd = .false.
               elseif (ldifx == 1) then
                  if (ifixx(1, j) == 0) then
                     isfixd = .true.
                  else
                     isfixd = .false.
                  end if
               else
                  isfixd = .false.
               end if

               if (isfixd) then
                  msgd(1 + lq + (j - 1)*q) = -1
               else
                  if (xplusd(nrow, j) == zero) then
                     if (tt(1, 1) < zero) then
                        typj = one/abs(tt(1, 1))
                     elseif (ldtt == 1) then
                        typj = one/tt(1, j)
                     else
                        typj = one/tt(nrow, j)
                     end if
                  else
                     typj = abs(xplusd(nrow, j))
                  end if

                  h0 = hstep(0, neta, nrow, j, stpd, ldstpd)
                  hc0 = hstep(1, neta, nrow, j, stpd, ldstpd)

                  ! Check derivative wrt the J-th column of DELTA at row NROW
                  call check_jac_value(fcn, &
                                       n, m, np, q, &
                                       betaj, xplusd, &
                                       ifixb, ifixx, ldifx, &
                                       eta, tol, nrow, epsmac, j, lq, typj, h0, hc0, &
                                       iswrtb, pv, fjacd(nrow, j, lq), &
                                       diffj, msgd1, msgd(2), istop, nfev, &
                                       wrk1, wrk2, wrk6, interval)
                  if (istop /= 0) then
                     msgd(1) = -1
                     return
                  else
                     diff(lq, np + j) = diffj
                  end if
               end if

            end do

         end if

      end do

      msgb(1) = msgb1
      msgd(1) = msgd1

   end subroutine check_jac