SUBROUTINE GAUI(NP,AMU,H,LDH,BETA)
C  This routine generates a random vector from a Gaussian distribution with
C  mean vector mu and precision matrix H, subject to the restriction that the
C  vector correspond to an invertible lag operator.
C  Inputs:
C  NP        Dimension of vector
C  AMU       Mean vector
C  H         Precision matrix H
C  LDH       Leading dimension of H
C  Output:
C  BETA      Random vecctor
      IMPLICIT REAL*8 (A-H,O-Z)
      LOGICAL LFLIP,LINOP
      PARAMETER(LD=100)
      DIMENSION AMU(NP),H(LDH,NP),BETA(NP)
      DIMENSION V1(LD),V2(LD)
      DATA NREP1,NREP2/250,250/
C  Make independent draws from normal with rejection.
      CALL GAU(NP,AMU,H,LDH,V1,V2)
      IF(LINOP(NP,V1))GO TO 20
      DO 10 IREP=1,NREP1
      CALL GAUA(NP,AMU,H,LDH,V1,V2)
   10 IF(LINOP(NP,V1))GO TO 20
      GO TO 100
   20 CALL DCOPY(NP,V1,1,BETA,1)
      RETURN
C  Make draws using independence Metropolis with uniform p.a.c.f. source.
C  Separate entry can be used to force independence Metropolis algorithm
C  (mostly for testing purposes).
      ENTRY GAUI1(NP,AMU,H,LDH,BETA)
  100 WTLOG0=DMACH(8)
      DO 120 IREP=1,NREP2
      DO 110 IP=1,NP
  110 V1(IP)=-1.0D0+2.0D0*DRNUNF()
      CALL PACFAR(NP,V1,V2)
      WTLOG=GAUK(NP,AMU,H,LDH,V2)+ARPALJ(NP,V1)
      IF(.NOT.LFLIP(WTLOG-WTLOG0))GO TO 120
      CALL DCOPY(NP,V2,1,BETA,1)
      WTLOG0=WTLOG
  120 CONTINUE
      RETURN
C  Auxilliary entry to provide a weighted sample (mostly for testing purposes)
C  Additional output:
C  WTL        Log weight, for importance sampling
      ENTRY GAUI2(NP,AMU,H,LDH,BETA,WTL)
      DO 210 IP=1,NP
  210 V1(IP)=-1.0D0+2.0D0*DRNUNF()
      CALL PACFAR(NP,V1,BETA)
      WTL=GAUK(NP,AMU,H,LDH,BETA)+ARPALJ(NP,V1)
      RETURN
      END 
C
      REAL*8 FUNCTION GAUIN(NP,AMU,H,LDH,SD,TIME)
C  This function evaluates the log normalizing constant for the log density
C  kernel of a random vector computed in GAUK, given an invertibility 
C  restriction.  This normalizing constant, when multiplied by the kernel
C  in GAUK, yields the p.d.f. of the random vector.
C  Inputs:
C  NP        Dimension of vector
C  AMU       Mean vector
C  H         Precision matrix H
C  LDH       Leading dimension of H
C  SD        Routine stops when standard deviation of approximation error
C            of log kernel reaches this point.
C  TIME      Routine stops when this much CPU time has been used.
C  Output:
C  GAUIN     Log normalizing constanat (approximate
C  SD        Standard error of approximation
C  TIME      CPU time used
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*4 CPSEC
      LOGICAL LINOP
      PARAMETER(LD=100)
      DIMENSION AMU(NP),H(LDH,NP)
      DIMENSION V1(LD),V2(LD)
      DATA MREP1,MREP2,JUMP1,JUMP2/500,1000,500,1000/
C  Measure efficiency of sampling from untruncated prior.
      CLOCK=CPSEC()
      AA1=0.0D0
      SS1=0.0D0
      DO 10 IREP=1,MREP1
      CALL GAU(NP,AMU,H,LDH,V1,V2)
      BB=0.0D0
      IF(LINOP(NP,V1))BB=BB+0.5D0
      IF(LINOP(NP,V2))BB=BB+0.5D0
      AA1=AA1+BB
   10 SS1=SS1+BB**2
      P=AA1/MREP1
      IF(P.EQ.0.0D0)GO TO 200
      IF(P.EQ.1.0D0)GO TO 300
      VARLOG=(SS1-MREP1*P**2)/(MREP1*(MREP1-1)*P**2)
      IF((CPSEC()-CLOCK).EQ.0)GO TO 100
      EFF1=1.0D0/(VARLOG*(CPSEC()-CLOCK))
C  Measure efficiency of sampling uniformly from p.a.c.f.
      CLOCK2=CPSEC()
      AA2=0.0D0
      SS2=0.0D0
      DO 30 IREP=1,MREP2
      DO 20 IP=1,NP
   20 V1(IP)=2.0D0*DRNUNF()-1.0D0
      CALL PACFAR(NP,V1,V2)
      BB=DEXP(ARPALJ(NP,V1)+GAUK(NP,AMU,H,LDH,V2))
      AA2=AA2+BB
   30 SS2=SS2+BB**2
      C=(2.0D0**NP)*DEXP(GAUN(NP,H,LDH))
      P=C*AA2/MREP2
      VARLOG=(SS2*C**2-MREP2*P**2)/(MREP2*(MREP2-1)*P**2)
      EFF2=1.0D0/(VARLOG*(CPSEC()-CLOCK2))
      IF(EFF2.GT.EFF1)GO TO 200
C  Approximate by sampling from untruncated prior.
  100 NREP1=MREP1
  110 DO 120 JUMP=1,JUMP1
      CALL GAU(NP,AMU,H,LDH,V1,V2)
      BB=0.0D0
      IF(LINOP(NP,V1))BB=BB+0.5D0
      IF(LINOP(NP,V2))BB=BB+0.5D0
      AA1=AA1+BB
  120 SS1=SS1+BB**2
      NREP1=NREP1+JUMP1
      P=AA1/NREP1
      SDLOG=DSQRT((SS1-NREP1*P**2)/(NREP1*(NREP1-1)*P**2))
      IF(SDLOG.LT.SD)GO TO 300
      IF(CPSEC()-CLOCK.GT.TIME)GO TO 300
      GO TO 110
C  Approximate by sampling uniformly from p.a.c.f.
  200 NREP2=MREP2
  210 DO 230 JUMP=1,JUMP2
      DO 220 IP=1,NP
  220 V1(IP)=2.0D0*DRNUNF()-1.0D0
      CALL PACFAR(NP,V1,V2)
      BB=DEXP(ARPALJ(NP,V1)+GAUK(NP,AMU,H,LDH,V2))
      AA2=AA2+BB
  230 SS2=SS2+BB**2
      NREP2=NREP2+JUMP2
      P=C*AA2/NREP2
      SDLOG=DSQRT((SS2*C**2-NREP2*P**2)/(NREP2*(NREP2-1)*P**2))
      IF(SDLOG.LT.SD)GO TO 300
      IF((CPSEC()-CLOCK).GT.TIME)GO TO 300
      GO TO 210
C  Return parameters and value.
  300 SD=SDLOG
      TIME=CPSEC()-CLOCK
      GAUIN=GAUN(NP,H,LDH)-DLOG(P)
      RETURN
      END