trust_region Subroutine

   subroutine trust_region &
      (n, m, np, q, npp, &
       f, fjacb, fjacd, &
       wd, ldwd, ld2wd, ss, tt, ldtt, delta, &
       alpha2, tau, epsfcn, isodr, &
       tfjacb, omega, u, qraux, jpvt, &
       s, t, nlms, rcond, irank, &
       wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc)
   !! Compute Levenberg-Marquardt parameter and steps `s` and `t` using analog of the
   !! trust-region Levenberg-Marquardt algorithm.

      use odrpack_kinds, only: zero

      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.
      integer, intent(in) :: npp
         !! Number of function parameters being estimated.
      real(wp), intent(in) :: f(n, q)
         !! Weighted estimated values of `epsilon`.
      real(wp), intent(in) :: fjacb(n, np, q)
         !! Jacobian with respect to `beta`.
      real(wp), intent(in) :: fjacd(n, m, q)
         !! Jacobian with respect to `delta`.
      real(wp), intent(in) :: wd(ldwd, ld2wd, m)
         !! `delta` weights.
      integer, intent(in) :: ldwd
         !! Leading dimension of array `wd`.
      integer, intent(in) :: ld2wd
         !! Second dimension of array `wd`.
      real(wp), intent(in) :: ss(np)
         !! Scaling values used for the unfixed `beta`s.
      real(wp), intent(in) :: tt(ldtt, m)
         !! Scale used for the `delta`s.
      integer, intent(in) :: ldtt
         !! Leading dimension of array `tt`.
      real(wp), intent(in) :: delta(n, m)
         !! Estimated errors in the explanatory variables.
      real(wp), intent(inout) :: alpha2
         !! Current Levenberg-Marquardt parameter.
      real(wp), intent(inout) :: tau
         !! Trust region diameter.
      real(wp), intent(in) :: epsfcn
         !! Function's precision.
      logical, intent(in) :: isodr
         !! Variable designating whether the solution is by ODR (`.true.`) or by OLS (`.false.`).
      real(wp), intent(out) :: tfjacb(n, q, np)
         !! Array `omega*fjacb`.
      real(wp), intent(out) :: omega(q, q)
         !! Array `(I-fjacd*inv(p)*trans(fjacd))**(-1/2)`.
      real(wp), intent(out) :: u(np)
         !! Approximate null vector for `tfjacb`.
      real(wp), intent(out) :: qraux(np)
         !! Array required to recover the orthogonal part of the QR decomposition.
      integer, intent(out) :: jpvt(np)
         !! Pivot vector.
      real(wp), intent(out) :: s(np)
         !! Step for `beta`.
      real(wp), intent(out) :: t(n, m)
         !! Step for `delta`.
      integer, intent(out) :: nlms
         !! Number of Levenberg-Marquardt steps taken.
      real(wp), intent(out) :: rcond
         !! Approximate reciprocal condition of `tfjacb`.
      integer, intent(out) :: irank
         !! Aank deficiency of the Jacobian wrt `beta`.
      real(wp), intent(out) :: wrk1(n, q, m)
         !! Work array of `(n, q, m)` 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) :: wrk4(m, m)
         !! Work array of `(m, m)` elements.
      real(wp), intent(out) :: wrk5(m)
         !! Work array of `(m)` elements.
      real(wp), intent(out) :: wrk(lwrk)
         !! Work array of `(lwrk)` elements, _equivalenced_ to `wrk1` and `wrk2`.
      integer, intent(in) :: lwrk
         !! Length of vector `wrk`.
      integer, intent(out) :: istopc
         !! Variable designating whether the computations were stopped due to some other
         !! numerical error detected within subroutine `dodstp`.

      ! Local scalars
      real(wp), parameter :: p001 = 0.001_wp, p1 = 0.1_wp
      real(wp) :: alpha1, alphan, bot, phi1, phi2, sa, top
      integer :: i, iwrk, j, k
      logical :: forvcv

      ! Variable Definitions (alphabetically)
      !  ALPHAN:  The new Levenberg-Marquardt parameter.
      !  ALPHA1:  The previous Levenberg-Marquardt parameter.
      !  BOT:     The lower limit for setting ALPHA.
      !  FORVCV:  The variable designating whether this subroutine was called to set up for the
      !           covariance matrix computations (TRUE) or not (FALSE).
      !  I:       An indexing variable.
      !  IWRK:    An indexing variable.
      !  J:       An indexing variable.
      !  K:       An indexing variable.
      !  PHI1:    The previous difference between the norm of the scaled step and the trust
      !           region diameter.
      !  PHI2:    The current difference between the norm of the scaled step and the trust region
      !           diameter.
      !  SA:      The scalar PHI2*(ALPHA1-ALPHA2)/(PHI1-PHI2).
      !  TOP:     The upper limit for setting ALPHA.

      forvcv = .false.
      istopc = 0

      ! Compute full Gauss-Newton step (ALPHA=0)
      alpha1 = zero
      call lcstep(n, m, np, q, npp, &
                  f, fjacb, fjacd, &
                  wd, ldwd, ld2wd, ss, tt, ldtt, delta, &
                  alpha1, epsfcn, isodr, &
                  tfjacb, omega, u, qraux, jpvt, &
                  s, t, phi1, irank, rcond, forvcv, &
                  wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc)
      if (istopc /= 0) then
         return
      end if

      ! Initialize TAU if necessary
      if (tau < zero) then
         tau = abs(tau)*phi1
      end if

      ! Check if full Gauss-Newton step is optimal
      if ((phi1 - tau) <= p1*tau) then
         nlms = 1
         alpha2 = zero
         return
      end if

      ! Full Gauss-Newton step is outside trust region - find locally constrained optimal step
      phi1 = phi1 - tau

      ! Initialize upper and lower bounds for ALPHA
      bot = zero

      do k = 1, npp
         tfjacb(1:n, 1:q, k) = fjacb(1:n, k, 1:q)
         wrk(k) = sum(tfjacb(:, :, k)*f)
      end do

      call scale_vec(npp, 1, ss, npp, wrk, npp, wrk, npp) ! work is input (as t) and output (as sclt)

      if (isodr) then
         call scale_mat(n, m, wd, ldwd, ld2wd, delta, wrk(npp + 1:npp + 1 + n*m - 1))
         iwrk = npp
         do j = 1, m
            do i = 1, n
               iwrk = iwrk + 1
               wrk(iwrk) = wrk(iwrk) + dot_product(fjacd(i, j, :), f(i, :))
            end do
         end do
         call scale_vec(n, m, tt, ldtt, wrk(npp + 1), n, wrk(npp + 1), n)
         top = norm2(wrk(1:npp + n*m))/tau
      else
         top = norm2(wrk(1:npp))/tau
      end if

      if (alpha2 > top .or. alpha2 == zero) then
         alpha2 = p001*top
      end if

      ! Main loop

      do i = 1, 10

         ! Compute locally constrained steps S and T and PHI(ALPHA) for current value of ALPHA
         call lcstep(n, m, np, q, npp, &
                     f, fjacb, fjacd, &
                     wd, ldwd, ld2wd, ss, tt, ldtt, delta, &
                     alpha2, epsfcn, isodr, &
                     tfjacb, omega, u, qraux, jpvt, &
                     s, t, phi2, irank, rcond, forvcv, &
                     wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc)
         if (istopc /= 0) then
            return
         end if
         phi2 = phi2 - tau

         ! Check whether current step is optimal
         if (abs(phi2) <= p1*tau .or. (alpha2 == bot .and. phi2 < zero)) then
            nlms = i + 1
            return
         end if

         ! Current step is not optimaL

         ! Update bounds for ALPHA and compute new ALPHA
         if (phi1 - phi2 == zero) then
            nlms = 12
            return
         end if
         sa = phi2*(alpha1 - alpha2)/(phi1 - phi2)
         if (phi2 < zero) then
            top = min(top, alpha2)
         else
            bot = max(bot, alpha2)
         end if
         if (phi1*phi2 > zero) then
            bot = max(bot, alpha2 - sa)
         else
            top = min(top, alpha2 - sa)
         end if

         alphan = alpha2 - sa*(phi1 + tau)/tau
         if (alphan >= top .or. alphan <= bot) then
            alphan = max(p001*top, sqrt(top*bot))
         end if

         ! Get ready for next iteration
         alpha1 = alpha2
         alpha2 = alphan
         phi1 = phi2

      end do

      ! Set NLMS to indicate an optimal step could not be found in 10 trys
      nlms = 12

   end subroutine trust_region