jac_cdiff Subroutine

   subroutine jac_cdiff &
      (fcn, &
       n, m, np, q, &
       beta, x, delta, xplusd, ifixb, ifixx, ldifx, &
       stpb, stpd, ldstpd, &
       ssf, tt, ldtt, neta, fn, stp, wrk1, wrk2, wrk3, wrk6, &
       fjacb, isodr, fjacd, nfev, istop, info, &
       lower, upper)
   !! Compute central difference approximations to the Jacobian wrt the estimated `beta`s and
   !! wrt the `delta`s.

      use odrpack_kinds, only: zero, one

      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(in) :: x(n, m)
         !! Explanatory variable.
      real(wp), intent(in) :: delta(n, m)
         !! Estimated errors in the explanatory variables.
      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)
         !! Relative step used for computing finite difference derivatives with respect to each `beta`.
      real(wp), intent(in) :: stpd(ldstpd, m)
         !! Relative step used for computing finite difference derivatives with respect to each `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`.
      integer, intent(in) :: neta
         !! Number of good digits in the function results.
      real(wp), intent(in) :: fn(n, q)
         !! New predicted values from the function. Used when parameter is on a boundary.
      real(wp), intent(out) :: stp(n)
         !! Step used for computing finite difference derivatives with respect to each `delta`.
      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) :: wrk3(np)
         !! A work array of `(np)` elements.
      real(wp), intent(out) :: wrk6(n, np, q)
         !! A work array of `(n, np, q)` elements.
      real(wp), intent(out) :: fjacb(n, np, q)
         !! Jacobian with respect to `beta`.
      logical, intent(in) :: isodr
         !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`).
      real(wp), intent(out) :: fjacd(n, m, q)
         !! Jacobian with respect to `delta`.
      integer, intent(inout) :: nfev
         !! Number of function evaluations.
      integer, intent(out) :: istop
         !! Variable designating whether there are problems computing the function at the
         !! current `beta` and `delta`.
      integer, intent(out) :: info
         !! Variable designating why the computations were stopped.
      real(wp), intent(in) :: lower(np)
         !! Lower bound on `beta`.
      real(wp), intent(in) :: upper(np)
         !! Upper bound on `beta`.

      ! Local scalars
      real(wp) :: betak, typj
      integer :: i, j, k, l
      logical :: doit, setzero

      ! Variable Definitions (alphabetically)
      !  BETAK:    The K-th function parameter.
      !  DOIT:     The variable designating whether the derivative wrt a given BETA or DELTA
      !            needs to be computed (TRUE) or not (FALSE).
      !  I:        An indexing variable.
      !  J:        An indexing variable.
      !  K:        An indexing variable.
      !  L:        An indexing variable.
      !  SETZERO:  The variable designating whether the derivative wrt some DELTA needs to be
      !            set to zero (TRUE) or not (FALSE).
      !  TYPJ:     The typical size of the J-th unknown BETA or DELTA.

      ! Compute the Jacobian wrt the estimated BETAS
      do k = 1, np
         if (ifixb(1) >= 0) then
            if (ifixb(k) == 0) then
               doit = .false.
            else
               doit = .true.
            end if
         else
            doit = .true.
         end if
         if (.not. doit) then
            fjacb(:, k, :) = zero
         else
            betak = beta(k)
            wrk3(k) = betak + derstep(1, k, betak, ssf, stpb, neta)
            wrk3(k) = wrk3(k) - betak

            beta(k) = betak + wrk3(k)
            if (beta(k) > upper(k)) then
               beta(k) = upper(k)
            elseif (beta(k) < lower(k)) then
               beta(k) = lower(k)
            end if
            if (beta(k) - 2*wrk3(k) < lower(k)) then
               beta(k) = lower(k) + 2*wrk3(k)
            elseif (beta(k) - 2*wrk3(k) > upper(k)) then
               beta(k) = upper(k) + 2*wrk3(k)
            end if
            if (beta(k) > upper(k) .or. beta(k) < lower(k)) then
               info = 60001
               return
            end if

            istop = 0
            if (beta(k) == betak) then
               wrk2 = fn
            else
               call fcn(beta, xplusd, ifixb, ifixx, 001, wrk2, wrk6, wrk1, istop)
               if (istop /= 0) then
                  return
               else
                  nfev = nfev + 1
               end if
            end if
            fjacb(:, k, :) = wrk2

            beta(k) = beta(k) - 2*wrk3(k)
            if (beta(k) > upper(k)) then
               info = 60001
               return
            end if
            if (beta(k) < lower(k)) then
               info = 60001
               return
            end if
            istop = 0
            if (beta(k) == betak) then
               wrk2 = fn
            else
               call fcn(beta, xplusd, ifixb, ifixx, 001, wrk2, wrk6, wrk1, istop)
               if (istop /= 0) then
                  return
               else
                  nfev = nfev + 1
               end if
            end if

            fjacb(:, k, :) = (fjacb(:, k, :) - wrk2)/(2*wrk3(k))
            beta(k) = betak

         end if
      end do

      ! Compute the Jacobian wrt the X's
      if (isodr) then
         do j = 1, m
            if (ifixx(1, 1) < 0) then
               doit = .true.
               setzero = .false.
            elseif (ldifx == 1) then
               if (ifixx(1, j) == 0) then
                  doit = .false.
               else
                  doit = .true.
               end if
               setzero = .false.
            else
               doit = any(ifixx(:, j) /= 0)
               setzero = any(ifixx(:, j) == 0)
            end if

            if (.not. doit) then
               fjacd(:, j, :) = zero
            else
               do i = 1, n

                  if (xplusd(i, j) == zero) then
                     if (tt(1, 1) < zero) then
                        typj = 1/abs(tt(1, 1))
                     elseif (ldtt == 1) then
                        typj = 1/tt(1, j)
                     else
                        typj = 1/tt(i, j)
                     end if
                  else
                     typj = abs(xplusd(i, j))
                  end if

                  stp(i) = xplusd(i, j) &
                           + sign(one, xplusd(i, j))*typj*hstep(1, neta, i, j, stpd, ldstpd)
                  stp(i) = stp(i) - xplusd(i, j)
                  xplusd(i, j) = xplusd(i, j) + stp(i)

               end do

               istop = 0
               call fcn(beta, xplusd, ifixb, ifixx, 001, wrk2, wrk6, wrk1, istop)
               if (istop /= 0) then
                  return
               else
                  nfev = nfev + 1
                  fjacd(:, j, :) = wrk2
               end if

               xplusd(:, j) = x(:, j) + delta(:, j) - stp

               istop = 0
               call fcn(beta, xplusd, ifixb, ifixx, 001, wrk2, wrk6, wrk1, istop)
               if (istop /= 0) then
                  return
               else
                  nfev = nfev + 1
               end if

               if (setzero) then
                  do i = 1, n
                     if (ifixx(i, j) == 0) then
                        fjacd(i, j, :) = zero
                     else
                        fjacd(i, j, :) = (fjacd(i, j, :) - wrk2(i, :))/(2*stp(i))
                     end if
                  end do
               else
                  do l = 1, q
                     fjacd(:, j, l) = (fjacd(:, j, l) - wrk2(:, l))/(2*stp(:))
                  end do
               end if

               xplusd(:, j) = x(:, j) + delta(:, j)

            end if

         end do
      end if

   end subroutine jac_cdiff