doddrv Subroutine

   impure subroutine doddrv &
      (head, fstitr, prtpen, &
       fcn, n, m, np, nq, beta, y, ldy, x, ldx, &
       we, ldwe, ld2we, wd, ldwd, ld2wd, ifixb, ifixx, ldifx, &
       job, ndigit, taufac, sstol, partol, maxit, &
       iprint, lunerr, lunrpt, &
       stpb, stpd, ldstpd, sclb, scld, ldscld, &
       work, lwork, tempret, iwork, liwork, &
       maxit1, tstimp, info, lower, upper)
   !! Perform error checking and initialization, and begin procedure for performing orthogonal
   !! distance regression (ODR) or ordinary linear or nonlinear least squares (OLS).
      ! Routines Called  FCN, DCOPY, DDOT, DETAF, DFCTRW, DFLAGS,
      !                  DINIWK, DIWINF, DJCK, DNRM2, DODCHK, DODMN,
      !                  DODPER, DPACK, DSETN, DUNPAC, DWGHT, DWINF, DERSTEP
      ! Date Written   860529   (YYMMDD)
      ! Revision Date  920619   (YYMMDD)

      use odrpack_kinds, only: zero, one, ten, p5 => half
      use odrpack_core, only: fcn_t, detaf, dfctrw, dflags, diniwk, diwinf, djck, dodchk, &
                              dpack, dsetn, dunpac, dwght, dwinf, derstep, mbfb
      use odrpack_reports, only: dodper
      use blas_interfaces, only: ddot, dnrm2, dcopy

      logical, intent(inout) :: head
         !! The variable designating whether the heading is to be printed (`head = .true.`)
         !! or not (`head = .false.`).
      logical, intent(inout) :: fstitr
         !! The variable designating whether this is the first iteration (`fstitr = .true.`)
         !! or not (`fstitr = .false.`).
      logical, intent(inout) :: prtpen
         !! The variable designating whether the penalty parameter is to be printed in the
         !! iteration report (`prtpen = .true.`) or not (`prtpen = .false.`).
      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) :: y(ldy, nq)
         !! The dependent variable. Unused when the model is implicit.
      integer, intent(in) :: ldy
         !! The leading dimension of array `y`.
      real(wp), intent(in) :: x(ldx, m)
         !! The explanatory variable.
      integer, intent(in) :: ldx
         !! The leading dimension of array `x`.
      real(wp), intent(inout) :: we(ldwe, ld2we, nq)
         !! The `epsilon` weights.
      integer, intent(in) :: ldwe
         !! The leading dimension of array `we`.
      integer, intent(in) :: ld2we
         !! The second dimension of array `we`.
      real(wp), intent(in) :: wd(ldwd, ld2wd, m)
         !! The `delta` weights.
      integer, intent(in) :: ldwd
         !! The leading dimension of array `wd`.
      integer, intent(in) :: ld2wd
         !! The second dimension of array `wd`.
      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`.
      integer, intent(inout) :: job
         !! The variable controlling problem initialization and computational method.
      integer, intent(in) :: ndigit
         !! The number of accurate digits in the function results, as supplied by the user.
      real(wp), intent(in) :: taufac
         !! The factor used to compute the initial trust region diameter.
      real(wp), intent(in) :: sstol
         !! The sum-of-squares convergence stopping tolerance.
      real(wp), intent(in) :: partol
         !! The parameter convergence stopping tolerance.
      integer, intent(in) :: maxit
         !! The maximum number of iterations allowed.
      integer, intent(in) :: iprint
         !! The print control variable.
      integer, intent(in) :: lunerr
         !! The logical unit number used for error messages.
      integer, intent(in) :: lunrpt
         !! The logical unit number used for computation reports.
      real(wp), intent(in) :: stpb(np)
         !! The step size for finite difference derivatives with respect to `beta`.
      real(wp), intent(in) :: stpd(ldstpd, m)
         !! The step size for finite difference derivatives with respect to `delta`.
      integer, intent(in) :: ldstpd
         !! The leading dimension of array `stpd`.
      real(wp), intent(in) :: sclb(np)
         !! The scaling values for `beta`.
      real(wp), intent(in) :: scld(ldscld, m)
         !! The scaling values for `delta`.
      integer, intent(in) :: ldscld
         !! The leading dimension of array `scld`.
      real(wp), intent(inout) :: work(lwork)
         !! The real work space.
      integer, intent(in) :: lwork
         !! The length of vector `work`.
      real(wp), intent(inout) :: tempret(:, :)
         !! Temporary work array for holding return values before copying to a lower rank array.
      integer, intent(inout) :: iwork(liwork)
         !! The integer work space.
      integer, intent(in) :: liwork
         !! The length of vector `iwork`.
      integer, intent(out) :: maxit1
         !! For implicit models, the iterations allowed for the next penalty parameter value.
      real(wp), intent(out) :: tstimp
         !! The relative change in the parameters between the initial values and the solution.
      integer, intent(inout) :: info
         !! The variable designating why the computations were stopped.
      real(wp), intent(in) :: lower(np)
         !! The lower bound for `beta`.
      real(wp), intent(in) :: upper(np)
         !! The upper bound for `beta`.

      ! Local scalars
      real(wp) :: epsmac, eta
      integer :: actrsi, alphai, betaci, betani, betasi, beta0i, boundi, deltai, &
                 deltni, deltsi, diffi, epsmai, etai, fi, fjacbi, fjacdi, fni, fsi, i, &
                 idfi, int2i, iprini, iranki, istop, istopi, jobi, jpvti, k, ldtt, &
                 ldtti, liwkmn, loweri, luneri, lunrpi, lwkmn, lwrk, maxiti, msgb, &
                 msgd, neta, netai, nfev, nfevi, niteri, njev, njevi, nnzw, nnzwi, &
                 npp, nppi, nrow, nrowi, ntol, ntoli, olmavi, omegai, partli, pnormi, &
                 prersi, qrauxi, rcondi, rnorsi, rvari, sdi, si, ssfi, ssi, sstoli, &
                 taufci, taui, ti, tti, ui, upperi, vcvi, we1i, wrk1i, wrk2i, wrk3i, &
                 wrk4i, wrk5i, wrk6i, wrk7i, wrk, wssi, wssdei, wssepi, xplusi
      logical :: anajac, cdjac, chkjac, dovcv, implct, initd, isodr, redoj, restrt

      ! Local arrays
      real(wp) :: betaj(np)
      integer :: interval(np)

      ! Variable Definitions (alphabetically)
      !  ACTRSI:   The location in array work of variable ACTRS.
      !  ALPHAI:   The location in array work of variable ALPHA.
      !  ANAJAC:   The variable designating whether the Jacobians are computed by finite
      !            differences (ANAJAC=FALSE) or not (ANAJAC=TRUE).
      !  BETA:     The function parameters.
      !  BETACI:   The starting location in array WORK of array BETAC.
      !  BETAJ:    The parameters to use in the derivative checking algorithm.
      !  BETANI:   The starting location in array WORK of array BETAN.
      !  BETASI:   The starting location in array WORK of array BETAS.
      !  BETA0I:   The starting location in array WORK of array BETA0.
      !  CDJAC:    The variable designating whether the Jacobians are computed by central
      !            differences (CDJAC=TRUE) or forward differences (CDJAC=FALSE).
      !  CHKJAC:   The variable designating whether the user supplied Jacobians are to be
      !            checked (CHKJAC=TRUE) or not (CHKJAC=FALSE).
      !  DELTAI:   The starting location in array WORK of array DELTA.
      !  DELTNI:   The starting location in array WORK of array DELTAN.
      !  DELTSI:   The starting location in array WORK of array DELTAS.
      !  DIFFI:    The starting location in array WORK of array DIFF.
      !  DOVCV:    The variable designating whether the covariance matrix is to be computed
      !            (DOVCV=TRUE) or not (DOVCV=FALSE).
      !  EPSMAI:   The location in array WORK of variable EPSMAC.
      !  ETA:      The relative noise in the function results.
      !  ETAI:     The location in array WORK of variable ETA.
      !  FCN:      THE USER SUPPLIED SUBROUTINE FOR EVALUATING THE MODEL.
      !  FI:       The starting location in array WORK of array F.
      !  FJACBI:   The starting location in array WORK of array FJACB.
      !  FJACDI:   The starting location in array WORK of array FJACD.
      !  FNI:      The starting location in array WORK of array FN.
      !  FSI:      The starting location in array WORK of array FS.
      !  FSTITR:   The variable designating whether this is the first iteration (FSTITR=TRUE)
      !            or not (FSTITR=FALSE).
      !  HEAD:     The variable designating whether the heading is to be printed (HEAD=TRUE)
      !            or not (HEAD=FALSE).
      !  I:        An index variable.
      !  IDFI:     The location in array iwork of variable IDF.
      !  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.
      !  IMPLCT:   The variable designating whether the solution is by implicit ODR (IMPLCT=TRUE)
      !            or explicit ODR (IMPLCT=FALSE).
      !  INFO:     The variable designating why the computations were stopped.
      !  INITD:    The variable designating whether DELTA is to be initialized to zero (INITD=TRUE)
      !            or to the values in the first N by M elements of array WORK (INITD=FALSE).
      !  INT2I:    The location in array IWORK of variable INT2.
      !  INTERVAL: Specifies which checks can be performed when checking derivatives based on the
      !            interval of the bound constraints.
      !  IPRINI:   The location in array iwork of variable IPRINT.
      !  IPRINT:   The print control variable.
      !  IRANKI:   The location in array IWORK of variable IRANK.
      !  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.
      !  ISTOPI:   The location in array IWORK of variable ISTOP.
      !  IWORK:    The integer work space.
      !  JOB:      The variable controling problem initialization and computational method.
      !  JOBI:     The location in array IWORK of variable JOB.
      !  JPVTI:    The starting location in array IWORK of array JPVT.
      !  K:        An index variable.
      !  LDIFX:    The leading dimension of array IFIXX.
      !  LDSCLD:   The leading dimension of array SCLD.
      !  LDSTPD:   The leading dimension of array STPD.
      !  LDTT:     The leading dimension of array TT.
      !  LDTTI:    The location in array IWORK of variable LDTT.
      !  LDWD:     The leading dimension of array WD.
      !  LDWE:     The leading dimension of array WE.
      !  LDX:      The leading dimension of array X.
      !  LDY:      The leading dimension of array Y.
      !  LD2WD:    The second dimension of array WD.
      !  LD2WE:    The second dimension of array WE.
      !  LIWKMN:   The minimum acceptable length of array IWORK.
      !  LIWORK:   The length of vector IWORK.
      !  LOWER:    The lower bound for BETA.
      !  LUNERI:   The location in array IWORK of variable LUNERR.
      !  LUNERR:   The logical unit number used for error messages.
      !  LUNRPI:   The location in array IWORK of variable LUNRPT.
      !  LUNRPT:   The logical unit number used for computation reports.
      !  LWKMN:    The minimum acceptable length of array WORK.
      !  LWORK:    The length of vector WORK.
      !  LWRK:     The length of vector WRK.
      !  M:        The number of columns of data in the explanatory variable.
      !  MAXIT:    The maximum number of iterations allowed.
      !  MAXIT1:   For implicit models, the iterations allowed for the next penalty parameter value.
      !  MAXITI:   The location in array IWORK of variable MAXIT.
      !  MSGB:     The starting location in array IWORK of array MSGB.
      !  MSGD:     The starting location in ARRAY IWORK of array MSGD.
      !  N:        The number of observations.
      !  NDIGIT:   The number of accurate digits in the function results, as supplied by the user.
      !  NETA:     The number of accurate digits in the function results.
      !  NETAI:    The location in array IWORK of variable NETA.
      !  NFEV:     The number of function evaluations.
      !  NFEVI:    The location in array IWORK of variable NFEV.
      !  NITERI:   The location in array IWORK of variable NITER.
      !  NJEV:     The number of Jacobian evaluations.
      !  NJEVI:    The location in array IWORK of variable NJEV.
      !  NNZW:     The number of nonzero observational error weights.
      !  NNZWI:    The location in array IWORK of variable NNZW.
      !  NP:       The number of function parameters.
      !  NPP:      The number of function parameters being estimated.
      !  NPPI:     The location in array IWORK of variable NPP.
      !  NQ:       The number of responses per observation.
      !  NROW:     The row number at which the derivative is to be checked.
      !  NROWI:    The location in array IWORK of variable NROW.
      !  NTOL:     The number of digits of agreement required between the numerical derivatives
      !            and the user supplied derivatives, set by DJCK.
      !  NTOLI:    The location in array IWORK of variable NTOL.
      !  OLMAVI:   The location in array WORK of variable OLMAVG.
      !  OMEGAI:   The starting location in array WORK of array OMEGA.
      !  PARTLI:   The location in array WORK of variable PARTOL.
      !  PARTOL:   The parameter convergence stopping tolerance.
      !  PNORM:    The norm of the scaled estimated parameters.
      !  PNORMI:   The location in array WORK of variable PNORM.
      !  PRERSI:   The location in array WORK of variable PRERS.
      !  PRTPEN:   The variable designating whether the penalty parameter is to be printed in
      !            the iteration report (PRTPEN=TRUE) or not (PRTPEN=FALSE).
      !  QRAUXI:   The starting location in array WORK of array QRAUX.
      !  RCONDI:   The location in array WORK of variable RCOND.
      !  REDOJ:    The variable designating whether the Jacobian matrix is to be recomputed for
      !            the computation of the covariance matrix (REDOJ=TRUE) or not (REDOJ=FALSE).
      !  RESTRT:   The variable designating whether the call is a restart (RESTRT=TRUE) or
      !            not (RESTRT=FALSE).
      !  RNORSI:   The location in array WORK of variable RNORMS.
      !  RVARI:    The location in array WORK of variable RVAR.
      !  SCLB:     The scaling values for BETA.
      !  SCLD:     The scaling values for DELTA.
      !  SDI:      The starting location in array WORK of array SD.
      !  SI:       The starting location in array WORK of array S.
      !  SSFI:     The starting location in array WORK of array SSF.
      !  SSI:      The starting location in array WORK of array SS.
      !  SSTOL:    The sum-of-squares convergence stopping tolerance.
      !  SSTOLI:   The location in array WORK of variable SSTOL.
      !  STPB:     The step size for finite difference derivatives wrt BETA.
      !  STPD:     The step size for finite difference derivatives wrt DELTA.
      !  TAUFAC:   The factor used to compute the initial trust region diameter.
      !  TAUFCI:   The location in array WORK of variable TAUFAC.
      !  TAUI:     The location in array WORK of variable TAU.
      !  TI:       The starting location in array WORK of array T.
      !  TSTIMP:   The relative change in the parameters between the initial values and
      !            the solution.
      !  TTI:      The starting location in array WORK of array TT.
      !  UI:       The starting location in array WORK of array U.
      !  UPPER:    The upper bound for BETA.
      !  VCVI:     The starting location in array WORK of array VCV.
      !  WD:       The DELTA weights.
      !  WE:       The EPSILON weights.
      !  WE1I:     The starting location in array WORK of array WE1.
      !  WORK:     The REAL (wp) work space.
      !  WRK:      The starting location in array WORK of array WRK, equivalenced to WRK1 and WRK2.
      !  WRK1I:    The starting location in array WORK of array WRK1.
      !  WRK2I:    The starting location in array WORK of array WRK2.
      !  WRK3I:    The starting location in array WORK of array WRK3.
      !  WRK4I:    The starting location in array WORK of array WRK4.
      !  WRK5I:    The starting location in array WORK of array WRK5.
      !  WRK6I:    The starting location in array WORK of array WRK6.
      !  WRK7I:    The starting location in array WORK of array WRK7.
      !  WSSI:     The location in array WORK of variable wss.
      !  WSSDEI:   The location in array WORK of variable WSSDEL.
      !  WSSEPI:   The location in array WORK of variable WSSEPS.
      !  X:        The explanatory variable.
      !  XPLUSI:   The starting location in array WORK of array XPLUSD.
      !  Y:        The dependent variable.  Unused when the model is implicit.

      ! Initialize necessary variables
      call dflags(job, restrt, initd, dovcv, redoj, anajac, cdjac, chkjac, isodr, implct)

      ! Set starting locations within integer workspace
      ! (invalid values of M, NP and/or NQ are handled reasonably by DIWINF)
      call diwinf(m, np, nq, &
                  msgb, msgd, jpvti, istopi, &
                  nnzwi, nppi, idfi, &
                  jobi, iprini, luneri, lunrpi, &
                  nrowi, ntoli, netai, &
                  maxiti, niteri, nfevi, njevi, int2i, iranki, ldtti, &
                  boundi, &
                  liwkmn)

      ! Set starting locations within REAL (wp) work space
      ! (invalid values of N, M, NP, NQ, LDWE and/or LD2WE
      ! are handled reasonably by DWINF)
      call dwinf(n, m, np, nq, ldwe, ld2we, isodr, &
                 deltai, fi, xplusi, fni, sdi, vcvi, &
                 rvari, wssi, wssdei, wssepi, rcondi, etai, &
                 olmavi, taui, alphai, actrsi, pnormi, rnorsi, prersi, &
                 partli, sstoli, taufci, epsmai, &
                 beta0i, betaci, betasi, betani, si, ssi, ssfi, qrauxi, ui, &
                 fsi, fjacbi, we1i, diffi, &
                 deltsi, deltni, ti, tti, omegai, fjacdi, &
                 wrk1i, wrk2i, wrk3i, wrk4i, wrk5i, wrk6i, wrk7i, &
                 loweri, upperi, &
                 lwkmn)

      if (isodr) then
         wrk = wrk1i
         lwrk = n*m*nq + n*nq
      else
         wrk = wrk2i
         lwrk = n*nq
      end if

      ! Update the penalty parameters
      ! (WE(1,1,1) is not a user supplied array in this case)
      if (restrt .and. implct) then
         we(1, 1, 1) = max(work(we1i)**2, abs(we(1, 1, 1)))
         work(we1i) = -sqrt(abs(we(1, 1, 1)))
      end if

      if (restrt) then

         ! Reset maximum number of iterations
         if (maxit >= 0) then
            iwork(maxiti) = iwork(niteri) + maxit
         else
            iwork(maxiti) = iwork(niteri) + 10
         end if

         if (iwork(niteri) < iwork(maxiti)) then
            info = 0
         end if

         if (job >= 0) iwork(jobi) = job
         if (iprint >= 0) iwork(iprini) = iprint
         if (partol >= zero .and. partol < one) work(partli) = partol
         if (sstol >= zero .and. sstol < one) work(sstoli) = sstol

         work(olmavi) = work(olmavi)*iwork(niteri)

         if (implct) then
            call dcopy(n*nq, work(fni), 1, work(fi), 1)
         else
            work(fi:fi + (n*nq - 1)) = &
               work(fni:fni + (n*nq - 1)) - reshape(y(1:n, :), shape=[n*nq])
         end if
         call dwght(n, nq, &
                    reshape(work(we1i:we1i + ldwe*ld2we*nq - 1), [ldwe, ld2we, nq]), &
                    ldwe, ld2we, &
                    reshape(work(fi:fi + n*nq - 1), [n, nq]), tempret(1:n, 1:nq))
         work(fi:fi + n*nq - 1) = reshape(tempret(1:n, 1:nq), [n*nq])
         work(wssepi) = ddot(n*nq, work(fi), 1, work(fi), 1)
         work(wssi) = work(wssepi) + work(wssdei)

      else

         ! Perform error checking
         info = 0
         call dodchk(n, m, np, nq, &
                     isodr, anajac, implct, &
                     beta, ifixb, &
                     ldx, ldifx, ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, &
                     ldy, &
                     lwork, lwkmn, liwork, liwkmn, &
                     sclb, scld, stpb, stpd, &
                     info, &
                     lower, upper)
         if (info > 0) then
            goto 50
         end if

         ! Initialize work vectors as necessary
         do i = n*m + n*nq + 1, lwork
            work(i) = zero
         end do
         do i = 1, liwork
            iwork(i) = 0
         end do

         call diniwk(n, m, np, &
                     work, lwork, iwork, liwork, &
                     x, ldx, ifixx, ldifx, scld, ldscld, &
                     beta, sclb, &
                     sstol, partol, maxit, taufac, &
                     job, iprint, lunerr, lunrpt, &
                     lower, upper, &
                     epsmai, sstoli, partli, maxiti, taufci, &
                     jobi, iprini, luneri, lunrpi, &
                     ssfi, tti, ldtti, deltai, &
                     loweri, upperi, boundi)

         iwork(msgb) = -1
         iwork(msgd) = -1
         work(taui) = -work(taufci)

         ! Set up for parameter estimation -
         ! Pull BETA's to be estimated and corresponding scale values
         ! and store in WORK(BETACI) and WORK(SSI), respectively
         call dpack(np, iwork(nppi), work(betaci), beta, ifixb)
         call dpack(np, iwork(nppi), work(ssi), work(ssfi), ifixb)
         npp = iwork(nppi)

         ! Check that WD is positive definite and WE is positive semidefinite,
         ! saving factorization of WE, and counting number of nonzero weights
         call dfctrw(n, m, nq, npp, &
                     isodr, &
                     we, ldwe, ld2we, wd, ldwd, ld2wd, &
                     work(wrk2i), work(wrk4i), &
                     work(we1i), nnzw, info)
         iwork(nnzwi) = nnzw

         if (info /= 0) then
            goto 50
         end if

         ! Evaluate the predicted values and weighted EPSILONS at the starting point
         call dunpac(np, work(betaci), beta, ifixb)
         work(xplusi:xplusi + (n*m - 1)) = &
            work(deltai:deltai + (n*m - 1)) + reshape(x(1:n, :), shape=[n*m])
         istop = 0
         call fcn(n, m, np, nq, &
                  n, m, np, &
                  beta, work(xplusi), &
                  ifixb, ifixx, ldifx, &
                  002, work(fni), work(wrk6i), work(wrk1i), &
                  istop)
         iwork(istopi) = istop
         if (istop == 0) then
            iwork(nfevi) = iwork(nfevi) + 1
            if (implct) then
               call dcopy(n*nq, work(fni), 1, work(fi), 1)
            else
               work(fi:fi + (n*nq - 1)) = &
                  work(fni:fni + (n*nq - 1)) - reshape(y(1:n, :), shape=[n*nq])
            end if
            call dwght(n, nq, &
                       reshape(work(we1i:we1i + ldwe*ld2we*nq - 1), [ldwe, ld2we, nq]), &
                       ldwe, ld2we, &
                       reshape(work(fi:fi + n*nq - 1), [n, nq]), tempret(1:n, 1:nq))
            work(fi:fi + n*nq - 1) = reshape(tempret(1:n, 1:nq), [n*nq])
         else
            info = 52000
            goto 50
         end if

         ! Compute norm of the initial estimates
         call dwght(npp, 1, &
                    reshape(work(ssi:ssi + npp - 1), [npp, 1, 1]), &
                    npp, 1, &
                    reshape(work(betaci:betaci + npp - 1), [npp, 1]), tempret(1:npp, 1:1))
         work(wrk:wrk + npp - 1) = tempret(1:npp, 1)
         if (isodr) then
            call dwght(n, m, &
                       reshape(work(tti:tti + iwork(ldtti)*1*m - 1), [iwork(ldtti), 1, m]), &
                       iwork(ldtti), 1, &
                       reshape(work(deltai:deltai + n*m - 1), [n, m]), tempret(1:n, 1:m))
            work(wrk + npp:wrk + npp + n*m - 1) = reshape(tempret(1:n, 1:m), [n*m])
            work(pnormi) = dnrm2(npp + n*m, work(wrk), 1)
         else
            work(pnormi) = dnrm2(npp, work(wrk), 1)
         end if

         ! Compute sum of squares of the weighted EPSILONS and weighted DELTAS
         work(wssepi) = ddot(n*nq, work(fi), 1, work(fi), 1)
         if (isodr) then
            call dwght(n, m, wd, ldwd, ld2wd, &
                       reshape(work(deltai:deltai + n*m), [n, m]), &
                       tempret(1:n, 1:m))
            work(wrk:wrk + n*m - 1) = reshape(tempret(1:n, 1:m), [n*m])
            work(wssdei) = ddot(n*m, work(deltai), 1, work(wrk), 1)
         else
            work(wssdei) = zero
         end if
         work(wssi) = work(wssepi) + work(wssdei)

         ! Select first row of X + DELTA that contains no zeros
         nrow = -1
         call dsetn(n, m, work(xplusi), n, nrow)
         iwork(nrowi) = nrow

         ! Set number of good digits in function results
         epsmac = work(epsmai)
         if (ndigit < 2) then
            iwork(netai) = -1
            nfev = iwork(nfevi)
            call detaf(fcn, &
                       n, m, np, nq, &
                       work(xplusi), beta, epsmac, nrow, &
                       work(betani), work(fni), &
                       ifixb, ifixx, ldifx, &
                       istop, nfev, eta, neta, &
                       work(wrk1i), work(wrk2i), work(wrk6i), work(wrk7i), &
                       info, &
                       lower, upper)
            iwork(istopi) = istop
            iwork(nfevi) = nfev
            if (istop /= 0 .or. info /= 0) then
               if (info == 0) then
                  info = 53000
               end if
               iwork(netai) = 0
               work(etai) = zero
               goto 50
            else
               iwork(netai) = -neta
               work(etai) = eta
            end if
         else
            iwork(netai) = min(ndigit, int(p5 - log10(epsmac)))
            work(etai) = max(epsmac, ten**(-ndigit))
         end if

         ! Check bounds are large enough for derivative calculations.
         if (.not. anajac .or. chkjac) then
            if (cdjac) then
               do k = 1, np
                  if (upper(k) - abs(2*derstep(1, k, upper(k), work(ssfi), stpb, neta)) &
                      < lower(k)) then
                     info = 90020
                     goto 50
                  end if
               end do
            else
               do k = 1, np
                  if (upper(k) - abs(2*derstep(0, k, upper(k), work(ssfi), stpb, neta)) &
                      < lower(k)) then
                     info = 90020
                     goto 50
                  end if
               end do
            end if
         end if

         ! Check derivatives if necessary
         if (chkjac .and. anajac) then
            ntol = -1
            nfev = iwork(nfevi)
            njev = iwork(njevi)
            neta = iwork(netai)
            ldtt = iwork(ldtti)
            eta = work(etai)
            epsmac = work(epsmai)

            ! Ensure beta is not too close to bounds for the derivative check
            betaj(:) = beta(:)
            call mbfb(np, betaj, lower, upper, work(ssfi), stpb, neta, eta, interval)

            ! Check the derivatives
            call djck(fcn, &
                      n, m, np, nq, &
                      beta, betaj, work(xplusi), &
                      ifixb, ifixx, ldifx, stpb, stpd, ldstpd, &
                      work(ssfi), work(tti), ldtt, &
                      eta, neta, ntol, nrow, isodr, epsmac, &
                      work(fni), work(fjacbi), work(fjacdi), &
                      iwork(msgb), iwork(msgd), work(diffi), &
                      istop, nfev, njev, &
                      work(wrk1i), work(wrk2i), work(wrk6i), &
                      interval)
            iwork(istopi) = istop
            iwork(nfevi) = nfev
            iwork(njevi) = njev
            iwork(ntoli) = ntol
            if (istop /= 0) then
               info = 54000
            elseif (iwork(msgb) /= 0 .or. iwork(msgd) /= 0) then
               info = 40000
            end if
         else
            ! Indicate user supplied derivatives were not checked
            iwork(msgb) = -1
            iwork(msgd) = -1
         end if

         ! Print appropriate error messages
50       if ((info /= 0) .or. &
             (iwork(msgb) /= -1)) then
            if (lunerr /= 0 .and. iprint /= 0) then
               call dodper &
                  (info, lunerr, &
                   n, m, np, nq, &
                   ldscld, ldstpd, ldwe, ld2we, ldwd, ld2wd, &
                   lwkmn, liwkmn, &
                   work(fjacbi), work(fjacdi), &
                   work(diffi), iwork(msgb), isodr, iwork(msgd), &
                   work(xplusi), iwork(nrowi), iwork(netai), iwork(ntoli))
            end if

            ! Set INFO to reflect errors in the user supplied Jacobians
            if (info == 40000) then
               if (iwork(msgb) == 2 .or. iwork(msgd) == 2) then
                  if (iwork(msgb) == 2) then
                     info = info + 1000
                  end if
                  if (iwork(msgd) == 2) then
                     info = info + 100
                  end if
               else
                  info = 0
               end if
            end if
            if (info /= 0) then
               return
            end if
         end if
      end if

      ! Save the initial values of BETA
      call dcopy(np, beta, 1, work(beta0i), 1)

      ! Find least squares solution
      call dcopy(n*nq, work(fni), 1, work(fsi), 1)
      ldtt = iwork(ldtti)
      call dodmn(head, fstitr, prtpen, &
                 fcn, n, m, np, nq, job, beta, y, ldy, x, ldx, &
                 we, work(we1i), ldwe, ld2we, wd, ldwd, ld2wd, &
                 ifixb, ifixx, ldifx, &
                 work(betaci), work(betani), work(betasi), work(si), &
                 work(deltai), work(deltni), work(deltsi), &
                 work(loweri), work(upperi), &
                 work(ti), work(fi), work(fni), work(fsi), &
                 work(fjacbi), iwork(msgb), work(fjacdi), iwork(msgd), &
                 work(ssfi), work(ssi), work(tti), ldtt, &
                 stpb, stpd, ldstpd, &
                 work(xplusi), work(wrk), lwrk, &
                 work, lwork, tempret, iwork, liwork, info, &
                 iwork(boundi))
      maxit1 = iwork(maxiti) - iwork(niteri)
      tstimp = zero
      do k = 1, np
         if (beta(k) == zero) then
            tstimp = max(tstimp, abs(beta(k) - work(beta0i - 1 + k))/work(ssfi - 1 + k))
         else
            tstimp = max(tstimp, abs(beta(k) - work(beta0i - 1 + k))/abs(beta(k)))
         end if
      end do

   end subroutine doddrv