init_work Subroutine

   pure subroutine init_work &
      (n, m, np, rwork, lrwork, iwork, liwork, &
       x, 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)
   !! Initialize work vectors as necessary.

      use odrpack_kinds, only: zero, one, two, three

      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.
      real(wp), intent(out) :: rwork(lrwork)
         !! Real work space.
      integer, intent(in) :: lrwork
         !! Length of vector `rwork`.
      integer, intent(out) :: iwork(liwork)
         !! Integer work space.
      integer, intent(in) :: liwork
         !! Length of vector `iwork`.
      real(wp), intent(in) :: x(n, m)
         !! Independent variable.
      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) :: scld(ldscld, m)
         !! Scaling values for `delta`.
      integer, intent(in) :: ldscld
         !! Leading dimension of array `scld`.
      real(wp), intent(in) :: beta(np)
         !! Function parameters.
      real(wp), intent(in) :: sclb(np)
         !! Scaling values for `beta`.
      real(wp), intent(in) :: sstol
         !! Sum-of-squares convergence stopping criteria.
      real(wp), intent(in) :: partol
         !! Parameter convergence stopping criteria.
      integer, intent(in) :: maxit
         !! Maximum number of iterations allowed.
      real(wp), intent(in) :: taufac
         !! Factor used to compute the initial trust region diameter.
      integer, intent(in) :: job
         !! Variable controlling problem initialization and computational method.
      integer, intent(in) :: iprint
         !! Print control variable.
      integer, intent(in) :: lunerr
         !! Logical unit number used for error messages.
      integer, intent(in) :: lunrpt
         !! Logical unit number used for computation reports.
      real(wp), intent(in) :: lower(np)
         !! Lower bounds for the function parameters.
      real(wp), intent(in) :: upper(np)
         !! Upper bounds for the function parameters.
      integer, intent(in) :: epsmai
         !! Location in array `rwork` of variable `epsmac`.
      integer, intent(in) :: sstoli
         !! Location in array `rwork` of variable `sstol`.
      integer, intent(in) :: partli
         !! Location in array `rwork` of variable `partol`.
      integer, intent(in) :: maxiti
         !! Location in array `iwork` of variable `maxit`.
      integer, intent(in) :: taufci
         !! Location in array `rwork` of variable `taufac`.
      integer, intent(in) :: jobi
         !! Location in array `iwork` of variable `job`.
      integer, intent(in) :: iprini
         !! Location in array `iwork` of variable `iprint`.
      integer, intent(in) :: luneri
         !! Location in array `iwork` of variable `lunerr`.
      integer, intent(in) :: lunrpi
         !! Location in array `iwork` of variable `lunrpt`.
      integer, intent(in) :: ssfi
         !! Starting location in array `rwork` of array `ssf`.
      integer, intent(in) :: tti
         !! Starting location in array `rwork` of the array `tt`.
      integer, intent(in) :: ldtti
         !! Leading dimension of array `tt`.
      integer, intent(in) :: deltai
         !! Starting location in array `rwork` of array `delta`.
      integer, intent(in) :: loweri
         !! Starting location in array `iwork` of array `lower`.
      integer, intent(in) :: upperi
         !! Starting location in array `iwork` of array `upper`.
      integer, intent(in) :: boundi
         !! Location in array `iwork` of variable `bound`.

      ! Local scalars
      integer :: i, j, istart
      logical :: anajac, cdjac, chkjac, dovcv, implicit, initd, isodr, redoj, restart

      ! Variable Definitions (alphabetically)
      !  ANAJAC:   The variable designating whether the Jacobians are computed by finite differences
      !            (FALSE) or not (TRUE).
      !  CDJAC:    The variable designating whether the Jacobians are computed by central differences
      !            (TRUE) or by forward differences (FALSE).
      !  CHKJAC:   The variable designating whether the user-supplied Jacobians are to be checked
      !            (TRUE) or not (FALSE).
      !  DOVCV:    The variable designating whether the covariance matrix is to be computed (TRUE)
      !            or not (FALSE).
      !  I:        An indexing variable.
      !  IMPLICIT: The variable designating whether the solution is by implicit ODR (TRUE) or
      !            explicit ODR (FALSE).
      !  INITD:    The variable designating whether DELTA is to be initialized to zero (TRUE) or
      !            to the values in the first N by M elements of array RWORK (FALSE).
      !  ISODR:    The variable designating whether the solution is by ODR (TRUE) or by OLS
      !            (FALSE).
      !  J:        An indexing variable.
      !  REDOJ:    The variable designating whether the Jacobian matrix is to be recomputed for the
      !            computation of the covariance matrix (TRUE) or not (FALSE).
      !  RESTART:  The variable designating whether the call is a restart (TRUE) or not
      !            (FALSE).

      call set_flags(job, restart, initd, dovcv, redoj, anajac, cdjac, chkjac, isodr, implicit)

      ! Store value of machine precision in work vector
      rwork(epsmai) = epsilon(zero)

      ! Set tolerance for stopping criteria based on the change in the parameters (see also
      ! subprogram DODCNT)
      if (partol < zero) then
         rwork(partli) = rwork(epsmai)**(two/three)
      else
         rwork(partli) = min(partol, one)
      end if

      ! Set tolerance for stopping criteria based on the change in the sum of squares of the
      ! weighted observational errors
      if (sstol < zero) then
         rwork(sstoli) = sqrt(rwork(epsmai))
      else
         rwork(sstoli) = min(sstol, one)
      end if

      ! Set factor for computing trust region diameter at first iteration
      if (taufac <= zero) then
         rwork(taufci) = one
      else
         rwork(taufci) = min(taufac, one)
      end if

      ! Set maximum number of iterations
      if (maxit < 0) then
         iwork(maxiti) = 50
      else
         iwork(maxiti) = maxit
      end if

      ! Store problem initialization and computational method control variable
      if (job <= 0) then
         iwork(jobi) = 0
      else
         iwork(jobi) = job
      end if

      ! Set print control
      if (iprint < 0) then
         iwork(iprini) = 2001
      else
         iwork(iprini) = iprint
      end if

      ! Set logical unit number for error messages
      iwork(luneri) = lunerr

      ! Set logical unit number for computation reports
      iwork(lunrpi) = lunrpt

      ! Compute scaling for BETA's and DELTA's
      if (sclb(1) <= zero) then
         call scale_beta(np, beta, rwork(ssfi))
      else
         rwork(ssfi:ssfi + (np - 1)) = sclb
      end if

      if (isodr) then
         if (scld(1, 1) <= zero) then
            iwork(ldtti) = n
            call scale_delta(n, m, x, rwork(tti), iwork(ldtti))
         else
            if (ldscld == 1) then
               iwork(ldtti) = 1
               rwork(tti:tti + (m - 1)) = scld(1:m, 1)
            else
               iwork(ldtti) = n
               do j = 1, m
                  istart = tti + (j - 1)*iwork(ldtti)
                  rwork(istart:istart + (n - 1)) = scld(1:n, j)
               end do
            end if
         end if
      end if

      ! Initialize DELTA's as necessary
      if (isodr) then
         if (initd) then
            rwork(deltai:deltai + (n*m - 1)) = zero
         else
            if (ifixx(1, 1) >= 0) then
               if (ldifx == 1) then
                  do j = 1, m
                     if (ifixx(1, j) == 0) then
                        istart = deltai + (j - 1)*n
                        rwork(istart:istart + (n - 1)) = zero
                     end if
                  end do
               else
                  do j = 1, m
                     do i = 1, n
                        if (ifixx(i, j) == 0) then
                           rwork(deltai - 1 + i + (j - 1)*n) = zero
                        end if
                     end do
                  end do
               end if
            end if
         end if
      else
         rwork(deltai:deltai + (n*m - 1)) = zero
      end if

      ! Copy bounds into RWORK
      rwork(loweri:loweri + np - 1) = lower
      rwork(upperi:upperi + np - 1) = upper

      ! Initialize parameters on bounds in IWORK
      iwork(boundi:boundi + np - 1) = 0

   end subroutine init_work