vcv_beta Subroutine

   subroutine vcv_beta &
      (n, m, np, q, npp, &
       f, fjacb, fjacd, &
       wd, ldwd, ld2wd, ssf, ss, tt, ldtt, delta, &
       epsfcn, isodr, &
       vcv, sd, &
       wrk6, omega, u, qraux, jpvt, &
       s, t, irank, rcond, rss, idf, rvar, ifixb, &
       wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc)
   !! Compute covariance matrix of estimated parameters.

      use odrpack_kinds, only: zero
      use linpack, only: dpodi

      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) :: ssf(np)
         !! Scaling values used for `beta`.
      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) :: 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) :: vcv(np, np)
         !! Covariance matrix of the estimated `beta`s.
      real(wp), intent(out) :: sd(np)
         !! Standard deviations of the estimated `beta`s.
      real(wp), intent(out) :: wrk6(n*q, np)
         !! A work array of `(n*q, np)` elements.
      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))`.
      real(wp), intent(out) :: u(np)
         !! Approximate null vector for `fjacb`.
      real(wp), intent(out) :: qraux(np)
         !! Array required to recover the orthogonal part of the Q-R 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) :: irank
         !! Rank deficiency of the Jacobian wrt `beta`.
      real(wp), intent(out) :: rcond
         !! Approximate reciprocal condition of `fjacb`.
      real(wp), intent(inout) :: rss
         !! Residual sum of squares.
      integer, intent(out) :: idf
         !! Degrees of freedom of the fit, equal to the number of observations with nonzero
         !! weighted derivatives minus the number of parameters being estimated.
      real(wp), intent(out) :: rvar
         !! Residual variance.
      integer, intent(in) :: ifixb(np)
         !! Values designating whether the elements of `beta` are fixed at their input values or not.
      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 `lwrk`.
      integer, intent(out) :: istopc
         !! Variable designating whether the computations were stoped due to a numerical
         !! error within subroutine `dodstp`.

      ! Local scalars
      real(wp) :: temp
      integer :: i, iunfix, j, junfix, kp
      logical :: forvcv

      ! Variable definitions (alphabetically)
      !  FORVCV:  The variable designating whether subroutine DODSTP is called to set up for
      !           the covariance matrix computations (TRUE) or not (FALSE).
      !  I:       An indexing variable.
      !  IUNFIX:  The index of the next unfixed parameter.
      !  J:       An indexing variable.
      !  JUNFIX:  The index of the next unfixed parameter.
      !  KP:      The rank of the Jacobian wrt BETA.
      !  TEMP:    A temporary storage location

      forvcv = .true.
      istopc = 0

      call lcstep(n, m, np, q, npp, &
                  f, fjacb, fjacd, &
                  wd, ldwd, ld2wd, ss, tt, ldtt, delta, &
                  zero, epsfcn, isodr, &
                  wrk6, omega, u, qraux, jpvt, &
                  s, t, temp, irank, rcond, forvcv, &
                  wrk1, wrk2, wrk3, wrk4, wrk5, wrk, lwrk, istopc)
      if (istopc /= 0) then
         return
      end if

      kp = npp - irank
      call dpodi(wrk6, n*q, kp, wrk3, 1)

      idf = 0
      do i = 1, n
         if (any(fjacb(i, :, :) /= zero)) then
            idf = idf + 1
            cycle
         end if
         if (isodr) then
            if (any(fjacd(i, :, :) /= zero)) then
               idf = idf + 1
               cycle
            end if
         end if
      end do

      if (idf > kp) then
         idf = idf - kp
         rvar = rss/idf
      else
         idf = 0
         rvar = rss
      end if

      ! Store variances in SD, restoring original order
      sd = zero
      do i = 1, kp
         sd(jpvt(i)) = wrk6(i, i)
      end do

      if (np > npp) then
         junfix = npp
         do j = np, 1, -1
            if (ifixb(j) == 0) then
               sd(j) = zero
            else
               sd(j) = sd(junfix)
               junfix = junfix - 1
            end if
         end do
      end if

      ! Store covariance matrix in VCV, restoring original order
      vcv = zero
      do i = 1, kp
         do j = i + 1, kp
            if (jpvt(i) > jpvt(j)) then
               vcv(jpvt(i), jpvt(j)) = wrk6(i, j)
            else
               vcv(jpvt(j), jpvt(i)) = wrk6(i, j)
            end if
         end do
      end do

      if (np > npp) then
         iunfix = npp
         do i = np, 1, -1
            if (ifixb(i) == 0) then
               do j = i, 1, -1
                  vcv(i, j) = zero
               end do
            else
               junfix = npp
               do j = np, 1, -1
                  if (ifixb(j) == 0) then
                     vcv(i, j) = zero
                  else
                     vcv(i, j) = vcv(iunfix, junfix)
                     junfix = junfix - 1
                  end if
               end do
               iunfix = iunfix - 1
            end if
         end do
      end if

      do i = 1, np
         vcv(i, i) = sd(i)
         sd(i) = sqrt(rvar*sd(i))
         do j = 1, i
            vcv(j, i) = vcv(i, j)
         end do
      end do

      ! Unscale standard errors and covariance matrix
      do i = 1, np
         if (ssf(1) > zero) then
            sd(i) = sd(i)/ssf(i)
         else
            sd(i) = sd(i)/abs(ssf(1))
         end if
         do j = 1, np
            if (ssf(1) > zero) then
               vcv(i, j) = vcv(i, j)/(ssf(i)*ssf(j))
            else
               vcv(i, j) = vcv(i, j)/(ssf(1)*ssf(1))
            end if
         end do
      end do

   end subroutine vcv_beta