SUBROUTINE chol (a, n, nn, u, nullty, ifault)

!    Algorithm AS6, Applied Statistics, vol.17, (1968)

!    Given a symmetric matrix order n as lower triangle in a( ),
!    calculates an upper triangle, u( ), such that uprime * u = a.
!    a must be positive semi-definite.  eta is set to multiplying
!    factor determining effective zero for pivot.

!    arguments:-
!    a()     = input, a +ve definite matrix stored in lower-triangular form.
!    n       = input, the order of a
!    nn      = input, the size of the a and u arrays >= n*(n+1)/2
!    u()     = output, a lower triangular matrix such that u*u' = a.
!              a & u may occupy the same locations.
!    nullty  = output, the rank deficiency of a.
!    ifault  = output, error indicator
!                    = 1 if n < 1
!                    = 2 if a is not +ve semi-definite
!                    = 3 if nn < n*(n+1)/2
!                    = 0 otherwise

!***********************************************************************

IMPLICIT NONE
INTEGER, PARAMETER     :: dp = SELECTED_REAL_KIND(15, 60)
REAL (dp), INTENT(IN)  :: a(:)
INTEGER, INTENT(IN)    :: n, nn
REAL (dp), INTENT(OUT) :: u(:)
INTEGER, INTENT(OUT)   :: nullty, ifault

! Local variables
REAL (dp)              :: eta2, x, w
INTEGER                :: i, icol, ii, irow, j, k, kk, l, m

!       The value of eta will depend on the word-length of the
!       computer being used.  See introductory text.

REAL (dp), PARAMETER   :: eta = 1.E-09_dp, zero = 0.0_dp

ifault = 1
IF (n <=  0) RETURN
ifault = 3
IF (nn < n*(n+1)/2) RETURN
ifault = 2
nullty = 0
j = 1
k = 0
eta2 = eta*eta
ii = 0

!       Factorize column by column, icol = column no.

DO icol = 1, n
  ii = ii + icol
  x = eta2*a(ii)
  l = 0
  kk = 0
  
!       IROW = row number within column ICOL
  
  DO irow = 1, icol
    kk = kk + irow
    k = k + 1
    w = a(k)
    m = j
    DO i = 1, irow
      l = l + 1
      IF (i == irow) EXIT
      w = w - u(l)*u(m)
      m = m + 1
    END DO

    IF (irow == icol) EXIT
    IF (u(l) /= zero) THEN
      u(k) = w / u(l)
    ELSE
      IF (w*w > ABS(x*a(kk))) RETURN
      u(k) = zero
    END IF
  END DO

  IF (ABS(w) > ABS(eta*a(k))) THEN
    IF (w < zero) RETURN
    u(k) = SQRT(w)
  ELSE
    u(k) = zero
    nullty = nullty + 1
  END IF
  j = j + icol
END DO

ifault = 0
RETURN
END SUBROUTINE chol




SUBROUTINE subchl (a, b, n, u, nullty, ifault, det)

!   REMARK ASR 44  APPL. STATIST. (1982) VOL. 31, NO. 3

!   A revised and enhanced version of
!   ALGORITHM AS 6  APPL. STATIST. (1968) VOL. 17, NO. 2

! N.B. Argument NDIM has been removed, and DET is an OPTIONAL argument

!   Given a symmetric matrix of order N as lower triangle in A(),
!   calculates an upper triangle, U(), such that U'U = the sub-matrix
!   of A whose rows and columns are specified in the integer array B().
!   U() may coincide with A().   A() must be +ve semi-definite.
!   ETA is set to multiplying factor determining effective zero for a pivot.
!   NULLTY is returned as number of effective zero pivots.
!   IFAULT is returned as 1 if N <= 0, 2 if A() is not +ve semi-
!   definite, otherwise 0 is returned.

IMPLICIT NONE
INTEGER, PARAMETER               :: dp = SELECTED_REAL_KIND(15, 60)
REAL (dp), INTENT(IN)            :: a(:)
REAL (dp), INTENT(OUT)           :: u(:)
REAL (dp), INTENT(OUT), OPTIONAL :: det
INTEGER, INTENT(IN)              :: b(:), n
INTEGER, INTENT(OUT)             :: nullty, ifault

!     Local variables

REAL (dp) :: w, eta2, x
INTEGER   :: i, icol, ii, ij, irow, j, jj, k, kk, l, m

!     The value of ETA below will depend upon the word length of the
!     computer being used.

REAL (dp), PARAMETER :: eta = 1.E-14_dp, zero = 0.0_dp, one = 1.0_dp

ifault = 1
IF (n <= 0) GO TO 90
ifault = 2
nullty = 0
IF (PRESENT(det)) det = one
j = 1
k = 0
eta2 = eta*eta
DO icol = 1, n
  ij = b(icol)*(b(icol) - 1)/2
  ii = ij + b(icol)
  x = eta2*a(ii)
  l = 0
  DO irow = 1, icol
    kk = b(irow)*(b(irow) + 1)/2
    k = k + 1
    jj = ij + b(irow)
    w = a(jj)
    m = j
    DO i = 1, irow
      l = l + 1
      IF (i == irow) GO TO 20
      w = w - u(l)*u(m)
      m = m + 1
    END DO
    20 IF (irow == icol) GO TO 50
    IF (u(l) == zero) GO TO 30
    u(k) = w/u(l)
    CYCLE
    30 IF (w*w > ABS(x*a(kk))) GO TO 90
    u(k) = zero
  END DO

  50 IF (ABS(w) <= ABS(eta*a(kk))) GO TO 60
  IF (w < zero) GO TO 90
  u(k) = SQRT(w)
  GO TO 70
  60 u(k) = zero
  nullty = nullty + 1
  70 j = j + icol
  IF (PRESENT(det)) det = det*u(k)*u(k)
END DO

ifault = 0
90 RETURN
END SUBROUTINE subchl