C This file contains Fortran subroutines called by the Splus function "rcwish".
C It should be saved as "rcwish.f", and compiled as a separate subroutine by
C typing "f77 -c rcwish.f" at the unix prompt (call may vary with systems). 
C  
C
C rscwish:      Returns N Standard Constrained Wishart matrices with
C ^^^^^^^       df degrees of freedom, and constrained by the vector d.
C               (see Everson & Morris (2000) JCGS for details).
C
C augment:      Augments a (t-1)x(t-1) acceptable candidate matrix (Tmat)
C ^^^^^^^       into a txt candidate matrix (possibly acceptable).
C
C checkcon:     Checks a condition on a txt candidate matrix to see if
C ^^^^^^^^      it is acceptable. If it is, an updated inverse is computed.
C
C load:		Loads an acceptable CWish matrix into the output array Ua. 
C ^^^^		
C
C mmult:        Matrix multiplication.
C ^^^^^
C
C mammult:      Rescales an array by pre- and post-multiplying each
C ^^^^^^^        component by a matrix.
C	
C rcchisq:      Generates a random Constrained Chi-square(df;d) variate.
C ^^^^^^^
C 	
C rchisq:       Generates a random (unconstrained) Chi-square(df;d) variate.
C ^^^^^^
C
C pchisq:       Evaluates the Chi-square(df) CDF.
C ^^^^^^
C
C qchisq:       Inverts the Chi-square(df) CDF to return quantiles.
C ^^^^^^
C
C gqchisq:      Same as qchisq, but with less expensive inputs.
C ^^^^^^^
C
C bisect:       Bisection routine for inverting the Chi-square(df) CDF. 
C ^^^^^^
C
C gammln:	Evaluates the log-gamma function.
C ^^^^^^ Taken from Press et. al. (1992) Numerical Recipes in Fortran, 2nd ed..
C  
C gammp:        Returns the incomplete gamma function P(a,x).
C ^^^^^  Taken from Press et. al. (1992) Numerical Recipes in Fortran, 2nd ed..
C
C gser:         Series representation of  P(a,x).
C ^^^^   Taken from Press et. al. (1992) Numerical Recipes in Fortran, 2nd ed..
C
C gcf:          Continued fraction representation of  P(a,x).
C ^^^    Taken from Press et. al. (1992) Numerical Recipes in Fortran, 2nd ed..
C
C rnorm:        Returns two random N(0,1) variates.
C ^^^^^
C
C ran2:         Returns a random U(0,1) variate.      
C ^^^^   Taken from Press et. al. (1992) Numerical Recipes in Fortran, 2nd ed..
C
C
C
C 
      SUBROUTINE rscwish(Ua, N, p, df, d, pd, chitab, seed, Tmat, 
     $		Deltinv, H11, H12, Zt, Zts, tempp, mcnt, drvcnt, 
     $		orvcnt, rejcnt)
C Returns N pxp Constrained Wishart(df, I; diag(d)) matrices in Ua. 
C		(see Everson & Morris (2000) JCGS for details).
      INTEGER N, p, seed, step, i, j, cnt, mcnt, drvcnt, orvcnt, accept
      INTEGER rejcnt(p)
      DOUBLE PRECISION Ua(p,p,N), df, d(p), pd(p), chitab(999)
      DOUBLE PRECISION Tmat(p,p), Ut, Zt(p), Zts(p), dt, pdt, Xts 
      DOUBLE PRECISION Uts, Deltinv(p,p), H11(p,p), H12(p), tempp(p)
      DOUBLE PRECISION lambda, lamsv(p), Usv(p)
      lambda=0.
      mcnt=0
      drvcnt=0
      orvcnt=0
C Keep track of how many matrices are begun (=mcnt) and how many diagonal 
C random variables (=drvcnt) and off-diagonal random variables (=orvcnt)
C are generated.
      cnt=0
C Keep track of how many acceptable matrices have been generated (=cnt).
      do while(cnt.lt.N)
       mcnt = mcnt + 1
       step=1
       pdt=pd(step)
       drvcnt = drvcnt + 1
       call rcchisq(Ut, df, pdt, chitab, seed)
       Tmat(1,1)=dsqrt(Ut)
       accept=1
       Deltinv(1,1) = 1./(d(1)-Ut)
       do while(step.lt.p.and.accept.eq.1)
        drvcnt = drvcnt + 1
        orvcnt = orvcnt + step
        step=step+1
        pdt=pd(step)
        call augment(Tmat, step, p, df, pdt, chitab, Zt, Zts, Xts, 
     $		Ut, Uts, seed)
        dt=d(step)
        if(dt.gt.0) then
         call checkcon(Tmat, step, p, dt, Zt, Ut, accept, Deltinv,
     $		H11, H12, lambda, tempp)
        endif 
       end do
       if(accept.eq.1) then
C An acceptable CWish matrix has been generated. Update cnt and load 
C U = Tmat%*%t(Tmat) into the next open position in Ua (pxpxN).
        cnt = cnt + 1
        call load(Ua, Tmat, p, N, cnt)
       else
	rejcnt(step)=rejcnt(step)+1
       endif
      end do
      return
      end




      SUBROUTINE augment(Tmat, step, p, df, pdt, chitab, Zt, Zts,
     $			Xts, Ut, Uts, seed)
C Augments a (t-1)x(t-1) Tmat to a txt (step = t), by simulating a  
C Constrained Chi-square(df; dt) and partitioning it with Betas to fill  
C in the t-th row of Tmat (t=2,..,p). 
C
C p:		physical dimension of Tmat
C pd:		P(U < dt), for U ~ Chi-square(df)
C chitab:	table of 999 evenly-spaced Chi-square(df) quantiles
C seed:		seed for random number generator
C
      INTEGER step, seed, p, i, j, k
      DOUBLE PRECISION Tmat(p,p), df, pdt, chitab(999), Zt(p), Zts(p)
      DOUBLE PRECISION Xts, Ut, Uts, z1, z2, ZZ, u1, u2, pi, ran2
      pi=3.14159265
      call rcchisq(Ut, df, pdt, chitab, seed)
      k=int((step-1)/2.)
C k is the largest integer less than or equal to (t-1)/2 (step=t)
      ZZ=0.0
      do 100 i=1,k
       call rnorm(z1, z2, seed)
       Zts(2*i-1)=z1
       Zts(2*i)=z2
       ZZ = ZZ + z1**2 + z2**2
 100  continue
      if(step-1.gt.2*k) then
C if t-1 is odd, get one more N(0,1) to fill out Zts:
       u1=ran2(seed)
       u2=ran2(seed)
       z1 = dcos(2.*pi*u2)*dsqrt((-2.)*dlog(u1))
       Zts(step-1)=z1
       ZZ = ZZ + z1**2
      endif
C Call the bisection routine to generate a Chi-square(df-step+1) rv Xts.
      call rchisq(Xts, df-step+1, seed)
C Zts[1:(t-1)] ~ iid N(0,1); Xts ~ Chi-square(df-t+1)
      Uts = ZZ + Xts
      do 300 j=1,step-1
       Zt(j) = Zts(j)*dsqrt(Ut/Uts)
       Tmat(step,j) = Zt(j)
 300  continue
      Tmat(step,step)=dsqrt(Xts*Ut/Uts)
      return
      end







      SUBROUTINE checkcon(Tmat,step,p,dt,Zt,Ut,accept,Deltinv,
     $ H11,H12,lambda, tempp)
C Checks to see if the t-th step (step=t) of the construction of Tmat  
C is acceptable. If not acceptable (i.e., if det(D-TT') is negative 
C after step t) then "accept" is returned as 0. Otherwise, "accept" 
C is returned as 1, and "Deltinv", which is input as the (t-1)x(t-1) 
C inverse of D-TT', is updated to contain the txt inverse of D-TT'.
C
      INTEGER step, p, i, j, k, accept
      DOUBLE PRECISION Tmat(p,p), dt, Zt(p), Ut, Deltinv(p,p) 
      DOUBLE PRECISION H11(p,p), H12(p), H22, lambda
      DOUBLE PRECISION tempp(p),  x
C
      call mmult(Tmat, Zt, p, p, 1, tempp, step-1, step-1, 1) 
      call mmult(Deltinv, tempp, p,p,1, H12, step-1, step-1, 1) 
      call mmult(tempp, H12, 1, p, 1, x, 1, step-1, 1)
      lambda = dt - Ut - x
C
      if(lambda.lt.0.) then
       accept = 0
C lambda < 0 at any step t means we have to start over with t=1.
      else 
C Update Deltinv from (t-1)x(t-1) to txt. Leave accept=1.
       if(step.lt.p) then
C Does not update Deltinv if t=p (step=t=p). That would indicate 
C that the matrix is finished and acceptable.
        H22 = 1./lambda
        do 200 i=1,step-1
         H12(i) = H12(i)/lambda
         do 100 j=1,step-1
          H11(i,j) = Deltinv(i,j) + H12(i)*H12(j)/H22
 100     continue
 200    continue
        do 400 i=1,step-1
         Deltinv(step,i) = H12(i)
         Deltinv(i,step) = H12(i)
         do 300 j=1,step-1
          Deltinv(i,j)=H11(i,j)
	  Deltinv(j,i)=H11(i,j)
 300     continue
 400    continue
        Deltinv(step,step)=H22
       endif
      endif
C
      return
      end

      SUBROUTINE load(Ua, Tmat, p, N, cnt)
C Computes Tmat%*%t(Tmat) and loads into Ua(,,cnt)
      INTEGER i, j, k, cnt, p, N
      DOUBLE PRECISION Ua(p,p,N), Tmat(p,p)
C
      Ua(1,1,cnt)=Tmat(1,1)**2
      do 300 i=2,p
       Ua(i,i,cnt)=Tmat(i,i)**2
       do 200 j = 1,i-1
        Ua(i,i,cnt) = Ua(i,i,cnt) + Tmat(i,j)**2
        Ua(i,j,cnt) = 0.0
        do 100 k=1,j
         Ua(i,j,cnt)=Ua(i,j,cnt)+Tmat(i,k)*Tmat(j,k)
 100    continue
        Ua(j,i,cnt)=Ua(i,j,cnt)
 200   continue
 300  continue
      return
      end

       
      SUBROUTINE mmult(A,B,p,q,r,C, pp,qq,rr)
C
C Returns the matrix product A%*%B = C. p,q,r are physical dimensions;
C pp,qq,rr are corresponding working dimensions.
      INTEGER p,q,r,pp,qq,rr,i,j,k
      DOUBLE PRECISION t, A(p,q), B(q,r), C(p,r)
      do 100 i=1,pp
       do 200 j=1,rr
        t=0.0
        do 300 k=1,qq
         t=t+A(i,k)*B(k,j)
 300    continue
        C(i,j)=t
 200   continue
 100  continue
      return
      end


      SUBROUTINE mammult(M, tM, A, p, N, U, V, MAM)
C
C Each of N pxp elements of A (pxpxN) is pre-multiplied by M and 
C post-multiplied by tM each pxp element of a pxpxn array (A)
C by a pxp matrix (M) and its transpose (Mt), respectively. 
C 
      INTEGER p, N, i, j, ind
      DOUBLE PRECISION M(p,p), tM(p,p), A(p,p,N)
      DOUBLE PRECISION MAM(p,p,N), U(p,p), V(p,p)
C
      do 500 ind=1,N
       do 200 i=1,p
        do 100 j=1,p
         U(i,j)=A(i,j,ind)
 100	 continue
 200	continue
	call mmult(M,U,p,p,p,V,p,p,p)
	call mmult(V,tM,p,p,p,U,p,p,p)
	do 400 i=1,p
	 do 300 j=1,p
	  MAM(i,j,ind)=U(i,j)
 300	 continue
 400    continue
 500   continue
       return
       end


      SUBROUTINE rcchisq(x, df, pd, chitab, seed)
C Returns one draw from the constrained Chi-square distribution:
C a Chi-square(df) draw that is constrained to be less than d,
C where pd is the probability that a Chi-square(df) is < d.
C "chitab" is a table of Chi-square(df) quantiles.
C
      INTEGER seed
      DOUBLE PRECISION x, df, pd, chitab(999),  p, u, gammp, ran2
      p = pd*ran2(seed)
      call qchisq(p, df, chitab, x)
      return
      end

      SUBROUTINE rchisq(x, df,  seed)
C Returns one draws from the (unconstrained) Chi-sqare(df) distribution.
C Calls gqchisq, and so does not require a table (chitab).
C
      INTEGER i, k, seed
      DOUBLE PRECISION x, df, p,  z1,z2,ran2
      if(df.gt.20.) then
C For larger df use the inverse CDF method.
       p = ran2(seed)
       call gqchisq(p,df,x)
      else
C For smaller df, add up squared N(0,1)'s to generate x.
       k = int(df/2.)
C k is the largest integer leq df/2
       x=0.0
       do 200 i=1,k
        call rnorm(z1, z2, seed)
        x = x + z1**2 + z2**2
 200   continue
       if(df.gt.k*2.) then
C if df is odd, get one more N(0,1) to finish x:
        call rnorm(z1, z2, seed)
        x = x + z1**2
       endif
      endif
      return
      end


      SUBROUTINE pchisq(xx,df,p)
C Call gammp to return p, the CDF of the Chi-square(df) distribution at x.
      DOUBLE PRECISION df, gammp, x, xx, p, a 
      a = df/2.
      x = xx/2.
      p = gammp(a, x)
      return
      end



      SUBROUTINE qchisq(p, df, chitab, q)
C Returns the pth quantile of the Chi-square(df) distribution as q.
C chitab is a table of equally spaced Chi-square(df) quantiles.
      INTEGER  ITMAX, cp, i
      DOUBLE PRECISION p, df, chitab(999), q, EPS, l, u, fl, fu, sd,gammp
      PARAMETER (ITMAX=1000, EPS=3.e-7)
      if(p.gt.0.9990) then
       l = chitab(999)
       fl= 0.9990
       sd = dsqrt(df*2.)
       u = l+sd
       call pchisq(u, df, fu)
       do while(fu.lt.p)
        l = u
        fl = fu
        u = u+sd
        call pchisq(u, df, fu)
       end do
      else
       if(p.lt.0.001) then
        l = EPS
        fl = EPS
        u=chitab(1)
        fu=0.0010
       else
	cp = int(1000.*p)
        fl = cp/1000.
	l = chitab(cp)
	fu = (cp+1)/1000.
	u = chitab(cp+1)
       endif
      endif
C now l and u bracket the desired quantile q, with l < q < u 
C  and fl < p < fu.
C
      q=(l+u)/2.
      i=0
      do while(u-l.gt.EPS*q.and.i.lt.ITMAX)
       i = i+1
       call bisect(p, df, l, q, u,fl, fu)
      end do
      return
      end


      SUBROUTINE gqchisq(p, df, q)
C Returns the pth quantile of the Chi-square(df) distribution as q.
C g = "general" version - does not require a quantile table. Intended
C for use with larger degrees of freedom (df > 40).
C
      INTEGER  ITMAX, cp, i
      DOUBLE PRECISION p, df, chitab(999), q, EPS, l, u, fl, fu, sd,gammp
      PARAMETER (ITMAX=1000, EPS=3.e-7)
      if(p.lt.0.5) then
       l=0.
       fl=0.
       u=df
       call pchisq(u, df, fu)
      else
       sd=dsqrt(df*2.)
       l=max(0.,df-sd)
       call pchisq(l,df,fl)
       do while(fl.gt.p)
	u=l
	fu=fl
	l=max(0.,l-sd)
	call pchisq(l,df,fl)
       end do
       u=df+sd
       call pchisq(u,df,fu)
       do while(fu.lt.p)
	l=u
	fl=fu
	u=u+sd
	call pchisq(u,df,fu)
       end do
      endif
C now l and u bracket the desired quantile q, with l < q < u 
C  and fl < p < fu.
C
      q=(l+u)/2.
      i=0
      do while(u-l.gt.EPS*q.and.i.lt.ITMAX)
       i = i+1
       call bisect(p, df, l, q, u,fl, fu)
      end do        
      return
      end


      SUBROUTINE bisect(p,df, l, mid, u, fl, fu)
C Bisects a lower and upper bound on the Chi-square(df) CDF
C to return a narrower bound.
      DOUBLE PRECISION df, p, a, x, l, u, fl, fu, mid, fm, gammp
      a = df/2.
      mid = (l+u)/2.
      x = mid/2.
      fm = gammp(a, x) 
      if(fm.lt.p) then
       l = mid
       fl = fm
      else
       u = mid
       fu = fm
      endif   
      return
      end




      DOUBLE PRECISION FUNCTION gammln(xx)
C Returns the value of Lanczos' approximation to ln(gamma(xx)), for n > 0.
C Taken from Numerical Recipes in Fortran, 2nd ed., Cambridge University Press.
C
      INTEGER j
      DOUBLE PRECISION xx, ser, stp, tmp, x, y, cof(6)
      SAVE cof, stp
      DATA cof, stp/76.18009172947146d0, -86.50532032941677d0,
     $ 24.01409824083091d0, -1.231739572450155d0,
     $ 0.1208650973866179d-2, -.5395239384953d-5,
     $ 2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*dlog(tmp)-tmp
      ser=1.000000000190015d0
      do 10 j=1,6
       y=y+1.d0
       ser=ser+cof(j)/y
 10   continue
      gammln=tmp+dlog(stp*ser/x)
      RETURN
      END


      DOUBLE PRECISION FUNCTION gammp(a,x)
C Returns the incomplete gamma function P(a,x).
C Taken from Numerical Recipes in Fortran, 2nd ed., Cambridge University Press.
      DOUBLE PRECISION a, x, gammcf, gamser,gln
      if(x.lt.a+1.) then
C use series representation
       call gser(gamser,a,x,gln)
       gammp=gamser
      else
C use continued fraction representation
       call gcf(gammcf,a,x,gln)
       gammp = 1.-gammcf  
      endif
      return
      end

      SUBROUTINE gser(gamser,a,x,gln)
C Taken from Numerical Recipes in Fortran, 2nd ed., Cambridge University Press.
      INTEGER ITMAX, n
      DOUBLE PRECISION a, gamser, gln, x, EPS, ap, del, sum, gammln
      PARAMETER (ITMAX=100, EPS=3.e-7)
      gln=gammln(a)
      ap=a
      sum=1./a
      del=sum
      do 10 n=1,ITMAX
       ap=ap+1.
       del=del*x/ap
       sum=sum+del
       if(abs(del).lt.abs(sum)*EPS) goto 20
 10   continue
 20   gamser=sum*exp(-x+a*dlog(x)-gln)
      return
      end

      SUBROUTINE gcf(gammcf, a,x,gln)
C Taken from Numerical Recipes in Fortran, 2nd ed., Cambridge University Press.
      INTEGER ITMAX, i
      DOUBLE PRECISION a, x, gammcf, gln, EPS, FPMIN,
     $ an,b,c,d,del,h,gammln 
      PARAMETER (ITMAX=100,EPS=3.e-7, FPMIN=1.e-30)
      gln=gammln(a)
      b=x+1.-a
      c=1./FPMIN
      d=1./b
      h=d
      do 10 i=1,ITMAX
       an=-i*(i-a)
       b=b+2.
       d=an*d+b
       if(abs(d).lt.FPMIN) d=FPMIN
       c=b+an/c
       if(abs(c).lt.FPMIN) c=FPMIN
       d=1./d
       del=d*c
       h=h*del
       if(abs(del-1.).lt.EPS) goto 20
 10   continue
 20   gammcf=exp(-x+a*dlog(x)-gln)*h
      return
      end
       
      SUBROUTINE rnorm(Z1,Z2,seed)
C Generates two independent N(0,1) random variables Z1 and Z2.
      double precision U1,U2,Z1,Z2,pi,ran2
      integer seed
      pi=3.14159265
      U1=ran2(seed)
      U2=ran2(seed)
      z1=dcos(2.*pi*u2)*dsqrt((-2.)*dlog(U1))
      Z2=dsin(2.*pi*U2)*dsqrt((-2.)*dlog(U1))
      return
      end

      DOUBLE PRECISION FUNCTION ran2(idum)
C Taken from Numerical Recipes in Fortran, 2nd ed., Cambridge University Press.
C
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      DOUBLE PRECISION AM, EPS, RNMX
C
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     $ IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,
     $ IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
C `Long period (> 2 x 10^18) random number generator of L'Ecuyer with
C Bays-Durham shuffle and added safeguards. Returns a uniform random 
C deviate between 0.0 and 1.0 (exclusive of the endpoint values). Call 
C with idum a negative integer to initialize; thereafter, do not alter 
C idum between successive deviates in sequence. '
C
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/ , iv/NTAB*0/, iy/0/
      if(idum.le.0) then
C initialize
       idum=max(idum,1)
C prevent idum=0
       idum2=idum
       do 100 j=NTAB+8,1,-1
C Load shuffle table (after 8 warm-ups)
        k=idum/IQ1      
        idum=IA1*(idum-k*IQ1)-k*IR1
        if(idum.lt.0) then
         idum=idum+IM1
        endif
        if(j.le.NTAB) then
         iv(j)=idum
        endif
 100   continue
       iy=iv(1)
      endif
C Start here when not initializing.
C Compute idum=mod(IA1*idum,IM1) without overflows by 
C Schrage's method
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if(idum.lt.0) then
       idum=idum+IM1
      endif
      k=idum2/IQ1
C Compute idum2=mode(IA2*idum2,IM2) likewise.
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if(idum2.lt.0) then
       idum2=idum2+IM2
      endif
      j=1+iy/NDIV
C j will be in the range 1:NTAB
C Here idum is shuffled, idum and idum2 are combined to 
C generate output.
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1) then
       iy=iy+IMM1
      endif
      ran2=min(AM*iy,RNMX)
C No endpoint values returned.
      return
      END