      subroutine antith(kmax,k,irsimBIG,it,v0,a)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c ANTITH is a routine that creates a "grid" of antithetic points S on the
c unit sphere, starting from a random basis provided by RNDBASE.
c These routines will ordinarily be called for each observation.  To
c avoid "chatter" when parameter values are changed, the seed for the
c random number generator for each observation should be stored so that
c the same outputs can be regenerated in iteration to parameter estimates.
c
ccc   k=m
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension v0(kmax,irsimBIG),a(kmax,2*kmax+4)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      do 10 i = 1,k
         do 10 j = 1,k
            v0ij = v0(i,j)
            a(i,j) = v0ij
            a(i,k+j) = -v0ij
 10   continue
      if (it .le. 1) then
         go to 999
      else
         xit=it
         xk=k
         j = 2 + intfloor((0.5d0*xit-xk)/(2.d0*xk*(xk-1.d0)))
         dj = j
         do 100 i = 1,k-1
            do 200 ii = i+1,k
               jj = 1
 300           continue
               djj = jj
               dj = djj / dj
               if (jj.lt.j) then
                  do 400 n = 1,k
                     v1n = v0(n,i)*dcos(1.5707963d0*dj)
                     v2n = v0(n,ii)*dsin(1.5707963d0*(1.0d0-dj))
                     a(n,k+1) = v1n + v2n
                     a(n,k+2) = v1n - v2n
                     a(n,k+3) = -v1n + v2n
                     a(n,k+4) = -v1n - v2n
 400              continue
                  jj = jj + 1
                  go to 300
               endif
 200        continue
 100     continue
      endif
 999  continue
      return
      end
      subroutine arse(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   P. Acceptance/Rejection Simulator (ARS)
c      The ARS simulator is given in a version using the truncated
c bilateral exponential comparison distribution (ARSE), and in a
c version using a recursive truncated comparison distribution (ARSR).
c PARM[.,1] is assumed to contain seeds for the uniform random number
c generator RNDUS called by the routines; there is one seed for each
c repetition R.  PARM[R+1,1] is assumed to contain a censoring limit
c for rejection.  Should only compute HC.
c
c The ARS method provides a mechanism for drawing from a conditional
c density when transformations from uniform or standard normal variates
c are not available.  We want to sample from f(x/a).  We generate x
c from g(x).  In our case g(x) = ARSE.  Generate tsi from uniform.
c Generate unconditional f(x) and  bound alpha. We accept x if
c                tsi < f(x)/(g(x)*alpha)
c x has conditional density f(x/a) = HC.  We set x = VD.
c
c NB: Upon exit, the first element of the BLANK common block, contains
c --  the proportion of non-censored draws 
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,irmax),parm(irmax+1,2),
     /          p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          e(mmaxloc,irmaxloc),vd(mmaxloc),
     /          bxmu(mmaxloc),axmu(mmaxloc),temp(mmaxloc),d(mmaxloc),
     /          pa(mmaxloc),pb(mmaxloc),xmatmm(mmaxloc,mmaxloc),
     /          tempt(1,mmaxloc)
      common/times/hseclim
c
      mp1=m+1
      r=ir
      mmaxp1 = mmax + 1
      rscale = 1.0d0/r
      q = 0.0d0
      do 100 j = 1,m
         do 101 i  = 1,m
            xmatmm(i,j) = wi(i,j)
 101     continue
 100  continue
      call invv(xmatmm,xmatmm,mmaxloc,mmaxloc,vd,m,npr,nzr,nnr,rcond,
     *          *999)
c
ccc-----note:vd contains the eigen values of wi in descending order-----
c
      call sorts(vd,m)
      alp = 0.0d0
      do 200 i = 1,m
         wii = w(i,i)
         si = dsqrt(wii+1.d-100)
         xmui = xmu(i)
         bxmui = b(i) - xmui
         axmui = a(i) - xmui
         tempi = -3.0d0 * si
         di = 3.0d0 * si
         t1 = 0.0d0
         t2 = 0.0d0
         t3 = 0.0d0
         t4 = 0.0d0
         if (bxmui.lt.tempi)  t1=1.0d0
         if (axmui.gt.di)  t2=1.0d0
         if (bxmui.lt.di)  t3=1.0d0
         if (axmui.gt.tempi)  t4=1.0d0
         t1 = t1*(bxmui-2.0d0*si)
         t1 = t1*(bxmui-2.0d0*si) + t2*axmui
     *        + (1.0d0-t1-t2)*(tempi+t4*(axmui+di))
         t2 = t1*bxmui + t2*(axmui+2.0d0*si)
     *        + (1.0d0-t1-t2)*(di+t3*(bxmui+tempi))
         bb = 0.0d0
         aa = 0.0d0
         if (t2.lt.0.0d0)  aa=1.0d0
         if (t1.gt.0.0d0)  bb=1.0d0
         si = 1.0d0 - aa - bb
         tempi = (t2-t1)*(0.5d0-si/4.0d0)/(wii+1.d-100)
         alp = alp + tempi**2.0d0
         bxmu(i) = bxmui
         axmu(i) = axmui
         temp(i) = tempi
 200  continue
      alp = alp/(2.0d0*vd(1)+1.d-100)
      call bexp(mmaxloc,m,axmu,temp,pa)
      call bexp(mmaxloc,m,bxmu,temp,pb)
      jl = parm(ir+1,1)
      do 400 j = 1,mp1
         do 401 i = 1,mp1
            hc(i,j) = 0.0d0
 401     continue
 400  continue
      timelim = hsec1990()
      do 500 i = 1,ir
         isd = idint(parm(i,1))
         do 600 j = 1,jl
            do 610 k = 1,m
cvahfix               uk = unifrm(isd)
               uk = randu(isd)
               d(k) = uk*pa(k) + (1.0d0-uk)*pb(k)
 610        continue
            call bexpi(mmaxloc,m,1,d,temp,pa,pb,vd,d)
cvahfix            uzscal = unifrm(isd)
            uzscal = randu(isd)
            call matptv(vd,mmaxloc,wi,mmax,tempt,1,m,1,m)
            scalz1 = 0.0d0
            scalz2 = 0.0d0
            do 620 ii = 1,m
               vdii = vd(ii)
               scalz1 = scalz1 + tempt(1,ii)*vdii
               scalz2 = scalz2 + temp(ii)*dabs(vdii)
 620        continue
            scale = -(scalz1/2.0d0) + scalz2 - alp
            if ( dlog(uzscal) .lt. scale ) then
               q = q + rscale
               jl = j
               go to 88
            endif
 600     continue
 88      continue
         call matpv(wi,mmax,vd,mmaxloc,vd,mmaxloc,m,m,1)
         do 700 j = 1,m
            tempt(1,j) = vd(j)
 700     continue
         call matpv(vd,mmaxloc,tempt,1,xmatmm,mmaxloc,m,1,m)
         hc(1,1) = hc(1,1)+1.d0
         do 800 j=1,m
            vdj = vd(j)
            hc(1,j+1) = hc(1,j+1)+vdj
            hc(j+1,1) = hc(j+1,1)+vdj
            do 801k=1,m
               hc(k+1,j+1) = hc(k+1,j+1)+(xmatmm(k,j) - wi(k,j))*0.5d0
 801        continue
 800     continue
         if ( (hsec1990()-timelim) .gt. hseclim ) then
c           give up
            r = i
            rscale = 1.0d0/r
            go to 99
         endif
 500  continue
 99   continue
      do 910 j = 1,m+1
         do 911 i = 1,m+1
            hc(i,j) = hc(i,j)*rscale
            hu(i,j) = -999.0d0
 911     continue
 910  continue
      p(1) = -999.0d0
ccc      print *, 'arse proportion non-censored', q
      e(1,1)=q
      go to 888
 999  continue
      print *, 'error in routine finding eigen values'
 888  continue
      return
      end
      subroutine arsr(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c see ARSE for documentation.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,irmax),parm(irmax+1,2),
     /          p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          uu(mmaxloc+1,irmaxloc),tt(mmaxloc),
     /          wic(mmaxloc,mmaxloc),v(mmaxloc),vt(1,mmaxloc),
     /          htmm(mmaxloc,mmaxloc)
      common/times/hseclim
c
      mp1 = m + 1
      mmaxp1 = mmax + 1
      mmaxlcp1 = mmaxloc + 1
      r=ir
      rscale = 1.0d0/r
      do 10 j = 1,m
         do 11 i = 1,m
            hc(i,j) = 0.0d0
 11      continue
 10   continue
      jl = parm(ir+1,1)
      chw11 = chw(1,1)
      xmu1 = xmu(1)
      q = 0.0d0
      timelim = hsec1990()
      do 100 i = 1,ir
         isd = idint(parm(i,1))
         do 200 j = 1,jl
            ta = pheev( (a(1)-xmu1) / (chw11+1.d-100) )
            tb = pheev( (b(1)-xmu1) / (chw11+1.d-100) )
            call mrunif(uu,mmaxlcp1,mp1,jl,isd)
            do 210 jk = 1,jl
               do 211 ik = 1,mp1
cvahfix                  uu(i,j) = unifrm(isd)
                  uu(i,j) = randu(isd)
 211           continue
 210        continue
            tt(1) = phinv( uu(1,j)*ta + (1.0d0 - uu(1,j))*tb )
            wgt = tb - ta
            alp = wgt
            do 300 k = 2,m
               km1 = k - 1
               ak = a(k)
               bk = b(k)
               xmuk = xmu(k)
               chwkk = chw(k,k)
               uukj = uu(k,j)
               chw1tt = 0.0d0
               do 301 jk = 1,km1
                  chw1tt = chw1tt + chw(k,jk)*tt(jk)
 301           continue
               ta = pheev( (ak-xmuk-chw1tt) / (chwkk+1.d-100) )
               tb = pheev( (bk-xmuk-chw1tt) / (chwkk+1.d-100) )
               tt(k) = phinv(uukj*ta+(1.0d0-uukj)*tb)
               wgt = wgt*(tb-ta)
               cdfn = pheev((bk-ak)/(2.0d0*chwkk+1.d-100))
               alp = alp*(2.0d0*cdfn - 1.0d0)
 300        continue
            if ( uu(mp1,j) .lt. wgt/alp ) then
               q = q + rscale
               jl = j
               go to 310
            endif
 200     continue
 310     continue
         call matpv(wi,mmax,chw,mmax,wic,mmaxloc,m,m,m)
         call matpv(wic,mmaxloc,tt,mmaxloc,v,mmaxloc,m,m,1)
         do 400 ik = 1,m
            vt(1,ik) = v(ik)
 400     continue
         call matpv(v,mmaxloc,vt,1,htmm,mmaxloc,m,1,m)
         hc(1,1) = hc(1,1)+1.d0
         do 500 j=1,m
            vj = v(j)
            hc(1,j+1) = hc(1,j+1)+vj
            hc(j+1,1) = hc(j+1,1)+vj
            do 501 k=1,m
               hc(k+1,j+1) = hc(k+1,j+1)+(htmm(k,j) - wi(k,j))*0.5d0
 501        continue
 500     continue
         if ( (hsec1990()-timelim) .gt. hseclim ) then
c           give up
            r=i
            rscale = 1.0d0/r
            goto 900
         endif
 100  continue
 900  continue
      do 910 j = 1,m+1
         do 911 i = 1,m+1
            hc(i,j) = hc(i,j)*rscale
            hu(i,j) = -999.0d0
 911     continue
 910  continue
      p(1) = -999.0d0
ccc      print *, 'arsr proportion non-censored', q
      uu(1,1)=q
      return
      end
      subroutine aus(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   S. Asymptotically Unbiased Simulator (AUS)
c      The AUS simulator is obtained as the product of an unbiased
c simulator of H and an asymptotically unbiased simulator of 1/P.  The
c proc is defined with the GHK simulator used for H and for the initial
c independent simulations of P(B)/P(A), where A is the event that the
c first component is in the rectangle bounds.  The proc assumes that
c PARM(.,1) contains a vector of R seeds for the normal random number
c generator RNDNS, PARM(1,2) contains the number of initial GHK
c simulators of P used in the simulation of 1/P, PARM(2,2) contains the
c maximum number of trials using the crude frequency simulator that are
c made before the process is censored.  If PARM(2,2) is made large (say
c 1.e50), then this is computationally the same as the Sequential
c Unbiased Simulator (SUS).  PARM(3,2) contains a seed for RNDUS.
c Gives HU and HC.
c
c This does continuous parameter estimation of 1/p.
c
c NB: Upon exit, the first element of the BLANK common block, contains
c --  the maximum number of CFS draws allowed in AUS.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,ir),parm(irmax+1,2),
     /          hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          f(mmaxloc-1),g(mmaxloc-1,mmaxloc-1),vt(mmaxloc-1),
     /          a2(mmaxloc-1),b2(mmaxloc-1)
c
      common/times/hseclim
c
      mm1 = m - 1
      mmaxlcm1 = mmaxloc - 1
      mp1 = m + 1

      chw11 = chw(1,1)
      xmu1  = xmu(1)

      do 101 j  = 2,m
         f(j-1) = w(j,1)/(w(1,1)+1.d-100)
         do 102 i = 2,m
            g(i-1,j-1) = chw(i,j)
 102     continue
 101  continue
      ta = pheev((a(1)-xmu1)/(chw11+1.d-100))
      tb = pheev((b(1)-xmu1)/(chw11+1.d-100))
      jlim = parm(2,2)
ccc      print *, 'in Aus: max. num. of cfs draws = ', jlim
      isd = idint(parm(3,2))
      ii = 2
      call fghk(ii,mmax,m,xmu,w,wi,chw,a,b,ir,isd,pa,pt,hu)
      pinv = 1.0d0
      q = 1.0d0
      ii = 1
      do 20 j = 1,parm(1,2)
         call fghk(ii,mmax,m,xmu,w,wi,chw,a,b,ir,isd,pa,rr,hc)
         q = q*(1.0d0-pt/(pa+1.d-100))
         pinv = pinv + q
 20   continue
      ppinv = 0.0d0
      timelim = hsec1990()
      do 30 i = 1,ir
         rr = 0.0d0
         isdn = idint(parm(i,1))
         do 40 j = 1,jlim
cvahfix            uu = unifrm(isd)
            uu = randu(isd)
            tt = chw11*phinv(uu*ta + (1.0d0-uu)*tb)
            do 49 k=1,mm1
               a2(k)=ranorm(0.d0,1.d0,isdn)
 49         continue
            do 50 k = 1,mm1
               temp=0.d0
               do 51 kk = 1,mm1
                  temp=temp+g(k,kk)*a2(kk)
 51            continue
               vt(k)=f(k)*tt+temp
 50         continue
            do 600 k = 1,mm1
               temp = vt(k)
               a2k = a(k+1) - xmu(k+1)
               b2k = b(k+1) - xmu(k+1)
               a2(k) = temp - a2k
               b2(k) = temp - b2k
 600        continue
            a2zmin = fmin(a2,mm1)
            b2zmax = fmax(b2,mm1)
ccc            if ( a2zmin .gt. 0.0d0 .and. b2zmax .lt. 0.0d0 ) then
            if (     a2zmin .gt. 0.0d0 .and. b2zmax .lt. 0.0d0
     *          .or. (hsec1990()-timelim) .gt. hseclim  ) then
               jlim = j
               go to 700
            else
               rr = rr + 1.0d0
            endif
 40      continue
 700     continue
         ppinv = ppinv + pinv + q*rr
ccc         print *,'ppinv',ppinv
         if ( (ppinv .gt. 1.d+50) .or.
     *        (hsec1990()-timelim) .gt. hseclim ) then
c           give up either to prevent overflow or because of time:
            ir = i
            go to 900
         endif
 30   continue
 900  continue
ccc      ppinv = ppinv/(ir*pa+1.d-100)
      temp=ir
      if(pa.gt.0.d0) then
         temp=temp*pa
         ppinv = ppinv/temp
         do 910 j = 1,mp1
            do 910 i = 1,mp1
               hc(i,j) = hu(i,j)*ppinv
 910     continue
      else
         do 911 j = 1,mp1
            do 911 i = 1,mp1
               hc(i,j) = -999.d0
 911     continue
      endif
      p = hu(1,1)
c     place jlim=parm(2,2) in common (must use double precision):
      f(1)=parm(2,2)
ccc      print *,'jlim parm(2,2) f(1)',jlim,parm(2,2),f(1)
      return
      end
      subroutine bexp(mmax,m,v,xlam,p)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c bilateral exponential cdf.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension v(mmax),xlam(mmax),p(mmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /         e(mmaxloc,irmaxloc)
c
      do 10 i = 1,m
         if (v(i) .gt. 0.0d0) then
            p(i) = 1.0d0 - 0.5d0*dexp(xlam(i)*(-dabs(v(i))))
         else
            p(i) = 0.5d0*dexp(xlam(i)*(-dabs(v(i))))
         endif
 10   continue
      return
      end
      subroutine bexpi(mmax,m,ir,p,xlam,pa,pb,x,d)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c inverse truncated bilateral exponential cdf.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension p(mmax,ir),xlam(mmax),pa(mmax,ir),pb(mmax,ir),
     /          x(mmax,ir),d(mmax,ir)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common // e(mmaxloc,irmaxloc)
c
      do 10 j = 1,ir
         do 20 i = 1,m
            pij = p(i,j)
            xlami = xlam(i)
            if ( pij .gt. 0.5d0 ) then
               phalf = 1.0d0
            else
               phalf = 0.0d0
            endif
            pp = pij + phalf*(1.0d0-2.0d0*pij)
            x(i,j) = (1.0d0-2.0d0*phalf)*dlog(2.0d0*pp+1.d-100)
     *               /(xlami+1.d-100)
            d(i,j) = xlami*pp /( pb(i,j)-pa(i,j)+1.d-100 )
 20      continue
 10   continue
      return
      end
      subroutine cf(irmax,m,ir,alp,bet,gam,y0,y1,y2)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c see PCF for documentation.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension alp(irmax),bet(irmax),gam(irmax),y0(irmax),
     /          y1(irmax),y2(irmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          t(irmaxloc),z0(irmaxloc),
     /          z1(irmaxloc),z2(irmaxloc)
c
      do 100 i = 1,ir
         alpi = alp(i) + 1.d-100
         gami = gam(i)
         beti = bet(i)
         dalp = dsqrt(alpi)
         ti = dexp(-0.5d0*alpi*(gami-beti/alpi)**2.0d0)
     *             /alpi
         z0i = 2.506628274631 * (pheev(beti/dalp)) / dalp
         y0i = -2.506628274631
     *           *(pheev((beti-alpi*gami)/dalp)) / dalp + z0i
         z1i = (dexp(-beti**2/(2.0d0*alpi))) / alpi
         y1i = z1i - ti + y0i*beti/alpi
         z1i = z1i + z0i*beti/alpi
         z2i = z1i*beti/alpi + z0i/alpi
         y2i = y1i*beti/alpi + y0i/alpi - gami*ti
         alp(i) = alpi
         t(i) = ti
         y0(i) = y0i
         y1(i) = y1i
         y2(i) = y2i
         z0(i) = z0i
         z1(i) = z1i
         z2(i) = z2i
 100  continue
      k = m-1
      if (k .gt. 0) then
         j = 2
 200     continue
         if (j .lt. k+2) then
            dj = j
            do 202 i = 1,ir
               alpi = alp(i)
               beti = bet(i)
               gji = gam(i)**j
               y0i = y1(i)
               y1i = y2(i)
               y2i = y1i*beti/alpi + y0i*dj/alpi
     *                   - gji*t(i)
               z0i = z1(i)
               z1i = z2(i)
               z2i = z1i*beti/alpi + z0i*dj/alpi
               y0(i) = y0i
               y1(i) = y1i
               y2(i) = y2i
               z0(i) = z0i
               z1(i) = z1i
               z2(i) = z2i
 202        continue
            j = j+1
            goto 200
         endif
      endif
      return
      end
      subroutine cfs(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   I. Crude Frequency Simulator (CFS)
c      This simulator assumes that U is an array of standard normal
c variates, independent within each column, and either independent or
c antithetic across columns.  It does not use PARM.  It returns both HU
c and HC.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,m),wi(mmax,m),chw(mmax,m),
     /          u(mmax,ir),parm(irmax+1,2),
     /          p(1),hu(mmax+1,m+1),hc(mmax+1,m+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          vd(mmaxloc,irmaxloc),
     /          f(irmaxloc),xmatmir(mmaxloc,irmaxloc),
     /          vdsq(mmaxloc,mmaxloc),vdt(irmaxloc,mmaxloc),
     /          vdd(mmaxloc)
c
      mp1 = m+1
      r = ir
      rscale = 1.0d0/r
      call matpv(chw,mmax,u,mmax,vd,mmaxloc,m,m,ir)
      do 10 j=1,ir
         do 11 i=1,m
            aa = 0.0d0
            bb = 0.0d0
            axmu = a(i)-xmu(i)
            bxmu = b(i)-xmu(i)
            if (vd(i,j).lt.bxmu) bb = 1.0d0
            if (vd(i,j).gt.axmu) aa = 1.0d0
            xmatmir(i,j) = aa*bb
 11      continue
 10   continue
      sum = 0.0d0
      do 20 j = 1,ir
         fj = fmin(xmatmir(1,j),m)
         f(j) = fj
         sum = sum + fj
         do 21 i = 1,m
            xmatmir(i,j) = vd(i,j)*fj
 21      continue
 20   continue
      p(1) = sum*rscale
      call matpv(wi,mmax,xmatmir,mmaxloc,vd,mmaxloc,m,m,ir)
      do 30 j = 1,m
         sum = 0.0d0
         do 31 i = 1,ir
            temp = vd(j,i)
            vdt(i,j) = temp
            sum = sum + temp
 31      continue
         vdd(j) = sum*rscale
 30   continue
      call matpv(vd,mmaxloc,vdt,irmaxloc,vdsq,mmaxloc,m,ir,m)
      hu(1,1) = p(1)
      do 40 j=1,m
         hu(1,j+1) = vdd(j)
         hu(j+1,1) = vdd(j)
         do 41 i=1,m
            hu(i+1,j+1) = (vdsq(i,j)*rscale - p(1)*wi(i,j))*0.5d0
 41      continue
 40   continue
      if (p(1) .gt. 0.0d0) then
         pscale = 1.0d0/p(1)
         do 50 j=1,mp1
            do 51 i=1,mp1
               hc(i,j)=hu(i,j)*pscale
 51         continue
 50      continue
      else
         do 60 j=1,mp1
            do 61 i=1,mp1
               hc(i,j)=-999.d0
 61         continue
 60      continue
      endif
      return
      end
      subroutine dcs(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   O. Deak Chi-Square Simulator (DCS)
c
c      The DCS assumes that U is an array whose columns are points in
c the unit sphere, antithetic across columns.
c The routine assumes that PARM(1,1) = DET(W).
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,irmax),parm(irmax+1,2),
     /          p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          aasp(mmaxloc+1,irmaxloc),
     /          bbsp(mmaxloc,irmaxloc),qu(irmaxloc),ql(irmaxloc),
     /          temp(1,mmaxloc),tempv(mmaxloc),htmm(mmaxloc,mmaxloc),
     /          tempm(mmaxloc,irmaxloc)
cvah     /          tempm(mmaxloc,irmaxloc),det(2),kpvt(mmaxloc)
c
      common/times/hseclim
c
cvah      dimension inert(3)
c
      dm = m
      mp1 = m + 1
      do 100 j = 1,m
         do 101 i = 1,m
            tempm(i,j) = w(i,j)
 101     continue
 100  continue
cvah  det(w) must be passed as parm(1,1)
      detw=parm(1,1)
cvah      call dsifa(tempm,mmaxloc,m,kpvt,info)
cvah      if(info.eq.0) then
cvah         info=010
cvah         call dsidi(tempm,mmaxloc,m,kpvt,det,inert,tempv,info)
cvah      else
cvah         print *,'W matrix singular: in DCS routine.  Aborting...'
cvah         stop
cvah      endif
cvah      detw = det(1)*10.0d0**det(2)
      dk = 1.0d0/dsqrt(detw + 1.d-100)
      gammo2 = dgamma(dm/2.0d0) + 1.d-100
      dk1 = dk * dgamma((dm+1.d0)/2.d0) * dsqrt(2.0d0) / gammo2
      dk2 = dk * 2.0d0 * dgamma((dm+2.d0)/2.d0) / gammo2
      do 200 j = 1,ir
         do 201 i = 1,m
            ai = a(i)
            bi = b(i)
            xmui = xmu(i)
            uij = u(i,j)
            if ( uij .eq. 0.0d0 ) then
               uu = uij + 1.d-100
            else
               uu = uij
            endif
            sn = 0.0d0
            sp = 0.0d0
            if (uij .lt. 0.0d0) sn = 1.0d0
            if (uij .gt. 0.0d0) sp = 1.0d0
            aasp(i,j) = ((ai-xmui)*sp+(bi-xmui)*sn) / (uu+1.d-100)
            bbsp(i,j) = ((bi-xmui)*sp+(ai-xmui)*sn) / (uu+1.d-100)
 201     continue
         aasp(m+1,j) = 0.0d0
         qlj = fmax(aasp(1,j),mp1)
         quj = fmin(bbsp(1,j),m)
         if (qlj.gt.quj) then
            qu(j) = qlj
         else
            qu(j) = quj
         endif
         ql(j) = qlj
 200  continue
      timelim = hsec1990()
      dfm = m
      dfmp1 = mp1
      dfmp2 = m+2
      do 500 j = 1,mp1
         do 501 i = 1,mp1
            hu(i,j) = 0.0d0
 501     continue
 500  continue
      do 600 j = 1,ir
         call matptv(u(1,j),mmax,wi,mmax,temp,1,m,1,m)
         t = 0.0d0
         do 610 i = 1,m
            t = t + temp(1,i)*u(i,j)
 610     continue
         tr = 1.0d0/(dsqrt(t)+1.d-100)
         tt = tr**m
         if ( qu(j) .gt. ql(j) ) then
            xzql = t*ql(j)**2
            xzqu = t*qu(j)**2
            d0zscale = - chicdf(xzql,dfm,ier)   + chicdf(xzqu,dfm,ier)
            d0 = dk*tt*d0zscale
            d1zscale = - chicdf(xzql,dfmp1,ier) + chicdf(xzqu,dfmp1,ier)
            d1 = dk1*tt*tr*d1zscale
            d2zscale = - chicdf(xzql,dfmp2,ier) + chicdf(xzqu,dfmp2,ier)
            d2 = dk2*tt*tr**2*d2zscale
         else
            d0 = 0.0d0
            d1 = 0.0d0
            d2 = 0.0d0
         endif
         call matpv(wi,mmax,u(1,j),mmax,tempv,mmaxloc,m,m,1)
         do 620 i = 1,m
            temp(1,i) = tempv(i)
 620     continue
         call matpv(tempv,mmaxloc,temp,1,tempm,mmaxloc,m,1,m)
         do 630 jj = 1,m
            do 631 k = 1,m
               tempm(k,jj) = (d2*tempm(k,jj)-d0*wi(k,jj)) / 2.0d0
 631        continue
            tempv(jj) = d1*tempv(jj)
 630     continue
         hu(1,1) = hu(1,1)+d0
         do 640 jj=1,m
            tempvjj = tempv(jj)
            hu(1,jj+1) = hu(1,jj+1)+tempvjj
            hu(jj+1,1) = hu(jj+1,1)+tempvjj
            do 641 i=1,m
               hu(i+1,jj+1) = hu(i+1,jj+1)+tempm(i,jj)
 641        continue
 640     continue
         if ( (hsec1990()-timelim) .gt. hseclim ) then
c           give up
            ir = j
            go to 999
         endif
 600  continue
 999  continue
      r=ir
      rscale= 1.0d0 / r
      do 700 j = 1,mp1
         do 701 i = 1,mp1
            hu(i,j) = hu(i,j)*rscale
 701     continue
 700  continue
      p(1) = hu(1,1)
      if (p(1) .gt. 0.0d0) then
         pscale = 1.0d0/p(1)
         do 800 j=1,mp1
            do 801 i=1,mp1
               hc(i,j)=hu(i,j)*pscale
 801        continue
 800      continue
      else
         do 900 j=1,mp1
            do 901 i=1,mp1
               hc(i,j)=-999.d0
 901        continue
 900     continue
      endif
      return
      end
      subroutine exaf(mmax,m,a,b,xmu,w,xl,p0,dm,dl)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c approximate answers by numerical quadrature.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension a(mmax),b(mmax),xmu(mmax),xl(mmax),w(mmax,mmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /         s(mmaxloc),ss(mmaxloc,mmaxloc),
     /         xmatmm(mmaxloc,mmaxloc),rmat(mmaxloc,mmaxloc),
     /         acf(mmaxloc),bcf(mmaxloc),
     /         xmup(mmaxloc),xlp(mmaxloc),wp(mmaxloc,mmaxloc)
ccc     /         xlzt(1,mmaxloc),
ccc     /         xlpzt(1,mmaxloc),xlmat(mmaxloc,mmaxloc),
ccc     /         xlpmat(mmaxloc,mmaxloc)

c
      mm1 = m - 1
      do 100 i = 1,mm1
         xmup(i) = xmu(i)
         xli = xl(i)
         xlp(i) = xli
ccc         xlzt(1,i) = xli
ccc         xlpzt(1,i) = xli
 100  continue
      xmup(m) = xmu(m) + 1.d-3
      xlm = xl(m)
ccc      xlzt(1,m) = xlm
      xlm = xlm + 1.d-3
      xlp(m) = xlm
ccc      xlpzt(1,m) = xlm
ccc      call matpv(xl,mmax,xlzt,1,xlmat,mmaxloc,m,1,m)
ccc      call matpv(xlp,mmaxloc,xlpzt,1,xlpmat,mmaxloc,m,1,m)
      do 201 j = 1,m
         do 202 i = 1,m
ccc            wp(i,j) = w(i,j) - xlmat(i,j) + xlpmat(i,j)
            wp(i,j) = w(i,j) - xl(i)*xl(j) + xlp(i)*xlp(j)
 202     continue
 201  continue
      call pf(mmaxloc,m,a,b,xmu ,w ,p0)
      call pf(mmaxloc,m,a,b,xmup,w ,p1m)
      call pf(mmaxloc,m,a,b,xmu ,wp,p1w)
      call pf(mmaxloc,m,a,b,xmup,w ,dm)
      call pf(mmaxloc,m,a,b,xmu ,wp,dl)
      if(p0.ne.-999.d0.and.p1m.ne.-999.d0) then
         dm = (dm-p0)/1.d-3
      else
         dm=-999.d0
      endif
      if(p0.ne.-999.d0.and.p1w.ne.-999.d0) then
         dl = (dl-p0)/1.d-3
      else
         dl=-999.d0
      endif
c     extra test inserted by vah, sep 22, 92, to allow for inaccuracies
c     in the fortran pheev3 code (used inside pf):
      if(p0.eq.0.d0.and.p1m.eq.0.d0) then
         dm=-999.d0
      endif
      if(p0.eq.0.d0.and.p1w.eq.0.d0) then
         dl=-999.d0
      endif
      return
      end
      subroutine ffact(mmax,m,x,a,b,xmu,d,xl,y0,y1,y2)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c integrand in one-factor model.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      double precision nrmpdf
      dimension x(40),a(mmax),b(mmax),xmu(mmax),d(mmax),xl(mmax),
     /          y0(40),y1(40),y2(40)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //    xt(1,40),xxx(mmaxloc,40)
c
      mm1 = m - 1
      do 10 i = 1,40
         xt(1,i) = x(i)
 10   continue
      call matpv(xl,mmax,xt,1,xxx,mmax,m,1,40)
      do 101 j = 1,40
         tt = 1.0d0
         do 102 i = 1,mm1
            xmui = xmu(i)
            xxxij = xxx(i,j)
            di = d(i)
            xhs = (b(i)-xmui-xxxij) / di
            rhs = (a(i)-xmui-xxxij) / di
            tt = tt*(pheev(xhs) - pheev(rhs))
 102     continue
         xmum = xmu(m)
         xxxmj = xxx(m,j)
         dm = d(m)
         xhs = (b(m)-xmum-xxxmj) / dm
         rhs = (a(m)-xmum-xxxmj) / dm
         y0(j) = tt*(pheev(xhs) - pheev(rhs))
         ttt = nrmpdf(xhs,0.0d0,1.0d0) - nrmpdf(rhs,0.0d0,1.0d0)
         y1(j) = -tt*ttt / dm
         y2(j) = y1(j)*x(j)
 101  continue
      return
      end
      subroutine fdl(mmax,m,h,xl,d)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c return specific elements from the H matrix.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension h(mmax+1,mmax+1),xl(mmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      mp1 = m + 1
      d = 0.0d0
      do 10 i = 1,m
         d = d + h(mp1,i+1)*xl(i)
 10   continue
      d = d * 2.0d0
      return
      end
      subroutine fghk(ii,mmax,m,xmu,w,wi,chw,a,b,ir,iseed,pm,p,hu)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c see GHK for documentation.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          p(1),hu(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          f(mmaxloc-1),g(mmaxloc-1,mmaxloc-1),vt(mmaxloc-1),
     /          a2(mmaxloc-1),b2(mmaxloc-1),
ccc--local variables for fghk, previous ones are from aus which calls it.
     /          u(mmaxloc,irmaxloc),
     /          tt(mmaxloc,irmaxloc),ttt(irmaxloc,mmaxloc),
     /          wgt(1,irmaxloc),chw1(1,mmaxloc),
     /          v(mmaxloc,irmaxloc),chw1tt(1,irmaxloc),
     /          xmatmm(mmaxloc,mmaxloc),vv(mmaxloc)
c
      mmaxp1 = mmax + 1
      mp1 = m + 1
      r = ir
      rscale = 1.0d0/r
      do 10 j = 1,ir
         do 11 i = 1,m
cvahfix            u(i,j) = unifrm(iseed)
            u(i,j) = randu(iseed)
 11      continue
 10   continue
      xmu1 = xmu(1)
      chw11 = chw(1,1)
      tascale = pheev((a(1)-xmu1)/(chw11+1.d-100))
      tbscale = pheev((b(1)-xmu1)/(chw11+1.d-100))
      do 20 j = 1,ir
         u1j = u(1,j)
ccc should never happen:
ccc         if(u1j*tascale+(1.0d0-u1j)*tbscale.le.0.d0) then
ccc            print *,'error! phinv in FGHK (AUS/SUS)'
ccc         endif
         tt(1,j) = phinv(u1j*tascale+(1.0d0-u1j)*tbscale)
         pa = tbscale - tascale
         wgt(1,j) = pa
 20   continue
      do 30 i = 2,m
         im1 = i - 1
         do 31 k = 1,im1
            chw1(1,k) = chw(i,k)
 31      continue
         call matpv(chw1,1,tt,mmaxloc,chw1tt,1,1,im1,ir)
         xmui = xmu(i)
         chwii = chw(i,i)
         do 32 j = 1,ir
            chw1tt1j = chw1tt(1,j)
            uij = u(i,j)
            ta = pheev((a(i)-xmui-chw1tt1j)/(chwii+1.d-100))
            tb = pheev((b(i)-xmui-chw1tt1j)/(chwii+1.d-100))
ccc should never happen:
ccc            if(uij*ta+(1.0d0-uij)*tb.le.0.d0) then
ccc               print *,'error! phinv in FGHK (AUS/SUS)'
ccc            endif
            tt(i,j) = phinv(uij*ta+(1.0d0-uij)*tb)
            wgt(1,j) = wgt(1,j)*(tb-ta)
 32      continue
 30   continue
      sum = 0.0d0
      do 40 i = 1,ir
         sum = sum + wgt(1,i)
 40   continue
      p(1) = sum/r
      if (ii .eq. 1) then
         do 45 j = 1,mp1
            do 46 i = 1,mp1
               hu(i,j) = 0.0d0
 46         continue
 45      continue
      else
         call matpv(chw,mmax,tt,mmaxloc,v,mmaxloc,m,m,ir)
         do 50 j=1,m
            sum = 0.0d0
            do 51 i=1,ir
              val = v(j,i)*wgt(1,i)
              sum = sum + val
              ttt(i,j) = val
 51         continue
            vv(j) = sum
 50      continue
         call matpv(wi,mmax,v,mmaxloc,tt,mmaxloc,m,m,ir)
         call matpv(ttt,irmaxloc,wi,mmax,ttt,irmaxloc,ir,m,m)
         call matpv(tt,mmaxloc,ttt,irmaxloc,xmatmm,mmaxloc,m,ir,m)
         call matpv(wi,mmax,vv,mmaxloc,vv,mmaxloc,m,m,1)
      hu(1,1) = p(1)
      do 60 j=1,m
         vvj = vv(j)
         hu(1,j+1) = vvj*rscale
         hu(j+1,1) = vvj*rscale
         do 61 i=1,m
            hu(i+1,j+1) = (xmatmm(i,j)*rscale - p(1)*wi(i,j))*0.5d0
 61      continue
 60   continue
      endif
      pm = pa
      return
      end
      subroutine fi(mmax,irmax,m,ir,a,b,dkap,fi0,fi1,fi2)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c FI is a proc that returns partial moments of the standard normal.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension a(mmax,irmax),b(mmax,irmax),dkap(mmax,irmax),
     /          fi0(mmax,irmax),fi1(mmax,irmax),fi2(mmax,irmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      do 101 j = 1,ir
         do 102 i = 1,m
            aij =  a(i,j)
            bij =  b(i,j)
            dkapij =  dkap(i,j)
            fi0ij = pheev(bij) - pheev(aij)
            fi1ij = -dkapij*fi0ij - fi0ij
            fi2ij = fi0ij - dkapij*fi1ij - (bij - dkapij)*pheev(bij)
     *              + (aij - dkapij)*pheev(aij)
            fi0(i,j) =  fi0ij
            fi1(i,j) =  fi1ij
            fi2(i,j) =  fi2ij
 102     continue
 101  continue
      return
      end
      subroutine ghk(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   M. Geweke-Hajivassiliou-Keane Simulator (GHK)
c      The GHK procedure assumes that U is an array of uniform [0,1]
c variates, independent within columns, and in general independent
c between columns. It does not use PARM.  It returns both HU and HC.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,m),wi(mmax,m),chw(mmax,m),
     /          u(mmax,ir),parm(irmax+1,2),
     /          p(1),hu(mmax+1,m+1),hc(mmax+1,m+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          tt(mmaxloc,irmaxloc),wgt(irmaxloc),
     /          chw1(1,mmaxloc),chw1tt(1,irmaxloc),
     /          ttt(irmaxloc,mmaxloc),v(mmaxloc,irmaxloc),
     /          xmatmm(mmaxloc,mmaxloc),vdd(mmaxloc)
c
      mp1 = m+1
      mmaxp1 = mmax + 1
      r = ir
      rscale = 1.0d0/r
      xmu1 = xmu(1)
      chw11 = chw(1,1)
      tascale = pheev((a(1)-xmu1)/(chw11+1.d-100))
      tbscale = pheev((b(1)-xmu1)/(chw11+1.d-100))
      do 10 j = 1,ir
         u1j = u(1,j)
ccc should never happen:
ccc         if(u1j*tascale+(1.0d0-u1j)*tbscale.le.0.d0) then
ccc            print *,'GHK: warning! argument to phinv LE 0'
ccc            print *,u1j,tascale,tbscale
ccc            print *,u1j*tascale+(1.0d0-u1j)*tbscale
ccc         endif
         tt(1,j) = phinv(u1j*tascale+(1.0d0-u1j)*tbscale)
         wgt(j) = tbscale - tascale
 10   continue
      do 20 i = 2, m
         im1 = i-1
         do 21 k = 1,im1
            chw1(1,k) = chw(i,k)
 21      continue
         call matpv(chw1,1,tt,mmaxloc,chw1tt,1,1,im1,ir)
         xmui = xmu(i)
         chwii = chw(i,i)
         do 22 j = 1,ir
            chw1tt1j = chw1tt(1,j)
            uij = u(i,j)
            ta = pheev((a(i)-xmui-chw1tt1j)/(chwii+1.d-100))
            tb = pheev((b(i)-xmui-chw1tt1j)/(chwii+1.d-100))
            tt(i,j) = phinv(uij*ta+(1.0d0-uij)*tb)
            wgt(j) = wgt(j) * (tb-ta)
 22      continue
 20   continue
      call matpv(chw,mmax,tt,mmaxloc,v,mmaxloc,m,m,ir)
      sum = 0.0d0
      do 25 j = 1,ir
         sum = sum + wgt(j)
 25   continue
      p(1) = sum*rscale
      do 30 j = 1,m
      sumt = 0.0d0
         do 31 i = 1,ir
            val = v(j,i)*wgt(i)
            sumt = sumt + val
            ttt(i,j) = val
 31      continue
         vdd(j) = sumt
 30   continue
      call matpv(wi,mmax,v,mmaxloc,tt,mmaxloc,m,m,ir)
      call matpv(ttt,irmaxloc,wi,mmax,ttt,irmaxloc,ir,m,m)
      call matpv(tt,mmaxloc,ttt,irmaxloc,xmatmm,mmaxloc,m,ir,m)
      call matpv(wi,mmax,vdd,mmaxloc,vdd,mmaxloc,m,m,1)
      hu(1,1) = p(1)
      do 40 j=1,m
         vddj = vdd(j)
         hu(1,j+1) = vddj*rscale
         hu(j+1,1) = vddj*rscale
         do 41 i=1,m
            hu(i+1,j+1) = (xmatmm(i,j)*rscale - p(1)*wi(i,j))*0.5d0
 41      continue
 40   continue
      if (p(1) .gt. 0.0d0) then
         pscale = 1.0d0/p(1)
         do 50 j=1,mp1
            do 51 i=1,mp1
               hc(i,j)=hu(i,j)*pscale
 51         continue
 50      continue
      else
         do 60 j=1,mp1
            do 61 i=1,mp1
               hc(i,j)=-999.d0
 61         continue
 60      continue
      endif
      return
      end
      subroutine gss(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   Q. Gibbs's Sampler Simulator (GSS)
c      The GSS simulates HC, but does not simulate P or HU.  It assumes
c that PARM[1,1] gives the number of rounds to be used in the 
c construction of the random draws (JL) and that U is an MxJL array of 
c uniform [0,1] random numbers.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,irmax),parm(1),
     /          p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter (mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /         v(mmaxloc),wimkk(mmaxloc-1),wmkk(mmaxloc-1),
     /         wmkmki(mmaxloc-1,mmaxloc-1),vbar(mmaxloc),
     /         vvpbar(mmaxloc,mmaxloc),condsdk(mmaxloc),
     /         condmukf(mmaxloc,mmaxloc-1),ui(mmaxloc,irmaxloc),
     /         mkind(mmaxloc-1)
c
      common/times/hseclim
c
      mp1 = m + 1
      mm1 = m - 1
      r = ir
      jl = parm(1)

      do 10 j = 1,jl
         do 11 i=1,m
            ui(i,j)=u(i,j)
 11      continue
 10   continue

      do 100 ii = 1,m
         v(ii) = (b(ii)+a(ii))/2.0d0 - xmu(ii)
 100  continue

      timestr = hsec1990()

      do 200 k=1,m
         if(k.eq.1) then
            do 201 n=1,mm1
               mkind(n)=n+1
 201        continue
         elseif(k.eq.m) then
            do 202 n=1,mm1
               mkind(n)=n
 202        continue
         else
            do 203 n=1,k-1
               mkind(n)=n
 203        continue
            do 204 n=k,mm1
               mkind(n)=n+1
 204        continue
         endif
         do 300 ni=1,mm1
            wmkk(ni)=w(mkind(ni),k)
            wimkk(ni)=wi(mkind(ni),k)
 300     continue
         do 301 nj=1,mm1
            do 302 ni=1,mm1
               wmkmki(ni,nj)=wi(mkind(ni),mkind(nj))
     *                     - wimkk(ni)*wimkk(nj)/(wi(k,k)+1.d-100)
 302        continue
 301     continue
         temp=0.d0
         do 303 nj=1,mm1
            do 304 ni=1,mm1
               temp=temp+wmkk(ni)*wmkmki(ni,nj)*wmkk(nj)
 304        continue
 303     continue
         condsdk(k)=dsqrt(w(k,k)-temp)
         do 305 nj=1,mm1
            temp=0.d0
            do 306 ni=1,mm1
               temp=temp+wmkk(ni)*wmkmki(ni,nj)
 306        continue
            condmukf(k,nj)=temp
 305     continue
 200  continue

      do 4001 nj=1,m
         vbar(nj)     =0.d0
         do 4002 ni=1,m
            vvpbar(ni,nj)=0.d0
 4002    continue
 4001 continue

      do 400 irepl=1,ir
         if (irepl.gt.1) then
            iseedgss=ui(1,1)*1.e+9
            do 401 nj=1,jl
               do 402 ni=1,m
                  ui(ni,nj)=randu(iseedgss)
 402           continue
 401        continue
         endif
         do 500 j=1,jl
            do 600 k=1,m
               if(k.eq.1) then
                  do 2001 n=1,mm1
                     mkind(n)=n+1
 2001             continue
               elseif(k.eq.m) then
                  do 2002 n=1,mm1
                     mkind(n)=n
 2002             continue
               else
                  do 2003 n=1,k-1
                     mkind(n)=n
 2003             continue
                  do 2004 n=k,mm1
                     mkind(n)=n+1
 2004             continue
               endif
               csdk=condsdk(k)
               cmuk=0.d0
               do 403 nj=1,mm1
                  mk=mkind(nj)
ccc         print *,'nj,mk,k,mmaxloc',nj,mk,k,mmaxloc
cerror???                  cmuk=cmuk+condmukf(k,mk)*v(mk)
                  cmuk=cmuk+condmukf(k,nj)*v(mk)
 403           continue
               ta=pheev((a(k)-xmu(k)-cmuk)/(csdk+1.d-100))
               tb=pheev((b(k)-xmu(k)-cmuk)/(csdk+1.d-100))
               uikj=ui(k,j)
               v(k)=cmuk+csdk*phinv(uikj*tb+(1.d0-uikj)*ta+1.d-100)
 600        continue
 500     continue
         do 501 nj=1,m
            vbar(nj)=vbar(nj)+v(nj)
            do 502 ni=1,m
               vvpbar(ni,nj)=vvpbar(ni,nj)+v(ni)*v(nj)
 502        continue
 501     continue
         if((hsec1990()-timestr) .gt. hseclim) then
c           give up because of time-limit
            ir=i
            r=ir
            goto 4000
         endif
 400  continue

 4000 continue

c     vah: potential to exploit SYMMETRY

      hc(1,1)=1.d0
      do 700 nj=1,m
         temp=0.d0
         do 701 i=1,m
            temp=temp+vbar(i)*wi(i,nj)
 701     continue
         temp=temp/r
         hc(1,nj+1)=temp
         hc(nj+1,1)=temp
         do 702 ni=1,m
            res=0.d0
            do 703 k=1,m
               do 704 j=1,m
                  res=res+wi(ni,j)*(vvpbar(j,k)/r-w(j,k))*wi(k,nj)/2.d0
 704           continue
 703        continue
            hc(ni+1,nj+1)=res
 702     continue
 700  continue

      do 800 nj=1,mp1
         do 801 ni=1,mp1
            hu(ni,nj)=-999.d0
  801    continue
  800 continue

      p(1)=-999.d0

      return
      end
      subroutine nden(mmax,irmax,m,ir,vd,wi,chw,q)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c NDEN is the multivariate normal density.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension vd(mmax,irmax),wi(mmax,mmax),chw(mmax,mmax),q(irmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //  e(mmaxloc,irmaxloc)
c
      dkk = 0.39894228040143 ** m
      prod = 1.0d0
      do 10 i = 1,m
         prod = prod * chw(i,i)
 10   continue
      dkk = dkk / prod
      call matpv(wi,mmax,vd,mmax,e,mmaxloc,m,m,ir)
      do 100 j = 1,ir
         sum = 0.0d0
         do 101 i = 1,m
            sum = sum + vd(i,j)*e(i,j)
 101     continue
         q(j)= dkk * dexp(-0.5d0*sum)
 100  continue
      return
      end
      subroutine ndenist(mmax,irmax,m,ir,vd,wi,chw,q)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension vd(mmax,irmax),wi(mmax,mmax),chw(mmax,mmax),q(irmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /     e(mmaxloc,irmaxloc),s(mmaxloc),vdc(mmaxloc,irmaxloc),
     /     ymat(mmaxloc,irmaxloc)
c
      dkk = 0.39894228040143 ** m
      prod = 1.0d0
      do 10 i = 1,m
         prod = prod * chw(i,i)
 10   continue
      dkk = dkk / (prod+1.d-100)
      call matpv(wi,mmax,vd,mmax,e,mmaxloc,m,m,ir)
      do 100 j = 1,ir
         sum = 0.0d0
         do 101 i = 1,m
            sum = sum + vd(i,j)*e(i,j)
 101     continue
         q(j)= dkk * dexp(-0.5d0*sum)
 100  continue
      return
      end
      subroutine nise(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   J. Normal Importance Sampling Simulator (NIS)
c      This simulator provides HC = HU/P.  It provides two alternative
c simulators, corresponding to the exponential (NISE) and independent
c truncated normal (NIST) comparison distributions.  These simulators
c assume U is an array of uniform random (0,1) variates.  These procs
c call the routine NDEN that returns a vector of values of the
c multivariate normal density.
c      Procedure NISE does not use C or PARM.  It calls the bilateral
c exponential CDF BEXP and its inverse BEXPI from LIB_A.G.  It
c guarantees P positive. Provides both HU and HC.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,ir),parm(irmax+1,2),
     /          p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /        e(mmaxloc,irmaxloc),s(mmaxloc),
     /        diagw(mmaxloc),bxmu(mmaxloc),axmu(mmaxloc),
     /        t1(mmaxloc),t2(mmaxloc),t3(mmaxloc),t4(mmaxloc),
     /        temp(mmaxloc),pa(mmaxloc,irmaxloc),pb(mmaxloc,irmaxloc),
     /        wgt(irmaxloc),xmatirm(irmaxloc,mmaxloc),q(irmaxloc),
     /        humm(mmaxloc,mmaxloc),vdd(mmaxloc),vd(mmaxloc,irmaxloc),
     /        d(mmaxloc,irmaxloc)
c
      mp1 = m + 1
      do 10 i = 1,m
         diagw(i) = w(i,i)
         s(i) = dsqrt(diagw(i))
         bxmu(i) = b(i) - xmu(i)
         axmu(i) = a(i) - xmu(i)
         temp(i) = -3.0d0 * s(i)
         q(i) = 3.0d0 * s(i)
         vdd(i) = 0.0d0
 10   continue
      do 20 i = 1,m
         t1(i) = 0.0d0
         t2(i) = 0.0d0
         t3(i) = 0.0d0
         t4(i) = 0.0d0
         if (bxmu(i).lt.temp(i))  t1(i)=1.0d0
         if (axmu(i).gt.q(i))  t2(i)=1.0d0
         if (bxmu(i).lt.q(i))  t3(i)=1.0d0
         if (axmu(i).gt.temp(i))  t4(i)=1.0d0
         t1(i) = t1(i)*(bxmu(i)-2.0d0*s(i))
     *             + t2(i)*axmu(i)
     *             + (1.0d0-t1(i)-t2(i))*
     *                        (temp(i)+t4(i)*(axmu(i)+q(i)))
         t2(i) = t1(i)*bxmu(i)
     *             + t2(i)*(axmu(i)+2.0d0*s(i))
     *             + (1.0d0-t1(i)-t2(i))*
     *                        (q(i)+t3(i)*(bxmu(i)+temp(i)))
 20   continue
      do 30 i = 1,m
         t3(i) = 0.0d0
         t4(i) = 0.0d0
         if (t2(i).lt.vdd(i))  t4(i)=1.0d0
         if (t1(i).gt.vdd(i))  t3(i)=1.0d0
         s(i) = 1.0d0 - t4(i) - t3(i)
         temp(i) = 1.d-100
     *            + (t2(i)-t1(i))*(0.5d0-s(i)/4.0d0)/(diagw(i)+1.d-100)
 30   continue
      call bexp(mmaxloc,m,axmu,temp,axmu)
      call bexp(mmaxloc,m,bxmu,temp,bxmu)
      do 40 j = 1,ir
         do 41 i =1,m
            pa(i,j) = axmu(i)
            pb(i,j) = bxmu(i)
            d(i,j) = u(i,j)*pa(i,j) + (1.0d0-u(i,j))*pb(i,j)
 41      continue
 40   continue
      r = ir
      rscale = 1.0d0/r
      call bexpi(mmaxloc,m,ir,d,temp,pa,pb,vd,d)
      call nden(mmaxloc,irmaxloc,m,ir,vd,wi,chw,q)
      sum = 0.0d0
      do 50 j = 1,ir
         prod = 1.0d0
         do 51 i = 1,m
            prod = prod*d(i,j)
 51      continue
         wgt(j) = q(j) / (prod + 1.0d-100)
         sum = sum + wgt(j)
         do 52 i = 1,m
            pa(i,j) = wgt(j)*vd(i,j)
 52      continue
 50   continue
      p(1) = sum*rscale
      call matpv(wi,mmax,pa,mmaxloc,pb,mmaxloc,m,m,ir)
      call matptv(vd,mmaxloc,wi,mmax,xmatirm,irmaxloc,m,ir,m)
      call matpv(pb,mmaxloc,xmatirm,irmaxloc,humm,mmaxloc,m,ir,m)
      do 60 j = 1, m
         do 61 i = 1,m
            humm(i,j) = humm(i,j)*rscale - p(1)*w(i,j)
 61      continue
 60   continue
      do 70 j = 1,m
         sum = 0.0d0
         do 71 i = 1,ir
            sum = sum + pa(j,i)
 71      continue
         vdd(j) = sum*rscale
 70   continue
      call matpv(wi,mmax,vdd,mmaxloc,vdd,mmaxloc,m,m,1)
      hu(1,1) = p(1)
      do 80 j=1,m
         hu(1,j+1) = vdd(j)
         hu(j+1,1) = vdd(j)
         do 81 i=1,m
            hu(i+1,j+1) = humm(i,j)*0.5d0
 81      continue
 80   continue
      if (p(1) .gt. 0.0d0) then
         pscale = 1.0d0/p(1)
         do 90 j=1,mp1
            do 91 i=1,mp1
               hc(i,j)=hu(i,j)*pscale
 91         continue
 90      continue
      else
         do 100 j=1,mp1
            do 101 i=1,mp1
               hc(i,j)=-999.d0
 101        continue
 100     continue
      endif
      return
      end
      subroutine nist(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c see NISE for documentation.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,irmax),parm(irmax+1,2),
     /          p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /         e(mmaxloc,irmaxloc),s(mmaxloc),vdc(mmaxloc,irmaxloc),
     /         ymat(mmaxloc,irmaxloc),xmatirm(irmaxloc,mmaxloc),
     /         y(irmaxloc),wgt(irmaxloc),vd(mmaxloc,irmaxloc),
     /         temp(mmaxloc),humm(mmaxloc,mmaxloc),vdd(mmaxloc)
c
      mp1 = m + 1
      r = ir
      rscale = 1.0d0/r
      do 10 i = 1,m
         temp(i) = a(i) - xmu(i)
         vdd(i) = b(i) - xmu(i)
 10   continue
      call rndt(mmax,irmax,m,ir,u,w,temp,vdd,vd,y)
      call ndenist(mmax,irmax,m,ir,vd,wi,chw,wgt)
      sum = 0.0d0
      do 20 j = 1,ir
         wgt(j) = wgt(j)/(y(j)+1.0d-100)
         sum = sum + wgt(j)
         do 21 i = 1,m
            ymat(i,j) = wgt(j)*vd(i,j)
 21      continue
 20   continue
      p(1) = sum*rscale
      do 30 j = 1,m
         sum = 0.0d0
         do 31 i = 1,ir
            sum = sum + ymat(j,i)
 31      continue
         vdd(j) = sum*rscale
 30   continue
      call matpv(wi,mmax,ymat,mmaxloc,ymat,mmaxloc,m,m,ir)
      call matptv(vd,mmaxloc,wi,mmax,xmatirm,irmaxloc,m,ir,m)
      call matpv(ymat,mmaxloc,xmatirm,irmaxloc,humm,mmaxloc,m,ir,m)
      do 40 j=1,m
         do 41 i=1,m
            humm(i,j) = humm(i,j)*rscale - wi(i,j)*p(1)
 41      continue
 40   continue
      call matpv(wi,mmax,vdd,mmaxloc,vdd,mmaxloc,m,m,1)
      hu(1,1) = p(1)
      do 50 j=1,m
         hu(1,j+1) = vdd(j)
         hu(j+1,1) = vdd(j)
         do 51 i=1,m
            hu(i+1,j+1) = humm(i,j)*0.5d0
 51      continue
 50   continue
      if (p(1) .gt. 0.0d0) then
         pscale = 1.0d0/p(1)
         do 60 j=1,mp1
            do 61 i=1,mp1
               hc(i,j)=hu(i,j)*pscale
 61         continue
 60      continue
      else
         do 70 j=1,mp1
            do 71 i=1,mp1
               hc(i,j)=-999.d0
 71         continue
 70      continue
      endif
      return
      end
      subroutine pcf(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   N. Parabolic Cylinder Function Simulator (PCF)
c      The PCF proc assumes that U is an array whose column vectors are
c points in the intersection of the unit sphere and the negative
c orthant.
c The routine assumes that PARM(1,1) = DET(W).
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,irmax),parm(irmax+1,2),
     /          p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          aa(mmaxloc,irmaxloc),bb(mmaxloc),
     /          vwv(irmaxloc),bwv(irmaxloc),
     /          svec(irmaxloc),ro(irmaxloc),
     /          y0(irmaxloc),y1(irmaxloc),y2(irmaxloc),tempv(mmaxloc)
c
      mp1 = m + 1
      r = ir
c
cvah  det(w) must be passed as parm(1,1)
      detw=parm(1,1)
c
c     dkp
      dm = m
      dkp = dgamma(dm/2.0d0) * (2.0d0**(dm/2.0d0-1.0d0)) * dsqrt(detw)
      dkp = 1.0d0/(dkp + 1.d-100)
c
c     bb=wi*(mu-b)   (mx1)
      do 20 i=1,m
         sum=0.d0
         do 21 k=1,m
            sum = sum + wi(i,k)*(xmu(k) - b(k))
 21      continue
         bb(i) = sum
 20   continue
c
c     aa=wi*u   (mxr)
      call matpv(wi,mmax,u,mmax,aa,mmaxloc,m,m,ir)
c
c     vwv=sumc(u.*aa)   (rx1)
      do 30 j = 1,ir
         sum = 0.0d0
         do 31 k = 1,m
            sum = sum + u(k,j)*aa(k,j)
 31      continue
         vwv(j) = sum
 30   continue
c
c     bwv=u'bb   (rx1)
      do 40 j=1,ir
         sum=0.d0
         do 41 k=1,m
            sum = sum + u(k,j)*bb(k)
 41      continue
         bwv(j) = sum
 40   continue
c
c     svec, ro
c     sval=-(xmu-b)'bb/2
      sval=0.d0
      do 50 k=1,m
         sval = sval - (xmu(k)-b(k))*bb(k)
 50   continue
      sval = sval/2.d0
      do 51 j=1,ir
c        svec
         temp = sval+bwv(j)**2/(2.d0*vwv(j)+1.d-100)
         svec(j) = dkp*dexp(temp)
c        ro
         valmin = (a(1)-b(1))/(u(1,j)+1.d-100)
         do 52 i=2,m
            temp = (a(i)-b(i))/(u(i,j)+1.d-100)
            if (temp.lt.valmin) valmin = temp
 52      continue
         ro(j) = valmin
 51   continue
c
c     y0, y1, y2
c     note slight difference from gauss: k=m-1 is calculated INSIDE cf:
      call cf(irmaxloc,m,ir,vwv,bwv,ro,y0,y1,y2)
c     p
      sum = 0.0d0
      do 55 i = 1,ir
         sum = sum + y0(i)*svec(i)
 55   continue
      p(1) = sum/r
c
c     tempv
      do 60 i = 1,m
         sum = 0.d0
         do 61 j = 1,m
            do 62 k = 1,ir
               sum = sum + wi(i,j)*u(j,k)*y1(k)*svec(k)
 62         continue
 61      continue
         tempv(i) = sum/r
 60   continue
c
c     hu
      twotom=2.d0**m
      twotomp1=twotom*2.d0
      hu(1,1) = p(1)/twotom
      do 70 j=1,m
         temp = (tempv(j) - p(1)*bb(j)) / twotom
         hu(1,j+1) = temp
         hu(j+1,1) = temp
c        should exploit SYMMETRY!!!
         do 71 i=1,m
            sum = 0.d0
            do 72 k=1,ir
               sum = sum + aa(i,k)*aa(k,i)*y2(k)*svec(k)
 72         continue
            temp = 
     * (- tempv(i)*bb(j) - tempv(j)*bb(i) 
     * + sum/r + p(1)*(bb(i)*bb(j) - wi(i,j)) ) / twotomp1
            hu(i+1,j+1) = temp
 71      continue
 70   continue
c
c     hc
      p(1) = hu(1,1)
      if ( p(1) .gt. 0.0d0 ) then
         do 80 j = 1,mp1
            do 81 i =1,mp1
               hc(i,j) = hu(i,j)/p(1)
 81        continue
 80     continue
      else
         do 90 j = 1,mp1
            do 91 i = 1,mp1
               hc(i,j) = -999.0d0
 91        continue
 90     continue
      endif
      return
      end
      subroutine pf(mmax,m,aa,bb,xmu,w,p)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c used by EXAF.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension aa(mmax),bb(mmax),xmu(mmax),w(mmax,mmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /         s(mmaxloc),ss(mmaxloc,mmaxloc),
     /         xmatmm(mmaxloc,mmaxloc),rmat(mmaxloc,mmaxloc),
     /         a(mmaxloc),b(mmaxloc)
c
      do 100 i = 1,m
         s(i) = dsqrt(w(i,i))
         ss(i,i) = 1.0d0/s(i)
 100  continue
      call matpv(ss,mmaxloc,w,mmax,xmatmm,mmaxloc,m,m,m)
      call matpv(xmatmm,mmaxloc,ss,mmaxloc,rmat,mmaxloc,m,m,m)
      do 110 i = 1,m
         a(i) = (aa(i)-xmu(i)) / s(i)
         b(i) = (bb(i)-xmu(i)) / s(i)
 110  continue
      if (m.eq.2) then
         r = rmat(1,2)
         a1 = a(1)
         a2 = a(2)
         b1 = b(1)
         b2 = b(2)
         p = pheev2(b1,b2,r) - pheev2(a1,b2,r)
     *       - pheev2(b1,a2,r) + pheev2(a1,a2,r)
      else if (m.eq.3) then
         r1 = rmat(1,2)
         r2 = rmat(2,3)
         r3 = rmat(3,1)
         a1 = a(1)
         a2 = a(2)
         a3 = a(3)
         b1 = b(1)
         b2 = b(2)
         b3 = b(3)
         p = pheev3(b1,b2,b3,r1,r2,r3) - pheev3(a1,b2,b3,r1,r2,r3)
         p = p - pheev3(b1,a2,b3,r1,r2,r3) - pheev3(b1,b2,a3,r1,r2,r3)
         p = p + pheev3(a1,a2,b3,r1,r2,r3) + pheev3(b1,a2,a3,r1,r2,r3)
         p = p + pheev3(a1,b2,a3,r1,r2,r3) - pheev3(a1,a2,a3,r1,r2,r3)
      else
         p = -999.0d0
      endif
      return
      end
      subroutine pker(z,mmax,irmax,m,ir,prodzwgt)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c see PKF for documentation.
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c     Polynomial kernel function; PKFS calls this routine.
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension z(mmax,irmax),prodzwgt(irmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
c      common //
c     /          onesmir(mmaxloc,irmaxloc),
c     /          zerosmir(mmaxloc,irmaxloc),xm1mir(mmaxloc,irmaxloc)
c
c      do 10 j = 1,ir
c         do 11 i = 1,m
c            onesmir(i,j) = 1.0d0
c            zerosmir(i,j) = 0.0d0
c            xm1mir(i,j) = - 1.0d0
c 11      continue
c 10   continue
c      call hadcomp(z,mmax,onesmir,mmaxloc,onesmir,mmaxloc,m,ir,'lt')
c      call hadcomp(z,mmax,zerosmir,mmaxloc,zerosmir,mmaxloc,m,ir,'lt')
c      call hadcomp(z,mmax,xm1mir,mmaxloc,xm1mir,mmaxloc,m,ir,'lt')
      do 101 j = 1,ir
         prod = 1.0d0
         do 102 i = 1,m
            z1 = 0.0d0
            z0 = 0.0d0
            zm1 = 0.0d0
            if (z(i,j).lt.(1.0d0)) z1=1.0d0
            if (z(i,j).lt.(0.0d0)) z0=1.0d0
            if (z(i,j).lt.(-1.0d0)) zm1=1.0d0
            factor = 0.5d0*(1.0d0-z(i,j)*(2.0d0-z(i,j)))*z1
     *               - z(i,j)**2*z0
     *               + 0.5d0*( 1.0d0 + z(i,j)*(2.0d0+z(i,j)) )*zm1
     *
            prod = prod*factor
 102     continue
        prodzwgt(j) = prod
 101  continue
      return
      end
      subroutine pkfs(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   K. Polynomial Kernel-Smoothed Frequency Simulator (PKF)
c      The KFS simulator uses the kernel (31), defined by the function
c PKER. It assumes that U is an array of standard normal variates,
c independent within rows.  It does not use W.  PARM[1,1] contains a
c window width parameter.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,irmax),parm(irmax+1,2),
     /          p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          vd(mmaxloc,irmaxloc),
     /          vdbb(mmaxloc,irmaxloc),vdaa(mmaxloc,irmaxloc),
     /          pkerbb(irmaxloc),pkeraa(irmaxloc),
     /          wivd(mmaxloc,irmaxloc),aat(irmaxloc,mmaxloc),
     /          humm(mmaxloc,mmaxloc),
     /          vdd(mmaxloc)
c
      mmaxlcp1 = mmaxloc + 1
      mp1 = m + 1
      mmaxp1 = mmax + 1
      r = ir
      rscale = 1.d0/r
c
      call matpv(chw,mmax,u,mmax,vd,mmaxloc,m,m,ir)
      do 10 j = 1,ir
         do 11 i = 1,m
            vdbb(i,j) = (vd(i,j)-(b(i)-xmu(i)))/parm(1,1)
            vdaa(i,j) = (vd(i,j)-(a(i)-xmu(i)))/parm(1,1)
 11      continue
 10   continue
      call pker(vdbb,mmaxloc,irmaxloc,m,ir,pkerbb)
      call pker(vdaa,mmaxloc,irmaxloc,m,ir,pkeraa)
      call matpv(wi,mmax,vd,mmax,wivd,mmaxloc,m,m,ir)
      sum = 0.0d0
      do 20 j = 1,ir
         fval = pkerbb(j) - pkeraa(j)
         sum = sum + fval
         do 21 i = 1,m
            aat(j,i) = wivd(i,j)*fval
 21      continue
 20   continue
      p(1) = sum*rscale
      call matpv(wivd,mmaxloc,aat,irmaxloc,humm,mmaxloc,m,ir,m)
      do 30 j=1,m
         sum = 0.0d0
         do 31 i=1,ir
            sum = sum + aat(i,j)
 31      continue
         vdd(j) = sum*rscale
 30   continue
      hu(1,1) = p(1)
      do 40 j=1,m
         hu(1,j+1) = vdd(j)
         hu(j+1,1) = vdd(j)
         do 41 i=1,m
            hu(i+1,j+1) = (humm(i,j)*rscale - wi(i,j)*p(1))*0.5d0
 41      continue
 40   continue
      if (p(1) .gt. 0.0d0) then
         pscale = 1.0d0/p(1)
         do 50 j = 1,mp1
            do 51 i = 1,mp1
               hc(i,j) = hu(i,j)*pscale
 51         continue
 50      continue
      else
         do 60 j = 1,mp1
            do 61 i = 1,mp1
               hc(i,j) = -999.0d0
 61         continue
 60      continue
      endif
      return
      end
      subroutine reswrite(cfil,line,irec)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c     direct access version
c     Reswrite writes the character string line to the end of the cfil.
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      character*(*) cfil,line
      open(10,file=cfil,access='direct',form='formatted',recl=130)
ccc      print *,'inside new: line',line
      write(10,'(a)',rec=irec) line(1:ltrim(line))
      close(10)
      end
      subroutine rndbase(kmax,k,iseed1,v,t)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c RNDBASE is a routine that draws a random basis in K space, where K is
c the dimension of the random coefficient vector.
c RNDBASE takes as input a dimension K and returns a K by K array V that
c is column orthonormal and random, bordered by a column of random
c lengths of K-dim. normal random vectors.  SEED1 is a seed for the
c random number generator.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension v(kmax,k),t(kmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /           vt(mmaxloc,mmaxloc),vvt(mmaxloc,mmaxloc)
c
      do 100 j=1,k
         do 101 i=1,k
            v(i,j) = ranorm(0.0d0,1.0d0,iseed)
            vt(j,i) = v(i,j)
 101     continue
 100  continue
      call matpv(v,kmax,vt,mmaxloc,vvt,mmaxloc,k,k,k)
 888  continue
      s = 0.0d0
      do 110 i=1,k
         s = s + v(i,1)**2
 110  continue
      if (s .lt. 1.d-16) then
         do 120 i=1,k
            v(i,1) = ranorm(0.0d0,1.0d0,iseed1)
 120     continue
         go to 888
      endif
      scale = 1.0d0 / (dsqrt(s)+1.d-100)
      do 200 i = 1,k
         v(i,1) = scale * v(i,1)
         t(i) = vvt(i,i)
 200  continue
      do 300 kk = 2,k
 999     continue
         do 310 i = 1,kk
            rhs = 0.0d0
            do 311 j = 1,k
               rhs = rhs + v(j,kk)*v(j,i)
               v(j,kk) = v(j,kk) - v(j,i)*rhs
 311        continue
 310     continue
         s = 0.0d0
         do 320 i = 1,k
            s = s + v(i,kk)**2
 320     continue
         if (s.lt. 1.d-16) then
            call mrnorm(v(1,kk),kmax,k,1,iseed1)
            go to 999
         endif
         scale = 1.0d0 / (dsqrt(s)+1.d-100)
         do 330 i = 1,k
            v(i,kk) = scale * v(i,kk)
 330     continue
 300  continue
      return
      end
      subroutine rndt(mmax,irmax,m,ir,u,w,a,b,v,y)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c RNDT is a proc for drawing an array from truncated independent
c normals that have the same means and variances as the multivariate
c density.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension a(mmax),b(mmax),w(mmax,mmax),u(mmax,irmax),
     /          v(mmax,irmax),y(irmax)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /         e(mmaxloc,irmaxloc),s(mmaxloc),vdc(mmaxloc,irmaxloc),
     /         ymat(mmaxloc,irmaxloc)
c
      do 10 j = 1,m
         do 10 i = 1,m
            if (i.eq.j) then
               s(i) = dsqrt(w(i,j)) + 1.d-100
            endif
 10   continue
      dkk = 0.39894228040143
      do 40 j = 1,ir
         prod = 1.0d0
         do 50 i = 1,m
            bs= pheev( b(i)/s(i) )
            as= pheev( a(i)/s(i) )
            vdc(i,j) = phinv( u(i,j)*bs+(1.0d0-u(i,j))*as )
            v(i,j) = s(i)*vdc(i,j)
ccc---truncated normal density
            sdk = dkk/s(i)
            ymat(i,j) = sdk*dexp(-0.5d0*vdc(i,j)*vdc(i,j))
     *                 /(bs-as+1.d-100)
            prod = prod * ymat(i,j)
 50      continue
         y(j) = prod
 40   continue
      return
      end
      subroutine sds(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   L. Stern Decomposition Simulator (SDS)
c      The SDS proc assumes that PARM[1,1] contains a scalar l > 0 such
c that W-lI is positive definite, and assumes that C is a Cholesky
c factor of W-lI. It assumes that U is an array of standard normal
c variates, independent within each column.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,irmax),parm(irmax+1,2),
     /          p(1),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      include 'hdim.inc'
      parameter(mmaxloc=mmaxim,irmaxloc=irmaxim)
c
      common //
     /          v(mmaxloc,irmaxloc),ta(mmaxloc,irmaxloc),
     /          tb(mmaxloc,irmaxloc),tk(mmaxloc,irmaxloc),
     /          fb0(mmaxloc,irmaxloc),fb1(mmaxloc,irmaxloc),
     /          fb2(mmaxloc,irmaxloc),qi1(mmaxloc,irmaxloc),
     /          qi2(mmaxloc,irmaxloc),s(mmaxloc,irmaxloc),
     /          dd(mmaxloc,mmaxloc),tt(mmaxloc),
     /          tempmm(mmaxloc,mmaxloc),ttt(1,mmaxloc)
c
      mp1 = m + 1
      mmaxp1 = mmax + 1
      mmaxlcp1 = mmaxloc + 1
      r = ir
      rscale = 1.0d0 / r
      xlambda = dsqrt(parm(1,1))
      call matpv(chw,mmax,u,mmax,v,mmaxloc,m,m,ir)
      do 101 j = 1,ir
         do 102 i = 1,m
            tb(i,j) = (b(i)-xmu(i)-v(i,j)) / xlambda
            ta(i,j) = (a(i)-xmu(i)-v(i,j)) / xlambda
            tk(i,j) = -2.0d0*xmu(i) / xlambda
 102     continue
 101  continue
      call fi(mmaxloc,irmaxloc,m,ir,ta,tb,tk,fb0,fb1,fb2)
      do 110 j = 1,ir
         do 120 i = 1,m
            if ( fb0(i,j) .eq. 0.0d0 ) then
               fb0(i,j) = .000000001
            endif
 120     continue
 110  continue
      do 201 j = 1,ir
         do 202 i = 1,m
            qi1(i,j) = xlambda*fb1(i,j)
            qi2(i,j) = xlambda**2*fb2(i,j)
            s(i,j) = qi1(i,j)/fb0(i,j)
 202     continue
 201  continue
      do 210 j = 1,mp1
         do 211 i = 1,mp1
            hu(i,j) = 0.0d0
 211     continue
 210  continue
      do 300 j = 1,ir
         qq = 1.0d0
         do 310 jk=1,m
            qq = qq*fb0(jk,j)
            do 311 ik=1,m
               if (ik.eq.jk) then
                  dd(ik,jk) = qi2(ik,j)/fb0(ik,j)
               else
                  dd(ik,jk) = 0.0d0
               endif
 311        continue
 310     continue
         call matpv(wi,mmax,s(1,j),mmaxloc,tt,mmaxloc,m,m,1)
         call matpv(wi,mmax,dd,mmax,tempmm,mmax,m,m,m)
         call matpv(tempmm,mmax,wi,mmax,tempmm,mmax,m,m,m)
         do 320 ik = 1,m
            ttt(1,ik) = tt(ik)
 320     continue
         call matpv(tt,mmaxloc,ttt,1,dd,mmaxloc,m,1,m)
         hu(1,1) = hu(1,1)+qq
         scale = qq*0.5d0
         do 330 jk=1,m
            hu(1,jk+1) = hu(1,jk+1)+tt(jk)*qq
            hu(jk+1,1) = hu(jk+1,1)+tt(jk)*qq
            do 331 i=1,m
               hu(i+1,jk+1) = hu(i+1,jk+1)
     *                      +(tempmm(i,jk)-wi(i,jk)+dd(i,jk))*scale
 331        continue
 330     continue
 300  continue
      do 400 j = 1,mp1
         do 401 i = 1,mp1
            hu(i,j) = hu(i,j)*rscale
 401     continue
 400  continue
      p(1) = hu(1,1)
      if (p(1) .gt. 0.0d0) then
         pscale = 1.0d0/p(1)
         do 500 j=1,mp1
            do 501 i=1,mp1
               hc(i,j)=hu(i,j)*pscale
 501        continue
 500     continue
      else
         do 600 j=1,mp1
            do 601 i=1,mp1
               hc(i,j)=-999.d0
 601        continue
 600     continue
      endif
      return
      end
      subroutine sim(mmax,m,xl,e1,p,hu,hc,ec,pp,dm,xlm,dl,xll)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c procedure to accumulate results.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension xl(mmax),hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
      pp = p
      dm = hu(m+1,1)
      xlm = hc(m+1,1)
      if ( hu(m+1,1) .ne. -999.0d0 ) then
         call fdl(mmax,m,hu,xl,dl)
      else
         dl = -999.0d0
      endif
      if ( hc(m+1,1) .ne. -999.0d0 ) then
         call fdl(mmax,m,hc,xl,xll)
      else
         xll = -999.0d0
      endif
      e2 = hsec1990()
      ec = e2 - e1
      e1 = e2
      return
      end
      subroutine sus(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c   R. Sequentially Unbiased Simulator (SUS)
c      The SUS simulator is obtained as the product of an unbiased
c simulator of H and an unbiased simulator of 1/P.
c Computationally, if PARM(2,2) is made large (say 1.e50), SUS is the
c same as the Asymptotically Unbiased Simulator (AUS).
c
c NB: Upon exit, the first element of the BLANK common block, contains
c --  the maximum number of CFS draws allowed in AUS
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension xmu(mmax),a(mmax),b(mmax),
     /          w(mmax,mmax),wi(mmax,mmax),chw(mmax,mmax),
     /          u(mmax,ir),parm(irmax+1,2),
     /          hu(mmax+1,mmax+1),hc(mmax+1,mmax+1)
c
c     simply calls aus with a very large PARM(2,2)

ccc      parm(2,2) = 1.0d50
      parm(2,2) = 1.0d7
      call aus(mmax,irmax,m,xmu,w,wi,chw,a,b,ir,u,parm,p,hu,hc)
      return
      end
