djacfd Subroutine

   subroutine djacfd &
      (fcn, &
       n, m, np, nq, &
       beta, x, ldx, 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 forward difference approximations to the Jacobian wrt the estimated `beta`s and
   !! wrt the `delta`s.
      ! Routines Called  FCN, DHSTEP, DERSTEP
      ! Date Written   860529   (YYMMDD)
      ! Revision Date  920619   (YYMMDD)

      use odrpack_kinds, only: zero, one

      procedure(fcn_t) :: fcn
         !! The user supplied subroutine for evaluating the model.
      integer, intent(in) :: n
         !! The number of observations.
      integer, intent(in) :: m
         !! The number of columns of data in the explanatory variable.
      integer, intent(in) :: np
         !! The number of function parameters.
      integer, intent(in) :: nq
         !! The number of responses per observation.
      real(wp), intent(inout) :: beta(np)
         !! The function parameters.
      real(wp), intent(in) :: x(ldx, m)
         !! The explanatory variable.
      integer, intent(in) :: ldx
         !! The leading dimension of array `x`.
      real(wp), intent(in) :: delta(n, m)
         !! The estimated errors in the explanatory variables.
      real(wp), intent(inout) :: xplusd(n, m)
         !! The values of `x + delta`.
      integer, intent(in) :: ifixb(np)
         !! The values designating whether the elements of `beta` are fixed at their input values or not.
      integer, intent(in) :: ifixx(ldifx, m)
         !! The values designating whether the elements of `x` are fixed at their input values or not.
      integer, intent(in) :: ldifx
         !! The leading dimension of array `ifixx`.
      real(wp), intent(in) :: stpb(np)
         !! The relative step used for computing finite difference derivatives with respect to each `beta`.
      real(wp), intent(in) :: stpd(ldstpd, m)
         !! The relative step used for computing finite difference derivatives with respect to each `delta`.
      integer, intent(in) :: ldstpd
         !! The leading dimension of array `stpd`.
      real(wp), intent(in) :: ssf(np)
         !! The scaling values used for `beta`.
      real(wp), intent(in) :: tt(ldtt, m)
         !! The scaling values used for `delta`.
      integer, intent(in) :: ldtt
         !! The leading dimension of array `tt`.
      integer, intent(in) :: neta
         !! The number of good digits in the function results.
      real(wp), intent(in) :: fn(n, nq)
         !! The new predicted values from the function. Used when parameter is on a boundary.
      real(wp), intent(out) :: stp(n)
         !! The step used for computing finite difference derivatives with respect to each `delta`.
      real(wp), intent(out) :: wrk1(n, m, nq)
         !! A work array of `(n, m, nq)` elements.
      real(wp), intent(out) :: wrk2(n, nq)
         !! A work array of `(n, nq)` elements.
      real(wp), intent(out) :: wrk3(np)
         !! A work array of `(np)` elements.
      real(wp), intent(out) :: wrk6(n, np, nq)
         !! A work array of `(n, np, nq)` elements.
      real(wp), intent(out) :: fjacb(n, np, nq)
         !! The Jacobian with respect to `beta`.
      logical, intent(in) :: isodr
         !! The variable designating whether the solution is by ODR (`isodr = .true.`) or
         !! by OLS (`isodr = .false.`).
      real(wp), intent(out) :: fjacd(n, m, nq)
         !! The Jacobian with respect to `delta`.
      integer, intent(inout) :: nfev
         !! The number of function evaluations.
      integer, intent(out) :: istop
         !! The variable designating whether there are problems computing the function at the
         !! current `beta` and `delta`.
      integer, intent(out) :: info
         !! The variable designating why the computations were stopped.
      real(wp), intent(in) :: lower(np)
         !! The lower bound on `beta`.
      real(wp), intent(in) :: upper(np)
         !! The upper bound on `beta`.

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

      ! Variable Definitions (alphabetically)
      !  BETA:    The function parameters.
      !  BETAK:   The K-th function parameter.
      !  DELTA:   The estimated errors in the explanatory variables.
      !  DOIT:    The variable designating whether the derivative wrt a given BETA or DELTA needs
      !           to be computed (DOIT=TRUE) or not (DOIT=FALSE).
      !  FCN:     The user supplied subroutine for evaluating the model.
      !  FJACB:   The Jacobian with respect to BETA.
      !  FJACD:   The Jacobian with respect to DELTA.
      !  FN:      The new predicted values from the function.
      !  I:       An indexing variable.
      !  IFIXB:   The values designating whether the elements of BETA are fixed at their input
      !           values or not.
      !  IFIXX:   The values designating whether the elements of X are fixed at their input values
      !           or not.
      !  ISODR:   The variable designating whether the solution is by ODR (ISODR=TRUE) or by
      !           OLS (ISODR=FALSE).
      !  ISTOP:   The variable designating whether there are problems computing the function at
      !           the current BETA and DELTA.
      !  J:       An indexing variable.
      !  K:       An indexing variable.
      !  L:       An indexing variable.
      !  LDIFX:   The leading dimension of array IFIXX.
      !  LDSTPD:  The leading dimension of array STPD.
      !  LDTT:    The leading dimension of array TT.
      !  LDX:     The leading dimension of array X.
      !  M:       The number of columns of data in the explanatory variable.
      !  N:       The number of observations.
      !  NETA:    The number of good digits in the function results.
      !  NFEV:    The number of function evaluations.
      !  NP:      The number of function parameters.
      !  SETZRO:  The variable designating whether the derivative wrt some DELTA needs to be set
      !           to zero (SETZRO=TRUE) or not (SETZRO=FALSE).
      !  SSF:     The scale used for the BETA'S.
      !  STP:     The step used for computing finite difference derivatives with respect to DELTA.
      !  STPB:    The relative step used for computing finite difference derivatives with respect
      !           to BETA.
      !  STPD:    The relative step used for computing finite difference derivatives with respect
      !           to DELTA.
      !  TT:      The scaling values used for DELTA.
      !  TYPJ:    The typical size of the J-th unknown BETA or DELTA.
      !  X:       The explanatory variable.
      !  XPLUSD:  The values of X + DELTA.
      !  WRK1:    A work array of (N by M by NQ) elements.
      !  WRK2:    A work array of (N BY NQ) elements.
      !  WRK3:    A work array of (NP) elements.
      !  WRK6:    A work array of (N BY NP BY NQ) elements.

      ! 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
            do l = 1, nq
               fjacb(1:n, k, l) = zero
            end do
         else
            betak = beta(k)
            step = derstep(0, k, betak, ssf, stpb, neta)
            wrk3(k) = betak + step
            wrk3(k) = wrk3(k) - betak
            beta(k) = betak + wrk3(k)
            if (beta(k) > upper(k)) then
               step = -step
               wrk3(k) = betak + step
               wrk3(k) = wrk3(k) - betak
               beta(k) = betak + wrk3(k)
            end if
            if (beta(k) < lower(k)) then
               step = -step
               wrk3(k) = betak + step
               wrk3(k) = wrk3(k) - betak
               beta(k) = betak + wrk3(k)
               if (beta(k) > upper(k)) then
                  info = 60001
                  return
               end if
            end if
            istop = 0
            call fcn(n, m, np, nq, &
                     n, m, np, &
                     beta, xplusd, &
                     ifixb, ifixx, ldifx, &
                     001, wrk2, wrk6, wrk1, &
                     istop)
            if (istop /= 0) then
               return
            else
               nfev = nfev + 1
            end if
            do l = 1, nq
               do i = 1, n
                  fjacb(i, k, l) = (wrk2(i, l) - fn(i, l))/wrk3(k)
               end do
            end do
            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.
               setzro = .false.
            elseif (ldifx == 1) then
               if (ifixx(1, j) == 0) then
                  doit = .false.
               else
                  doit = .true.
               end if
               setzro = .false.
            else
               doit = .false.
               setzro = .false.
               do i = 1, n
                  if (ifixx(i, j) /= 0) then
                     doit = .true.
                  else
                     setzro = .true.
                  end if
               end do
            end if
            if (.not. doit) then
               do l = 1, nq
                  fjacd(1:n, j, l) = zero
               end do
            else
               do i = 1, n
                  if (xplusd(i, 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(i, j)
                     end if
                  else
                     typj = abs(xplusd(i, j))
                  end if

                  stp(i) = xplusd(i, j) &
                           + sign(one, xplusd(i, j)) &
                           *typj*dhstep(0, 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(n, m, np, nq, &
                        n, m, np, &
                        beta, xplusd, &
                        ifixb, ifixx, ldifx, &
                        001, wrk2, wrk6, wrk1, &
                        istop)
               if (istop /= 0) then
                  return
               else
                  nfev = nfev + 1
                  do l = 1, nq
                     do i = 1, n
                        fjacd(i, j, l) = wrk2(i, l)
                     end do
                  end do

               end if

               if (setzro) then
                  do i = 1, n
                     if (ifixx(i, j) == 0) then
                        do l = 1, nq
                           fjacd(i, j, l) = zero
                        end do
                     else
                        do l = 1, nq
                           fjacd(i, j, l) = (fjacd(i, j, l) - fn(i, l))/ &
                                            stp(i)
                        end do
                     end if
                  end do
               else
                  do l = 1, nq
                     do i = 1, n
                        fjacd(i, j, l) = (fjacd(i, j, l) - fn(i, l))/stp(i)
                     end do
                  end do
               end if
               do i = 1, n
                  xplusd(i, j) = x(i, j) + delta(i, j)
               end do
            end if
         end do
      end if

   end subroutine djacfd