SUBROUTINE WISH(M,SIG,LDSIG,ANU,A,LDA)
C  This routine generates a pseudo-random, mxm matrix A from a Wishart
C  distribution with matrix parameter Sigma and degrees of freedom
C  parameter Nu. 
C  Inputs:
C  M       Order of matrix, m
C  SIG     Matrix parameter, Sigma
C  LDSIG   Leading dimension of SIG
C  ANU     Degrees of freedom parameter, Nu
C  LDA     Leading dimension of A
C  Output:
C  A       Random matrix A
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(LD=100)
      COMMON/SCRA/V1(LD),V2(LD),A1(LD,LD),A2(LD,LD)
      DIMENSION SIG(LDSIG,M),A(LDA,M),R(LDR,M)
      CALL ULFTDS(M,SIG,LDSIG,A1,LD)
C  Entry WISHA: Matrix parameter already factored
C  *** Can be no interveneing calls that change A1 !! ***
      ENTRY WISHA(M,SIG,LDSIG,ANU,A,LDA)
    5 DF=ANU
      DO 20 I=1,M
      DO 10 J=1,I
      A(J,I)=0.0D0
   10 A(I,J)=DRNNOF()
      CALL DRNCHI(1,DF,AA)
      A(I,I)=DSQRT(AA)
   20 DF=DF-1.0D0
      CALL DMXTYF(M,M,A1,LD,M,M,A,LDA,M,M,A2,LD)
      CALL DMXYTF(M,M,A2,LD,M,M,A2,LD,M,M,A,LDA)
      RETURN
C  Entry WISHB:  Coming in with matrix parameter already factored
C  R: R'R=Matrix parameter, R upper triangular with zeros set below
      ENTRY WISHB(M,R,LDR,ANU,A,LDA)
      CALL UCRGRG(M,M,R,LDR,A1,LD)
      GO TO 5
      END
      REAL*8 FUNCTION WISHK(M,P,LDP,ANU,A,LDA)
C  This function evaluates the log density kernel of a Wishart random
C  matrix A with matrix parameter Sigma and degrees of freedom parameter
C  Nu.  The log density kernel is .5(Nu+1-m)logdet(A)-.5tr[(Sigma^-1)A].
C  Inputs:
C  M       Order of matrix, m
C  P       Square matrix P: PP'=Sigma^-1
C  LDP     Leading dimension of P
C  ANU     Degrees of freedom parameter, Nu
C  A       Matrix point of evaluation of log density kernel
C  LDA     Leading dimension of A
C  Output:
C  WISHK   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)
      CALL DLFTDS(M,A,LDA,A1,LD)
      CALL UXTAXF(M,M,P,LDP,A,LDA,A2,LD)
      AA=0.0D0
      BB=0.0D0
      DO 10 I=1,M
      AA=AA+DLOG(A1(I,I))
   10 BB=BB+A2(I,I)
      WISHK=(ANU-1-M)*AA-0.5D0*BB
      RETURN
      END
      REAL*8 FUNCTION WISHN(M,SIG,LDSIG,ANU,P,LDP)
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 matrix computed in WISHK.  This 
C  normalizing constant, when multiplied by the kernel in WISHK,
C  yields the p.d.f. of the random matrix.
C  Inputs:
C  M       Order of matrix, m
C  SIG     Matrix parameter, Sigma
C  LDSIG   Leading dimension of SIG
C  ANU     Degrees of freedom parameter, Nu
C  LDP     Leading dimension of P
C  Outputs:
C  WISHN   Log normalizing factor
C  P       Upper triangular mxm matrix: PP'=Sigma^-1
      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)
      CALL DLFTDS(M,SIG,LDSIG,A1,LD)
      AA=0.0D0
      BB=0.0D0
      ANU1=ANU+1.0D0
      DO 10 I=1,M
      AA=AA+DLOG(A1(I,I))
   10 BB=BB+DLNGAM(0.5D0*(ANU1-I))
      WISHN=-0.5D0*M*(ANU*DLOG2+0.5D0*(M-1)*DLGPI)-ANU*AA-BB
      CALL DLINRT(M,A1,LD,2,P,LDP)
      RETURN
      END