fctr Subroutine

   pure subroutine fctr(oksemi, a, lda, n, info)
   !! Factor the positive (semi)definite matrix `a` using a modified Cholesky factorization.
      ! Adapted from LINPACK subroutine DPOFA.

      use odrpack_kinds, only: zero

      logical, intent(in) :: oksemi
         !! Flag indicating whether the factored array can be positive semidefinite
         !! (`.true.`) or whether it must be found to be positive definite (`.false.`).
      real(wp), intent(inout) :: a(lda, n)
         !! Array to be factored. Upon return, `a` contains the upper triangular matrix
         !! `r` so that `a = trans(r)*r` where the strict lower triangle is set to zero.
         !! If `info /= 0`, the factorization is not complete.
      integer, intent(in) :: lda
         !! Leading dimension of array `a`.
      integer, intent(in) :: n
         !! Number of rows and columns of data in array `a`.
      integer, intent(out) :: info
         !! Output flag.
         !! `0`: Factorization was completed.
         !! `k`: Error condition. The leading minor of order `k` is not positive (semi)definite.

      ! Local scalars
      real(wp) :: xi, s, t
      integer j, k

      ! Variable Definitions (alphabetically)
      !  J:       An indexing variable.
      !  K:       An indexing variable.
      !  XI:      A value used to test for non positive semidefiniteness.

      ! Set relative tolerance for detecting non positive semidefiniteness.
      xi = -10*epsilon(zero)

      ! Compute factorization, storing in upper triangular portion of A
      do j = 1, n
         info = j
         s = zero
         do k = 1, j - 1
            if (a(k, k) == zero) then
               t = zero
            else
               t = a(k, j) - dot_product(a(1:k - 1, k), a(1:k - 1, j))
               t = t/a(k, k)
            end if
            a(k, j) = t
            s = s + t**2
         end do
         s = a(j, j) - s

         ! Exit
         if (a(j, j) < zero .or. s < xi*abs(a(j, j))) then
            return
         elseif (.not. oksemi .and. s <= zero) then
            return
         elseif (s <= zero) then
            a(j, j) = zero
         else
            a(j, j) = sqrt(s)
         end if
      end do
      info = 0

      ! Zero out lower portion of A
      do j = 1, n - 1
         a(j + 1:n, j) = zero
      end do

   end subroutine fctr