eval_jac Subroutine

   subroutine eval_jac &
      (fcn, &
       anajac, cdjac, &
       n, m, np, q, &
       betac, beta, stpb, &
       ifixb, ifixx, ldifx, &
       x, delta, xplusd, stpd, ldstpd, &
       ssf, tt, ldtt, neta, fn, &
       stp, wrk1, wrk2, wrk3, wrk6, tempret, &
       fjacb, isodr, fjacd, we1, ldwe, ld2we, &
       njev, nfev, istop, info, &
       lower, upper)
   !! Compute the weighted Jacobians wrt `beta` and `delta`.

      use odrpack_kinds, only: zero

      procedure(fcn_t) :: fcn
         !! User-supplied subroutine for evaluating the model.
      logical, intent(in) :: anajac
         !! Variable designating whether the Jacobians are computed by finite differences
         !! (`.false.`) or not (`.true.`).
      logical, intent(in) :: cdjac
         !! Variable designating whether the Jacobians are computed by central differences
         !! (`.true.`) or by forward differences (`.false.`).
      integer, intent(in) :: n
         !! Number of observations.
      integer, intent(in) :: m
         !! Number of columns of data in the independent variable.
      integer, intent(in) :: np
         !! Number of function parameters.
      integer, intent(in) :: q
         !! Number of responses per observation.
      real(wp), intent(in) :: betac(np)
         !! Current estimated values of the unfixed `beta`s.
      real(wp), intent(out) :: beta(np)
         !! Function parameters.
      real(wp), intent(in) :: stpb(np)
         !! Relative step used for computing finite difference derivatives with respect to `beta`.
      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 `delta` are fixed at their input values or not.
      integer, intent(in) :: ldifx
         !! Leading dimension of array `ifixx`.
      real(wp), intent(in) :: x(n, m)
         !! Independent variable.
      real(wp), intent(in) :: delta(n, m)
         !! Estimated values of `delta`.
      real(wp), intent(out) :: xplusd(n, m)
         !! Values of `x + delta`.
      real(wp), intent(in) :: stpd(ldstpd, m)
         !! Relative step used for computing finite difference derivatives with respect to `delta`.
      integer, intent(in) :: ldstpd
         !! Leading dimension of array `stpd`.
      real(wp), intent(in) :: ssf(np)
         !! Scale used for the `beta`s.
      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 accurate digits in the function results.
      real(wp), intent(in) :: fn(n, q)
         !! Predicted values of the function at the current point.
      real(wp), intent(out) :: stp(n)
         !! Step used for computing finite difference derivatives with respect to `delta`.
      real(wp), intent(out) :: wrk1(n, m, q)
         !! Work array of `(n, m, q)` elements.
      real(wp), intent(out) :: wrk2(n, q)
         !! Work array of `(n, q)` elements.
      real(wp), intent(out) :: wrk3(np)
         !! Work array of `(np)` elements.
      real(wp), intent(out) :: wrk6(n, np, q)
         !! Work array of `(n, np, q)` elements.
      real(wp), intent(inout) :: tempret(:, :)
         !! Temporary work array for holding return values before copying to a lower rank array.
      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`.
      real(wp), intent(in) :: we1(ldwe, ld2we, q)
         !! Square roots of the `epsilon` weights in array `we`.
      integer, intent(in) :: ldwe
         !! Leading dimension of arrays `we` and `we1`.
      integer, intent(in) :: ld2we
         !! Second dimension of arrays `we` and `we1`.
      integer, intent(inout) :: njev
         !! Number of Jacobian evaluations.
      integer, intent(inout) :: nfev
         !! Number of function evaluations.
      integer, intent(out) :: istop
         !! Variable designating that the user wishes the computations stopped.
      integer, intent(out) :: info
         !! Variable designating why the computations were stopped.
      real(wp), intent(in) :: lower(np)
         !! Lower bound of `beta`.
      real(wp), intent(in) :: upper(np)
         !! Upper bound of `beta`.

      ! Local scalars
      integer :: ideval, j, j1
      logical :: ferror

      ! Variable Definitions (alphabetically)
      !  FERROR:  The variable designating whether ODRPACK95 detected nonzero values in array DELTA
      !           in the OLS case, and thus whether the user may have overwritten important information
      !           by computing FJACD in the OLS case.
      !  IDEVAL:  The variable designating what computations are to be performed by user-supplied
      !           subroutine FCN.
      !  J:       An indexing variable.
      !  J1:      An indexing variable.

      ! Insert current unfixed BETA estimates into BETA
      call unpack_vec(np, betac, beta, ifixb)

      ! Compute XPLUSD = X + DELTA
      xplusd = x + delta

      ! Compute the Jacobian wrt the estimated BETAS (FJACB) and the Jacobian wrt DELTA (FJACD)
      istop = 0
      if (isodr) then
         ideval = 110
      else
         ideval = 010
      end if
      if (anajac) then
         call fcn(beta, xplusd, ifixb, ifixx, ideval, wrk2, fjacb, fjacd, istop)
         if (istop /= 0) then
            return
         else
            njev = njev + 1
         end if
         ! Make sure fixed elements of FJACD are zero
         if (isodr) then
            do j = 1, q
               call set_ifix(n, m, ifixx, ldifx, fjacd(1, 1, j), n, fjacd(1, 1, j), n)
            end do
         end if
      elseif (cdjac) then
         call 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)
      else
         call jac_fwdiff(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)
      end if
      if (istop < 0 .or. info >= 10000) then
         return
      elseif (.not. isodr) then
         ! Try to detect whether the user has computed JFACD within FCN in the OLS case
         ferror = sum(delta**2) /= zero
         if (ferror) then
            info = 50300
            return
         end if
      end if

      ! Weight the Jacobian wrt the estimated BETAS
      if (ifixb(1) < 0) then
         do j = 1, np
            call scale_mat(n, q, we1, ldwe, ld2we, &
                           reshape(fjacb(:, j, :), [n, q]), &
                           tempret(1:n, 1:q))
            fjacb(:, j, :) = tempret(1:n, 1:q)
         end do
      else
         j1 = 0
         do j = 1, np
            if (ifixb(j) >= 1) then
               j1 = j1 + 1
               call scale_mat(n, q, we1, ldwe, ld2we, &
                              reshape(fjacb(:, j, :), [n, q]), &
                              tempret(1:n, 1:q))
               fjacb(:, j1, :) = tempret(1:n, 1:q)
            end if
         end do
      end if

      ! Weight the Jacobian's wrt DELTA as appropriate
      if (isodr) then
         do j = 1, m
            call scale_mat(n, q, we1, ldwe, ld2we, &
                           reshape(fjacd(:, j, :), [n, q]), &
                           tempret(1:n, 1:q))
            fjacd(:, j, :) = tempret(1:n, 1:q)
         end do
      end if

   end subroutine eval_jac