lcstep Subroutine

   subroutine lcstep &
      (n, m, np, q, npp, &
       f, fjacb, fjacd, &
       wd, ldwd, ld2wd, ss, tt, ldtt, delta, &
       alpha, epsfcn, isodr, &
       tfjacb, omega, u, qraux, kpvt, &
       s, t, phi, irank, rcond, forvcv, &
       wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc)
   !! Compute locally constrained steps `s` and `t`, and `phi(alpha)`.
      ! @note: This is one of the most time-consuming subroutines in ODRPACK (~25% of total).

      use odrpack_kinds, only: zero, one
      use linpack, only: dchex, dqrdc, dqrsl, dtrco, dtrsl
      use blas_interfaces, only: drot, drotg

      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)
         !! Squared `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 for the unfixed `beta`s.
      real(wp), intent(in) :: tt(ldtt, m)
         !! Scaling values for `delta`.
      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(in) :: alpha
         !! Levenberg-Marquardt parameter.
      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 defined such that:
         !! `omega*trans(omega) = inv(I + fjacd*inv(e)*trans(fjacd))
         !! = (I - fjacd*inv(p)*trans(fjacd))`
         !! where `e = d**2 + alpha*tt**2` and
         !! `p = trans(fjacd)*fjacd + d**2 + alpha*tt**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 Q-R decomposition.
      integer, intent(out) :: kpvt(np)
         !! Pivot vector.
      real(wp), intent(out) :: s(np)
         !! Step for `beta`.
      real(wp), intent(out) :: t(n, m)
         !! Step for `delta`.
      real(wp), intent(out) :: phi
         !! Difference between the norm of the scaled step and the trust region diameter.
      integer, intent(out) :: irank
         !! Rank deficiency of the Jacobian wrt `beta`.
      real(wp), intent(out) :: rcond
         !! Approximate reciprocal condition number of `tfjacb`.
      logical, intent(in) :: forvcv
         !! Variable designating whether this subroutine was called to set up for the
         !! covariance matrix computations (`forvcv = .true.`) or not (`forvcv = .false.`).
      real(wp), intent(out) :: wrk1(n, q, m)
         !! A work array of `(n, q, m)` elements.
      real(wp), intent(out) :: wrk2(n, q)
         !! A work array of `(n, q)` elements.
      real(wp), intent(out) :: wrk3(np)
         !! A work array of `(np)` elements.
      real(wp), intent(out) :: wrk4(m, m)
         !! A work array of `(m, m)` elements.
      real(wp), intent(out) :: wrk5(m)
         !! A work array of `(m)` elements.
      real(wp), intent(out) :: wrk(lwrk)
         !! A work array of `(lwrk)` elements, _equivalenced_ to `wrk1` and `wrk2`.
      integer, intent(in) :: lwrk
         !! Length of vector `wrk`.
      integer, intent(inout) :: istopc
         !! Variable designating whether the computations were stopped due to a numerical
         !! error within subroutine `dodstp`.

      ! Local scalars
      real(wp) :: co, si, temp
      integer :: i, imax, inf, ipvt, j, k, k1, k2, kp, l
      logical :: elim

      ! Local arrays
      real(wp) :: dum(2)

      ! Variable definitions (alphabetically)
      !  CO:      The cosine from the plane rotation.
      !  DUM:     A dummy array.
      !  ELIM:    The variable designating whether columns of the Jacobian
      !  I:       An indexing variable.
      !  IMAX:    The index of the element of U having the largest absolute value.
      !  INF:     The return code from LINPACK routines.
      !  IPVT:    The variable designating whether pivoting is to be done.
      !  J:       An indexing variable.
      !  K:       An indexing variable.
      !  K1:      An indexing variable.
      !  K2:      An indexing variable.
      !  KP:      The rank of the Jacobian wrt BETA.
      !  L:       An indexing variable.
      !  SI:      The sine from the plane rotation.
      !  TEMP:    A temporary storage LOCATION.

      ! Compute loop parameters which depend on weight structure

      !  Set up KPVT if ALPHA = 0
      if (alpha == zero) then
         kp = npp
         do k = 1, np
            kpvt(k) = k
         end do
      else
         if (npp >= 1) then
            kp = npp - irank
         else
            kp = npp
         end if
      end if

      if (isodr) then
         ! T = WD * DELTA = D*G2
         call scale_mat(n, m, wd, ldwd, ld2wd, delta, t)

         do i = 1, n

            !  Compute WRK4, such that TRANS(WRK4)*WRK4 = E = (D**2 + ALPHA*TT**2)
            call esubi(n, m, wd, ldwd, ld2wd, alpha, tt, ldtt, i, wrk4)
            call fctr(.false., wrk4, m, m, inf)
            if (inf /= 0) then
               istopc = 60000
               return
            end if
            ! Compute OMEGA, such that
            ! trans(OMEGA)*OMEGA = I+FJACD*inv(E)*trans(FJACD)
            ! inv(trans(OMEGA)*OMEGA) = I-FJACD*inv(P)*trans(FJACD)
            call vevtr(m, q, i, fjacd, n, m, wrk4, m, wrk1, n, q, omega, q, wrk5)
            do l = 1, q
               omega(l, l) = one + omega(l, l)
            end do
            call fctr(.false., omega, q, q, inf)
            if (inf /= 0) then
               istopc = 60000
               return
            end if
            ! Compute WRK1 = trans(FJACD)*(I-FJACD*inv(P)*trans(JFACD))
            !              = trans(FJACD)*inv(trans(OMEGA)*OMEGA)
            do j = 1, m
               wrk1(i, :, j) = fjacd(i, j, :)
               call solve_trl(q, omega, q, wrk1(i, 1:q, j), 4)
               call solve_trl(q, omega, q, wrk1(i, 1:q, j), 2)
            end do

            ! Compute WRK5 = inv(E)*D*G2
            wrk5 = t(i, :)
            call solve_trl(m, wrk4, m, wrk5, 4)
            call solve_trl(m, wrk4, m, wrk5, 2)

            ! Compute TFJACB = inv(trans(OMEGA))*FJACB
            do k = 1, kp
               tfjacb(i, :, k) = fjacb(i, kpvt(k), :)
               call solve_trl(q, omega, q, tfjacb(i, 1:q, k), 4)
               if (ss(1) > zero) then
                  tfjacb(i, :, k) = tfjacb(i, :, k)/ss(kpvt(k))
               else
                  tfjacb(i, :, k) = tfjacb(i, :, k)/abs(ss(1))
               end if
            end do

            ! Compute WRK2 = (V*inv(E)*D**2*G2 - G1)
            do l = 1, q
               wrk2(i, l) = dot_product(fjacd(i, :, l), wrk5) - f(i, l)
            end do

            ! Compute WRK2 = inv(trans(OMEGA))*(V*inv(E)*D**2*G2 - G1)
            call solve_trl(q, omega, q, wrk2(i, 1:q), 4)

         end do

      else
         do l = 1, q
            do i = 1, n
               do k = 1, kp
                  tfjacb(i, l, k) = fjacb(i, kpvt(k), l)
                  if (ss(1) > zero) then
                     tfjacb(i, l, k) = tfjacb(i, l, k)/ss(kpvt(k))
                  else
                     tfjacb(i, l, k) = tfjacb(i, l, k)/abs(ss(1))
                  end if
               end do
               wrk2(i, l) = -f(i, l)
            end do
         end do
      end if

      ! Compute S
      ! Do QR factorization (with column pivoting of TFJACB if ALPHA = 0)
      if (alpha == zero) then
         ipvt = 1
         kpvt = 0
      else
         ipvt = 0
      end if

      call dqrdc(tfjacb, n*q, n*q, kp, qraux, kpvt, wrk3, ipvt)
      call dqrsl(tfjacb, n*q, n*q, kp, qraux, wrk2, dum, wrk2, dum, dum, dum, 1000, inf)
      if (inf /= 0) then
         istopc = 60000
         return
      end if

      ! Eliminate alpha part using givens rotations
      if (alpha /= zero) then
         s(1:npp) = zero
         do k1 = 1, kp
            wrk3(1:kp) = zero
            wrk3(k1) = sqrt(alpha)
            do k2 = k1, kp
               call drotg(tfjacb(k2, 1, k2), wrk3(k2), co, si)
               if (kp - k2 >= 1) then
                  call drot(kp - k2, tfjacb(k2, 1, k2 + 1), n*q, wrk3(k2 + 1), 1, co, si)
               end if
               temp = co*wrk2(k2, 1) + si*s(kpvt(k1))
               s(kpvt(k1)) = -si*wrk2(k2, 1) + co*s(kpvt(k1))
               wrk2(k2, 1) = temp
            end do
         end do
      end if

      ! Compute solution - eliminate variables if necessary
      if (npp >= 1) then
         if (alpha == zero) then
            kp = npp
            elim = .true.
            do while (elim .and. kp >= 1)
               ! Estimate RCOND - U will contain approx null vector
               call dtrco(tfjacb, n*q, kp, rcond, u, 1)
               if (rcond <= epsfcn) then
                  elim = .true.
                  imax = maxloc(u(1:kp), dim=1)
                  ! IMAX is the column to remove - use DCHEX and fix KPVT
                  if (imax /= kp) then
                     call dchex(tfjacb, n*q, kp, imax, kp, wrk2, n*q, 1, qraux, wrk3, 2)
                     k = kpvt(imax)
                     do i = imax, kp - 1
                        kpvt(i) = kpvt(i + 1)
                     end do
                     kpvt(kp) = k
                  end if
                  kp = kp - 1
               else
                  elim = .false.
               end if
            end do
            irank = npp - kp
         end if
      end if

      if (forvcv) return

      ! Backsolve and unscramble
      if (npp >= 1) then
         do i = kp + 1, npp
            wrk2(i, 1) = zero
         end do
         if (kp >= 1) then
            call dtrsl(tfjacb, n*q, kp, wrk2, 01, inf)
            if (inf /= 0) then
               istopc = 60000
               return
            end if
         end if
         do i = 1, npp
            if (ss(1) > zero) then
               s(kpvt(i)) = wrk2(i, 1)/ss(kpvt(i))
            else
               s(kpvt(i)) = wrk2(i, 1)/abs(ss(1))
            end if
         end do
      end if

      if (isodr) then
         !  NOTE: T and WRK1 have been initialized above,
         !          where T    = WD * DELTA = D*G2
         !          WRK1 = trans(FJACD)*(I-FJACD*inv(P)*trans(JFACD))
         do i = 1, n

            ! Compute WRK4, such that trans(WRK4)*WRK4 = E = (D**2 + ALPHA*TT**2)
            call esubi(n, m, wd, ldwd, ld2wd, alpha, tt, ldtt, i, wrk4)
            call fctr(.false., wrk4, m, m, inf)
            if (inf /= 0) then
               istopc = 60000
               return
            end if

            ! Compute WRK5 = inv(E)*D*G2
            wrk5 = t(i, :)
            call solve_trl(m, wrk4, m, wrk5, 4)
            call solve_trl(m, wrk4, m, wrk5, 2)

            do l = 1, q
               wrk2(i, l) = f(i, l) &
                            + dot_product(fjacb(i, 1:npp, l), s(1:npp)) &
                            - dot_product(fjacd(i, :, l), wrk5)
            end do

            do j = 1, m
               wrk5(j) = dot_product(wrk1(i, :, j), wrk2(i, :))
               t(i, j) = -(wrk5(j) + t(i, j))
            end do

            call solve_trl(m, wrk4, m, t(i, 1:m), 4)
            call solve_trl(m, wrk4, m, t(i, 1:m), 2)

         end do

      end if

      ! Compute PHI(ALPHA) from scaled S and T
      call scale_mat(npp, 1, ss, npp, 1, s, wrk(1:npp))
      if (isodr) then
         call scale_mat(n, m, tt, ldtt, 1, t, wrk(npp + 1:npp + 1 + n*m - 1))
         phi = norm2(wrk(1:npp + n*m))
      else
         phi = norm2(wrk(1:npp))
      end if

   end subroutine lcstep