SUBROUTINE GAU(K,AMU,H,LDH,X1,X2)
C  This routine generates two antithetic pseudo-random vectors x1 and x2 
C  from a Gaussian distribution with mean vector mu and precision 
C  matrix H.
C  Inputs:
C  K       Dimension of x
C  AMU     Mean vector mu (kx1)
C  H       Precision matrix H (kxk)
C  LDH     Leading dimension of H
C  Outputs:
C  X1,X2   Random draws
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=100)
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      DIMENSION AMU(K),H(LDH,K),X1(K),X2(K)
C  Entry GAU:  Factor precision
      CALL DLFTDS(K,H,LDH,A1,LD)
      CALL DLINRT(K,A1,LD,2,A1,LD)
C  Entry GAUA: Precision already factored;
C  *** Can be no intervening calls that change A1 !! ***
      ENTRY GAUA(K,AMU,H,LDH,X1,X2)
   10 CALL URNMVN(K,A1,LD,V1)
      CALL UVADD(K,AMU,V1,X1)
      CALL UVSUB(K,AMU,V1,X2)
      RETURN
C  Entry GAUB:  Coming in with R: RR'=Variance, R upper triangular
      ENTRY GAUB(K,AMU,R,LDR,X1,X2)
      CALL UCRGRG(K,K,R,LDR,A1,LD)
      GO TO 10
      END
      REAL*8 FUNCTION GAUK(K,AMU,H,LDH,X)
C  This function evaluates the log density kernel of a multivariate
C  normal random vector x with mean mu and precision matrix H.  The
C  kernel density is exp(-.5(x-mu)'H(x-mu))
C  Inputs:
C  K       Dimension of x
C  AMU     Mean vector mu (kx1)
C  H       Precision matrix h (kxk)
C  LDH     Leading dimension of H
C  X       Point of evaluation, x (kx1)
C  Output:
C  GAUK    Evaluation of log kernel density
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=100)
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      DIMENSION H(LDH,K),X(K)
      CALL UVSUB(K,X,AMU,V1)
      GAUK=-0.5D0*UQUAF(K,V1,H,LDH)
      RETURN
      END
      REAL*8 FUNCTION GAUN(K,H,LDH)
C***********************************************************************
C  In the calling routine a call to INIT must be made prior to calling *
C  this routine for the first time.                                    *
C***********************************************************************    
C  This routine evaluates the log normalizing constant for the log 
C  density kernel of a random vector computed in GAUK.  This 
C  normalizing constant, when multiplied by the kernel in GAUK,
C  yields the p.d.f. of the random vector.
C  Inputs:
C  K       Dimension of random vector
C  H       Precision matrix H (kxk)
C  LDH     Leading dimension of H
C  Output:
C  GAUN    Log normalizing factor
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=100)
      COMMON/CONST/DLOG2,DLGPI,HLG2PI,PI,PII
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      DIMENSION H(LDH,K)
      CALL DLFTDS(K,H,LDH,A1,LD)
      AA=0.0D0
      DO 10 I=1,K
   10 AA=AA+DLOG(A1(I,I))
      GAUN=-K*HLG2PI+AA
      RETURN
      END