random_mvnorm Subroutine

SUBROUTINE random_mvnorm(n, h, d, f, first, x, ier)
  ! Adapted from Fortran 77 code from the book:
  !     Dagpunar, J. 'Principles of random variate generation'
  !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  ! N.B. An extra argument, ier, has been added to Dagpunar's routine
  !     SUBROUTINE GENERATES AN N VARIATE RANDOM NORMAL
  !     VECTOR USING A CHOLESKY DECOMPOSITION.
  ! ARGUMENTS:
  !        N = NUMBER OF VARIATES IN VECTOR
  !           (INPUT,INTEGER >= 1)
  !     H(J) = J'TH ELEMENT OF VECTOR OF MEANS
  !           (INPUT,REAL)
  !     X(J) = J'TH ELEMENT OF DELIVERED VECTOR
  !           (OUTPUT,REAL)
  !
  !    D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I)
  !            (INPUT,REAL)
  !    F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR
  !           DECOMPOSITION OF VARIANCE MATRIX (J <= I)
  !            (OUTPUT,REAL)
  !    FIRST = .TRUE. IF THIS IS THE FIRST CALL OF THE ROUTINE
  !    OR IF THE DISTRIBUTION HAS CHANGED SINCE THE LAST CALL OF THE ROUTINE.
  !    OTHERWISE SET TO .FALSE.
  !            (INPUT,LOGICAL)
  !    ier = 1 if the input covariance matrix is not +ve definite
  !        = 0 otherwise
  INTEGER, INTENT(IN)   :: n
  REAL, INTENT(IN)      :: h(:), d(:)   ! d(n*(n+1)/2)
  REAL, INTENT(IN OUT)  :: f(:)         ! f(n*(n+1)/2)
  REAL, INTENT(OUT)     :: x(:)
  LOGICAL, INTENT(IN)   :: first
  INTEGER, INTENT(OUT)  :: ier
  !     Local variables
  INTEGER       :: j, i, m
  REAL          :: y, v
  INTEGER, SAVE :: n2
  IF (n < 1) THEN
     WRITE(*, *) 'SIZE OF VECTOR IS NON POSITIVE'
     STOP
  END IF
  ier = 0
  IF (first) THEN                        ! Initialization, if necessary
     n2 = 2*n
     IF (d(1) < zero) THEN
        ier = 1
        RETURN
     END IF
     f(1) = SQRT(d(1))
     y = one/f(1)
     DO j = 2,n
        f(j) = d(1+j*(j-1)/2) * y
     END DO
     DO i = 2,n
        v = d(i*(i-1)/2+i)
        DO m = 1,i-1
           v = v - f((m-1)*(n2-m)/2+i)**2
        END DO
        IF (v < zero) THEN
           ier = 1
           RETURN
        END IF
        v = SQRT(v)
        y = one/v
        f((i-1)*(n2-i)/2+i) = v
        DO j = i+1,n
           v = d(j*(j-1)/2+i)
           DO m = 1,i-1
              v = v - f((m-1)*(n2-m)/2+i)*f((m-1)*(n2-m)/2 + j)
           END DO ! m = 1,i-1
           f((i-1)*(n2-i)/2 + j) = v*y
        END DO ! j = i+1,n
     END DO ! i = 2,n
  END IF
  x(1:n) = h(1:n)
  DO j = 1,n
     y = random_normal()
     DO i = j,n
        x(i) = x(i) + f((j-1)*(n2-j)/2 + i) * y
     END DO ! i = j,n
  END DO ! j = 1,n
  RETURN
END SUBROUTINE random_mvnorm