      program harness
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   A. Overview
c      This section contains a series of FORTRAN procedures  for
c simulation  of multivariate normal rectangle probabilities, and  the
c derivatives  of  these probabilities with respect to the mean and
c covariances of the normal  distribution.  Table 2 indexes the
c procedures in this section.
c
c                           TABLE 2.  FORTRAN PROCEDURES
c
c      SECTION   MNEUMONIC*             DESCRIPTION
c
c         H        CFS          Crude Frequency Simulator
c         I        NIS          Normal Importance Sampling Simulator
c                               Procs NISE and NIST use exponential and
c                               truncated normal comparison
c                               distributions, respectively
c         J        KFS          Polynomial Kernel Smoothed Frequency
c                               Simulator.  The proc is PKFS
c         K        SDS          Stern Decomposition Simulator
c         L        GHK          Geweke-Hajivassiliou-Keane Simulator
c         M        PCF          Parabolic Cylinder Function Simulator
c         N        DCS          Deak chi-square Simulator
c         O        ARS          Acceptance/Rejection Simulator
c                               Procs ARSE and ARSR use exponential and
c                               recursive truncated normal comparison
c                               distributions, respectively
c         P        GSS          Gibbs Sampler Simulator
c         Q        AUS          Approximate (OR Sequential) Unbiased
c                               Simulator
c         R        SUS          Sequential Unbiased Simulator
c         S        LIB_A        Statistical Functions
c         T        LIB_B        Spherical Transformations and
c                               Antithetics
c         U        LIB_C        Miscellaneous Utilities
c
c * The procs are located in libraries named *.G, where '*' denotes the
c mneumonic.
c
c      Each of the simulator procedures requires the following standard
c inputs; the interpretation of some inputs may vary from routine to
c routine, and not all are used in some routines:
c
c      M         Dimension of the multivariate normal
c      MU        Mean of multivariate normal, a M x 1 vector
c      W         Covariance matrix of multivariate normal, a M x M
c                array
c      WI        Inverse of W
c      C         Lower triangular Cholesky factor of W, a M x M array
c      A         Lower bound of rectangle, a M x 1 vector [When the
c                lower bound is -infinity, set A = (-1.0E10)*ONES(M,1)]
c      B         Upper bound of rectangle, a M x 1 vector [When the
c                upper bound is +infinity, set B = (+1.0E10)*ONES(M,1)]
c      R         Number of repetitions
c      U         Random variates, a M x R array
c      PARM      Parameters and constants for the simulation
c
c The simulators all return {P,HU,HC}, where P is the scalar rectangle
c probability, HU is the (M+1) x (M+1) array of unconditional partial
c moments (6), and HC is the (M+1) x (M+1) array of conditional moments
c (7).  Parts of the output not provided by a simulator are set to
c -999.
c
c Interpretations of input parameters
c -----------------------------------
c ARSR,ARSE
c ---------
c These procs assume that PARM[.,1] contains seeds for the uniform
c random number generator RNDUS called by the routines; there is one
c seed for each repetition R.  PARM[R+1,1] is assumed to contain a
c censoring limit for rejection.
c ARSR and ARSE use  U = UUBIG;                              (uniforms)
c
c AUS
c ---
c The proc assumes that PARM(.,1) contains a vector of R seeds for the
c normal random number generator RNDNS, PARM(1,2) contains the number
c of initial GHK simulators of P used in the simulation of 1/P,
c PARM(2,2) contains the maximum number of trials using the crude
c frequency simulator that are made before the process is censored.  If
c PARM(2,2) is made large (say 1.e50), then this is computationally the
c same as the Sequential Unbiased Simulator (SUS).  PARM(3,2) contains
c a seed for RNDUS.
c AUS uses  U = UUBIG[.,1:R];                                (uniforms)
c
c CFS
c ---
c This proc does not use PARM.
c CFS uses  U = UNBIG[.,1:R]                         (standard normals)
c
c DCS
c ---
c The routine assumes that PARM(1,1) = DET(W).
c DCS uses  U = USBIG[.,1:R];       (random normal grid on unit sphere)
c
c GHK
c ---
c This proc does not use PARM.
c GHK uses  U = UUBIG[.,1:R];                                (uniforms)
c
c GSS
c ---
c It assumes that PARM[1,1] gives the number of rounds to be used in
c the construction of the random draws.
c GSS uses  U = UUBIG[.,1:R];                                (uniforms)
c
c NISE,NIST
c ---------
c These procs do not use PARM.
c NISE and NIST use  U = UUBIG[.,1:R];                       (uniforms)
c
c PCF
c ---
c The routine assumes that PARM(1,1) = DET(W).
c PCF uses  U = USSBIG[.,1:R]; (random normal grid on unit sphere and
c                              orthant)
c
c PKF
c ---
c This proc does not use W.
c It assumes that PARM[1,1] contains a window width parameter.
c KFS uses  U = UNBIG[.,1:R];                        (standard normals)
c
c SDS
c ---
c The SDS proc assumes that PARM[1,1] contains a scalar l > 0 such that
c W-lI is positive definite.
c SDS uses  U = UNBIG[.,1:R];                        (standard normals)
c Non-standard feature:
c --------------------
c SDS assumes that C is a Cholesky factor of W-lI, not of W.
c
c SUS
c ---
c This procedure is the same the Asymptotically Unbiased Simulator
c (AUS) with PARM(2,2) made large (1.e50).
c The proc assumes that PARM(.,1) contains a vector of R seeds for the
c normal random number generator RNDNS, PARM(1,2) contains the number
c of initial GHK simulators of P used in the simulation of 1/P,
c PARM(2,2) contains the maximum number of trials using the crude
c frequency simulator that are made before the process is censored.
c This is made very large.
c PARM(3,2) contains a seed for RNDUS.
c SUS uses  U = UUBIG[.,1:R];                               (uniforms)
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c      Test Harness
c      ------------
c      The following program provides a simple harness for testing the
c procs above.  A one-factor model is used to generate probabilities,
c with gaussian quadrature used for evaluation.  The accuracy of this
c quadrature can be checked against symmetric cases where the
c probabilities are known.
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      double precision nrmpdf
      character*2   resfil
      character*130 results
      character*12  cfilCFS,cfilNISE,cfilNIST,cfilKFS,cfilSDS,
     /              cfilGHK,cfilPCF,cfilDCS,cfilARSE,cfilARSR,
     /              cfilGSS,cfilAUS,cfilSUS
c
      include 'hdim.inc'
c
      parameter(
     * larse = mmaxim*irmaxim+mmaxim*mmaxim+8*mmaxim,
     * larsr = mmaxim*irmaxim+irmaxim+2*mmaxim*mmaxim+3*mmaxim,
     * laus  = 4*mmaxim*irmaxim+2*irmaxim+2*mmaxim*mmaxim-2*mmaxim+3,
     * lcfs  = 3*mmaxim*irmaxim+irmaxim+mmaxim*mmaxim+mmaxim,
     * ldcs  = 3*mmaxim*irmaxim+3*irmaxim+mmaxim*mmaxim+3*mmaxim+2,
     * lghk  = 3*mmaxim*irmaxim+2*irmaxim+mmaxim*mmaxim+2*mmaxim,
     * lgss  = mmaxim*irmaxim+3*mmaxim*mmaxim+3*mmaxim-2,
     * lnise = 6*mmaxim*irmaxim+9*mmaxim+2*irmaxim+mmaxim*mmaxim,
     * lnist = 5*mmaxim*irmaxim+2*irmaxim+mmaxim*mmaxim+3*mmaxim,
     * lpcf  = 4*mmaxim*irmaxim+11*irmaxim+mmaxim*mmaxim+4*mmaxim+2,
     * lpkf  = 5*mmaxim*irmaxim+2*irmaxim+mmaxim*mmaxim+mmaxim,
     * lsds  = 10*mmaxim*irmaxim+2*mmaxim*mmaxim+2*mmaxim,
     * lsus  = 4*mmaxim*irmaxim+2*irmaxim+2*mmaxim*mmaxim-2*mmaxim+3)
      parameter(
     * lcom=lsds)
cLaheyBug     * lcom  = max(larse,larsr,laus ,lcfs ,ldcs ,lghk ,lgss ,
cLaheyBug     *             lnise,lnist,lpcf ,lpkf ,lsds ,lsus))
      common // dummy(lcom)

      common /times/ hseclim
c
      dimension hu(mmaxim+1,mmaxim+1),hc(mmaxim+1,mmaxim+1),
     /          rsim(13),parm(irmaxim+1,2),
     /          xmu(mmaxim),a(mmaxim),b(mmaxim),xl(mmaxim),d(mmaxim),
     /          s(mmaxim),p(1),pvec(mmaxim),
     /          w(mmaxim,mmaxim),wi(mmaxim,mmaxim),chw(mmaxim,mmaxim),
     /          tempmxm(mmaxim,mmaxim),tempmx1(mmaxim),
     /          temp1xm(1,mmaxim),
     /          t(mmaxim),tmm1(mmaxim-1),
     /          uuBIG(mmaxim,irmaxim),
     /          us(mmaxim,mmaxim),uss(mmaxim,2*mmaxim+4),
     /          unBIG(mmaxim,irmaxim),unBIG2(mmaxim,irmaxim),
     /          unzsum(irmaxim),unzsumt(1,irmaxim),
     /          usBIG(mmaxim,irmaxim),ussBIG(mmaxim,irmaxim),
     /          aw(mmaxim,mmaxim),
     /          p40(40),t40(40),tt40(40),
     /          pp0(16),dm0(16),xlm0(16),dl0(16),xll0(16),rep0(16),
     /          el0(16),aa(16,7),res(7),
     /          qw(20,2),qw2(40,2),el(16),rep(16),pp(16),dm(16),
     /          xlm(16),dl(16),xll(16),vec40(40),vec40zl(40)
      dimension index(mmaxim)
      data qw/
     *     0.048613160848618d0,0.146006911993027d0,0.243907555937767d0,
     *     0.342668563127518d0,0.442667961120606d0,0.544320702552795d0,
     *     0.648093938827515d0,0.754526972770691d0,0.864258110523224d0,
     *     0.978062629699707d0,1.096907615661621d0,1.222035169601440d0,
     *     1.355093598365784d0,1.498357176780701d0,1.655115008354187d0,
     *     1.830419778823853d0,2.032696485519409d0,2.277773618698120d0,
     *     2.601793527603149d0,3.127621173858643d0,
     *     0.038752973079681d0,0.038519907742739d0,0.038055181503296d0,
     *     0.037361584603786d0,0.036443289369345d0,0.035305824130774d0,
     *     0.033956021070480d0,0.032402005046606d0,0.030653120949864d0,
     *     0.028719885274768d0,0.026613922789693d0,0.024347903206944d0,
     *     0.021935453638434d0,0.019391084089875d0,0.016730098053813d0,
     *     0.013968503102660d0,0.011122924275696d0,0.008210529573262d0,
     *     0.005249142181128d0,0.002260638633743d0/
c
      write(6,'(a/)')
     *' *-----------------------------------------------------------*',
     *' |                  Simulation Algorithms                    |',
     *' |    for Orthant Normal Probabilities and Derivatives:      |',
     *' |               A Harness Test Program Run                  |',
ccc     *' |              version 1.0, March 22, 1992                  |',
     *' |              version 1.1, July 11, 1992                   |',
     *' |           Written by Vassilis Hajivassiliou               |',
     *' |              with the aid of Yoosoon Chang                |',
     *' |            This Program is Part of the Paper:             |',
     *' |         Simulation of Normal Orthant Probabilities:       |',
     *' |                   Methods and Programs                    |',
     *' |          by Hajivassiliou, McFadden, and Ruud             |',
ccc     *' |                     February 1992                         |',
     *' |                       July 1992                           |',
     *' *-----------------------------------------------------------*'
      call machine(results)
      write(6,'(/a)')
     *' --------------------------'
      write(6,'(a)')
     *' Computer Used for this Run: '//results(1:ltrim(results))
      write(6,'(a/)')
     *' --------------------------'
      call when(6)
c
c     The following allows for the first line in the input file
c     required by the  gauss program, giving the number of numeric
c     parameters in the input file. (See the GAUSS program.)
      print *, 'Hit <Return> when ready:'
      read(*,'(a)') results
      print *, ' '
c
 5    print *, 'enter input dimension m : '
      read(*,*) temp
      m=temp
      print *,'***',m
      if (m.lt.2) then
         print *, 'm should be >= 2 '
         goto 5
      elseif (m.gt.mmaxim) then
         print *, 'm cannot exceed ',mmaxim
         print *, 'which was the limit at compilation time'
         goto 5
      endif
      print *, 'enter seed for random number generator: '
      read(*,*) temp
      iseed1=temp
      print *,'***',iseed1
      iseed = iseed1
c
      print *,'Enter total number of Monte Carlo repetitions: '
      read(*,*) temp
      mcrepet=temp
      print *,'***',mcrepet
c
      write(6,1)
 1    format(' enter 1 to calculate ENDOGENOUSLY number of',
     /       ' simulations for each method,' /
     /       ' enter 0 to specify number for each method yourself: ')
      read(*,*) temp
      iendog=temp
      print *,'***',iendog
c
      if ( iendog .eq. 1 ) then
 10      print *, 'Num. of simul. for CFS ',
     /            '(will calculate endogenously for other methods): '
         read(*,*) temp
         irsimCFS=temp
         print *,'***',irsimCFS
         if (irsimCFS .lt. m+2 ) then
            print *,'Num. of simul. for CFS should be >= m+2 '
            goto 10
         endif
         print *, 'Num. of simul. for endogenous timing of methods ',
     /            '(irsimTRY): '
         read(*,*) temp
         irsimTRY=temp
         print *,'***',irsimTRY
         if ( irsimTRY .gt. irsimCFS ) then
            print *, 'ERROR!! irsimTRY cannot exceed irsimCFS'
            goto 10
         endif
         irsimBIG = irsimCFS
         rsimBIG = irsimBIG
      else
 15      print *, 'Number of repetitions for CFS:                  '
         read(*,*) temp
         irsimCFS=temp
         print *,'***',irsimCFS
         rsim(1) = irsimCFS
         print *, 'Number of repetitions for NIS_exponential:      '
         read(*,*) temp
         irsimNSE=temp
         print *,'***',irsimNSE
         rsim(2) = irsimNSE
         print *, 'Number of repetitions for NIS_truncated_normal: '
         read(*,*) temp
         irsimNST=temp
         print *,'***',irsimNST
         rsim(3) = irsimNST
         print *, 'Number of repetitions for KFS:                  '
         read(*,*) temp
         irsimKFS=temp
         print *,'***',irsimKFS
         rsim(4) = irsimKFS
         print *, 'Number of repetitions for SDS:                  '
         read(*,*) temp
         irsimSDS=temp
         print *,'***',irsimSDS
         rsim(5) = irsimSDS
         print *, 'Number of repetitions for GHK:                  '
         read(*,*) temp
         irsimGHK=temp
         print *,'***',irsimGHK
         rsim(6) = irsimGHK
         print *, 'Number of repetitions for PCF:                  '
         read(*,*) temp
         irsimPCF=temp
         print *,'***',irsimPCF
         rsim(7) = irsimPCF
         print *, 'Number of repetitions for DCS:                  '
         read(*,*) temp
         irsimDCS=temp
         print *,'***',irsimDCS
         rsim(8) = irsimDCS
         print *, 'Number of repetitions for ARS_exponential:      '
         read(*,*) temp
         irsimARE=temp
         print *,'***',irsimARE
         rsim(9) = irsimARE
         print *, 'Number of repetitions for ',
     /            'ARS_recursive_truncated_normal:                 '
         read(*,*) temp
         irsimARR=temp
         print *,'***',irsimARR
         rsim(10) = irsimARR
         print *, 'Number of repetitions for GSS:                  '
         read(*,*) temp
         irsimGSS=temp
         print *,'***',irsimGSS
         rsim(11) = irsimGSS
         print *, 'Number of repetitions for AUS:                  '
         read(*,*) temp
         irsimAUS=temp
         print *,'***',irsimAUS
         rsim(12) = irsimAUS
         print *, 'Number of repetitions for SUS:                  '
         read(*,*) temp
         irsimSUS=temp
         print *,'***',irsimSUS
         rsim(13) = irsimSUS
         print *, '  '
         irsimBIG = fmax(rsim,13)
         if (irsimBIG .lt. m+2 ) then
            print *, 'For at least one method, num. of repetition ',
     *               'should be >= m+2 '
            goto 15
         elseif (irsimBIG .gt. irmaxim) then
            print *,
     *      '*********************************************',
     *      '* For at least one method, you have requested',
     *      '* TOO LARGE a number of simulations (',irsimBIG,')',
     *      '* which exceeds the LIMIT at compilation (',irmaxim,')',
     *      '*********************************************'
            goto 15
         endif
         rsimBIG = irsimBIG
      endif
      print *, 'Time limit (in seconds) for single run of any method: '
      read(*,*) seconds
      print '(1x,a,f10.3)','***',seconds
      hseclim = seconds * 100.d0
c
      write(6,2)
 2    format(' Do you want a special case with p analytic? ' /
     /       3x,'enter 0 = No' /
     /       3x,'      1 = Factor loading on m only' /
     /       3x,'      2 = Orthant prob. with full symmetry'/
     /       3x,'      3 = Arbitrary omega')
      read(*,*) temp
      icase=temp
      print *,'***',icase
c
      if ( icase .eq. 2 ) then
         do 40 i = 1,m
            xmu(i) = 0.0d0
            b(i) = 0.0d0
            a(i) = -25.0d0
            xl(i) = 1.0d0
            d(i) = 1.0d0
 40      continue
      else
         write(6,3)
 3       format(' Do you want automatic generation of data? ' /
     /          ' Enter 1=Yes, 0=No')
         read(*,*) temp
         iautogen=temp
         print *,'***',iautogen
         if ( iautogen .eq. 1 ) then
            call mrnorm(xmu,mmaxim,m,1,iseed)
            call fillm(b,mmaxim,1,0.0d0)
            do 20 i = 1,m
               a(i) = -2.0d0 - dabs(xmu(i))
 20         continue
            if ( icase .eq. 0 ) then
               call mrnorm(xl,mmaxim,m,1,iseed)
               call mdscal(xl,mmaxim,xl,mmaxim,m,1,-3.0d0)
            else if (icase .eq. 1) then
               call fillm(xl,mmaxim,1,0.0d0)
               xl(m) = -3.0d0 * ranorm(1.0d0,0.0d0,iseed)
            else if (icase .eq. 3) then
               call fillm(xl,mmaxim,1,-999.0d0)
            endif
            if (icase .eq. 3) then
               call fillm(d,mmaxim,1,-999.0d0)
            else
               call mrnorm(d,mmaxim,m,1,iseed)
               do 30 i = 1,m
                  d(i) = 0.5d0 + 0.3d0*dabs(d(i))
 30            continue
            endif
         else
            print *, 'Enter means: '
            read(*,*) (xmu(i), i=1,m)
            print *, 'Enter lower bounds: (A)'
            read(*,*) (a(i), i=1,m)
            print *, 'Enter upper bounds: (B)'
            read(*,*) (b(i), i=1,m)
            if (icase.eq.0 .or. icase.eq.1) then
               print *, 'Enter indep. std. dev.: '
               read(*,*) (d(i), i=1,m)
            endif
            call hadcomp(b,mmaxim,a,mmaxim,tempmx1,mmaxim,m,1,'gt')
            do 70 i = 1,m
               b(i) = a(i) + 1.0d0 + (b(i)-1.0d0-a(i))*tempmx1(i)
 70         continue
            if (icase .ne. 3) then
               call fillm(tempmx1,mmaxim,1,0.0d0)
               call hadcomp(d,mmaxim,tempmx1,mmaxim,
     *                      tempmx1,mmaxim,m,1,'gt')
               do 80 i = 1,m
                  d(i) = 1.0d0 + (d(i)-1.0d0)*tempmx1(i)
 80            continue
            else
               call fillm(d,mmaxim,1,-999.0d0)
            endif
            if ( icase .eq. 0 ) then
               print *, 'Enter Factor Loadings: '
               read(*,*) (xl(i), i=1,m )
            else if (icase .eq. 1) then
               print *, 'Enter Factor Loading Component m: '
               read(*,*) temp
               m=temp
               call fillm(xl,mmaxim,1,1.0d0)
               xl(m) = m
            else if (icase .eq. 3) then
               call fillm(xl,mmaxim,1,-999.0d0)
               print *,'Enter lower triange of OMEGA, ',
     /                 'one row at a time:'
               do 81 i=1,m
                  read(*,*) (w(i,j), j=1,i)
                  print *, (w(i,j), j=1,i)
 81            continue
            endif
         endif
      endif
c
      print *, 'Enter 1 to print FULL    output on the screen'
      print *, '      0 to print SUMMARY output'
      read(*,*) temp
      ifullout=temp
      print *,'***',ifullout
c
c     The following allows for the line in the input file
c     required by the gauss program, giving the number of string
c     parameters in the input file. (See the GAUSS program.)
      print *, 'Hit <Return> when ready:'
      read(*,'(a)') results
      print *, ' '
c
c
      print *, ' '
      print *, 'Basic Name of file to save results in: ',
     /         '(MUST be 2-letters long)'
      results(3:3)=' '
      read(*,'(a)') results
      if(results(3:3).eq.' ') then
         resfil=results(1:2)
         print *, ' '//resfil
         print *, ' '
      else
         print *,'Error!!!  Name must be EXACTLY 2-letters long'
         goto 999
      endif

c--assign a name to output file to each methods.

      cfilCFS  = resfil//'CFS.asd'
      cfilNISE = resfil//'NISE.asd'
      cfilNIST = resfil//'NIST.asd'
      cfilKFS  = resfil//'KFS.asd'
      cfilSDS  = resfil//'SDS.asd'
      cfilGHK  = resfil//'GHK.asd'
      cfilPCF  = resfil//'PCF.asd'
      cfilDCS  = resfil//'DCS.asd'
      cfilARSE = resfil//'ARSE.asd'
      cfilARSR = resfil//'ARSR.asd'
      cfilGSS  = resfil//'GSS.asd'
      cfilAUS  = resfil//'AUS.asd'
      cfilSUS  = resfil//'SUS.asd'

  700 continue
      call fsize(cfilSUS,itemp)
      if(mcrepet.gt.0) then
         if(itemp.gt.0) then
c        *.asd files must NOT exist already
            print '(1x,a/1x,a)',
     *      'Error!!!',
     *      'When MCrepet is Positive, *.asd files must NOT exist.'
            goto 999
         else
            imcrepet=1
         endif
      elseif(mcrepet.lt.0) then
         if(itemp.le.0) then
c        *.asd files MUST exist already
            print '(1x,a/1x,a)',
     *      'Error!!!',
     *      'When MCrepet is Negative, *.asd files MUST exist.'
            goto 999
         else
c        must find actual completed repetition in files
            open(10,file=cfilSUS,
     *      access='direct',form='formatted',recl=130)
            irec=1
  777       continue
ccc            read(10,'(a)',rec=irec,end=778) results
            read(10,'(a)',rec=irec,err=778) results
            irec=irec+1
            goto 777
  778       continue
            close(10)
            imcrepet=irec-2
            print *,'******************************************'
            print *,'*** Will proceed from existing results ***'
            print *,'*** Completed repetitions: ',imcrepet
            print *,'******************************************'
            imcrepet=imcrepet+1
            mcrepet=-mcrepet
         endif
      else
         print *,'The number of MC repetitions cannot be 0'
         print *,'Will carry out 1 MC repetition'
         mcrepet=1
         goto 700
      endif
      write(6,84)
 84   format(' The inputs for this exercise are:',/1x/
     *       5x,'MEANS',10x,'A',14x,'B',9x,'SD',4x,'FACTOR LOADINGS'/
     *       ' -------------------------------------------------------',
     *       ' ------------')
      do 85 i = 1,m
         write(6,'(1x,5g13.5)') xmu(i),a(i),b(i),d(i),xl(i)
 85   continue
      print *,' '
      print *,'Initial Seed: ',iseed1
      do 90 i = 1,m
         s(i) = dsqrt( d(i)**2 + xl(i)**2 )
         pvec(i) =
     /   pheev((b(i)-xmu(i))/s(i)) - pheev((a(i)-xmu(i))/s(i))
 90   continue
      index(1)=1
      call srtrnk(pvec,m,index)
      call matcopy(tempmx1,mmaxim,d  ,mmaxim,1,m,1,1)
      do 701 i=1,m
         d(i)=tempmx1(index(i))
 701  continue
      call matcopy(tempmx1,mmaxim,xl ,mmaxim,1,m,1,1)
      do 702 i=1,m
         xl(i)=tempmx1(index(i))
 702  continue
      call matcopy(tempmx1,mmaxim,a  ,mmaxim,1,m,1,1)
      do 703 i=1,m
         a(i)=tempmx1(index(i))
 703  continue
      call matcopy(tempmx1,mmaxim,b  ,mmaxim,1,m,1,1)
      do 704 i=1,m
         b(i)=tempmx1(index(i))
 704  continue
      call matcopy(tempmx1,mmaxim,xmu,mmaxim,1,m,1,1)
      do 705 i=1,m
         xmu(i)=tempmx1(index(i))
 705  continue
c
      write(6,706)
 706  format(1x/
     *  ' After sorting to put low probability in first component:'/2x/
     *  5x,'Means',10x,'A',12x,'B',11x,'Sd',4x,'Factor loadings'/
     *  ' --------------------------------',
     *  '---------------------------------')
      do 100 i = 1,m
         write(6,'(1x,5g13.5)') xmu(i),a(i),b(i),d(i),xl(i)
 100  continue
      if (icase .ne. 3) then
         do 110 i = 1,m
            tempmx1(i) = d(i)**2
 110     continue
         call fillm(w,mmaxim,m,0.0d0)
         call diagrv(w,mmaxim,m,m,tempmx1,mmaxim)
         call transps(xl,mmaxim,m,1,temp1xm,1)
         call matpv(xl,mmaxim,temp1xm,1,tempmxm,mmaxim,m,1,m)
         call hadoper(w,mmaxim,tempmxm,mmaxim,w,mmaxim,m,m,1)
      else
         do 111 i=2,m
            do 112 j=1,i-1
               w(m-i+1,m-j+1) = w(i,j)
 112        continue
 111     continue
      endif
c
c     calculate chw=CHOL(W)
      call matcopy(tempmxm,mmaxim,w,mmaxim,1,m,1,m)
      call chol(tempmxm,mmaxim,m,chw,mmaxim)
c
c     calculate wi=INV(W)
      call matcopy(wi,mmaxim,w,mmaxim,1,m,1,m)
      call invv(wi,tempmxm,mmaxim,mmaxim,tempmx1,m,
     *          npr,nzr,nnr,condnum,*888)
c
c     calculate detw=DET(W)
      call prodc(tempmx1,mmaxim,m,1,t(1),1)
      detw=t(1)
c
      print *,'  '
      print *,'Implied Variance-Covariance Matrix'
      print *,'----------------------------------'
      call mprint(w,mmaxim,m,m,'omega',6,6)
      print '(/1x,a/1x,a,e15.7,a,e15.7)',
     *      'Var-Cov matrix inverted and its Cholesky factor '//
     *      'calculated.','Condition Number = ', condnum,
     *      '  Determinant = ',detw
c
      do 120 i = 1,20
         qw2(i,1) = qw(i,1)
         qw2(i,2) = qw(i,2)
         qw2(i+20,1) = -qw(i,1)
         qw2(i+20,2) = qw(i,2)
 120  continue
c
ccc---initialize to "Not available"
c
      gresampl = -999.d0
c
      do 125 i = 1,16
         pp(i) = -999.0d0
         dm(i) = -999.0d0
         xlm(i) = -999.0d0
         dl(i) = -999.0d0
         xll(i) = -999.0d0
 125  continue
c
      do 130 i = 1,2
         el(i) = 1.0d0
         rep(i) = 1.d0
 130  continue

c--if omega is arbitrary, skip the quadratures

      if (icase.eq.3) goto 1000
c                          DRAWRAND
c
ccc======================= quadrature 1 ============================ccc
c
      mmaximm1 = mmaxim - 1
      mp1 = m + 1
      mm1 = m - 1
c
      if (icase.eq.1) then
         call diag(w,mmaxim,m,s,mmaxim)
         do 140 i = 1,m
            s(i) = dsqrt(s(i))
            xhs = pheev( (b(i)-xmu(i)) / s(i) )
            rhs = pheev( (a(i)-xmu(i)) / s(i) )
            t(i) = xhs - rhs
 140     continue
         call prodc(t,mmaxim,m,1,tempmx1(1),1)
         pp(1) = tempmx1(1)
         bm = (b(m)-xmu(m)) / s(m)
         am = (a(m)-xmu(m)) / s(m)
         q = nrmpdf(bm,0.0d0,1.0d0) - nrmpdf(am,0.0d0,1.0d0)
         qq = bm*nrmpdf(bm,0.0d0,1.0d0) - am*nrmpdf(am,0.0d0,1.0d0)
         call matcopy(tmm1,mmaximm1,t,mmaxim,1,mm1,1,1)
         call prodc(tmm1,mmaximm1,mm1,1,dlp,1)
         dl(1) = -dlp
         dm(1) = dl(1)*q / s(m)
         xlm(1) = dm(1) / pp(1)
         dl(1) = dl(1)*qq*xl(m) / s(m)**3
         xll(1) = dl(1) / pp(1)
      endif
c
      if (icase.eq.2) then
         pp(1) = 1.0d0 / dble(float(mp1))
         do 150 i = 1,40
            cdf = pheev(qw2(i,1)**mm1)
            pdf = nrmpdf(qw2(i,1),0.0d0,1.0d0)
            vec40(i) = cdf*pdf*qw2(i,2)
            vec40zl(i) = cdf*pdf*qw2(i,1)*qw2(i,2)
 150     continue
         call sumc(vec40,40,40,1,dmp,1)
         dm(1) = -dmp
         xlm(1) = dm(1) / pp(1)
         call sumc(vec40zl,40,40,1,dlp,1)
         dl(1) = dlp
         xll(1) = dl(1) / pp(1)
      endif
c
ccc======================== quadrature 2 =========================ccc
c
      if (m.lt.3) then
         call exaf(mmaxim,m,a,b,xmu,w,xl,pzscal,tzscal,ttzscal)
         pp(2) = pzscal
         dm(2) = tzscal
         dl(2) = ttzscal
         xlm(2) = tzscal / pzscal
         xll(2) = ttzscal / pzscal
      else
         pp(2) = -999.0d0
         dm(2) = -999.0d0
         dl(2) = -999.0d0
         xlm(2) = -999.0d0
         xll(2) = -999.0d0
      endif
c
 1000 print *, '  '
c     So that for a given iseed1, the same iseed will be used
c     at a given imcrepet:
      iseed=iseednew(iseed1,imcrepet)
ccc      print *, 'DRAWRAND'
ccc      print *, 'Current seed: ', iseed
      curseed = iseed
      call mrunif(uuBIG,mmaxim,m,irsimBIG,iseed)
      call rndbase(mmaxim,m,iseed,us,t)
      it = idnint(dble( float(irsimBIG-m-1) / float(m-1) ))
      call antith(mmaxim,m,irsimBIG,it,us,uss)
c-----note: us[1=m,1=m], uss[1=m,1=2*m+4]
      do 161 j = 1, 2*m+4
         do 162 i = 1,m
            uss(i,j) = -dabs(uss(i,j))
 162     continue
 161  continue
      call mrnorm(unBIG,mmaxim,m,irsimBIG,iseed)
      call hadoper(unBIG,mmaxim,unBIG,mmaxim,unBIG2,mmaxim,m,irsimBIG,3)
      call sumc(unBIG2,mmaxim,m,irsimBIG,unzsum,irmaxim)
      call transps(unzsum,irmaxim,irsimBIG,1,unzsumt,1)
      do 171 j = 1,irsimBIG
         do 172 i = 1,m
            usBIG(i,j) = unBIG(i,j) / dsqrt(unzsumt(1,j))
            ussBIG(i,j) = -dabs(usBIG(i,j))
            if (ussBIG(i,j) .eq. 0.0d0) then
               ussBIG(i,j) = ussBIG(i,j) + 1.d-10
            endif
 172     continue
 171  continue

c vah
      if(imcrepet.gt.1) goto 9995
c                            ENDTIM
c
ccc=========================== Exact ==============================ccc
c
ccc   obtain exact results only for first mcrepetitions
ccc   Use as "true" GHK results with 10*irsimBIG repetitions.

      if (icase.eq.3) then
         do 174 i=1,2
            pp(i)=-999.0d0
            dm(i)=-999.0d0
            xlm(i)=-999.0d0
            dl(i)=-999.0d0
            xll(i)=-999.0d0
            rep(i)=-999.0d0
            el(i)=-999.0d0
 174     continue
c
         print *,'*** As TRUE, will use GHK results with 10*irsimBIG ',
     /        'simulations ***'
         print *,'   '
c
ccc------GHK
c
c vah         e1 = hsec1990()
         e1top = hsec1990()
         ir = irsimBIG
         spp = 0
         sdm = 0
         sxlm = 0
         sdl = 0
         sxll = 0
         do 173 ipart=1,10
            print *,'Part ',ipart,' of 10'
            call mrunif(uuBIG,mmaxim,m,ir,iseed)
            print *,'back from mrunif'
            call ghk(mmaxim,irmaxim,
     /               m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
            print *,'back from ghk: p=',p
            call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /               dlloc,xllloc)
            print *,'back from sim: pploc:',pploc
            spp = spp + pploc
            sdm = sdm + dmloc
            sxlm = sxlm + xlmloc
            sdl = sdl + dlloc
            sxll = sxll + xllloc
 173     continue
         pp(3) = spp/10.d0
         dm(3) = sdm/10.d0
         xlm(3) = sxlm/10.d0
         dl(3) = sdl/10.d0
         xll(3) = sxll/10.d0
         rep(3) = 1.0d0
         el(3) = hsec1990() - e1top
         goto 7000
c             HAVEXACT
      endif
c
      e1 = hsec1990()
      do 177 i = 1,40
         vec40(i) = qw2(i,1)
         vec40zl(i) = qw2(i,2)
 177  continue
c
ccc---note: calling ffact here is different from that in gauss
ccc   program.(the gauss call does not include xmu as an argument)
c
      call ffact(mmaxim,m,vec40,a,b,xmu,d,xl,p40,t40,tt40)
      call matptv(p40,40,vec40zl,40,pp(3),1,40,1,1)
      call matptv(t40,40,vec40zl,40,dm(3),1,40,1,1)
      call matptv(tt40,40,vec40zl,40,dl(3),1,40,1,1)
      if(pp(3).gt.0.d0) then
         xlm(3) = dm(3) / pp(3)
         xll(3) = dl(3) / pp(3)
      else
         xlm(3)=-999.d0
         xll(3)=-999.d0
      endif
      if (dl(1) .eq. 0.0d0) then
         dl(1) = dl(3)
      endif
      if (xll(1) .eq. 0.0d0) then
         xll(1) = xll(3)
      endif
c     vah 22 sep 92: added to guard against meanigless derivative zeros
c     returned from ffact:
      if(dm(3).eq.0.d0.and.dl(3).eq.0.d0) then
         dm(3)=-999.d0
         xlm(3)=-999.d0
         dl(3)=-999.d0
         xll(3)=-999.d0
      endif
      write(6,181)
 181  format(1x/'Quadrature1, Quadrature2, and Exact Values'/
     *       ' ---------------------------------------',
     *       '----------------------------------------------')
      write(6,186)
 186  format(20x,'p',8x,'dp/dmu(m)',4x,'dln(p)/dmu(m)',5x,'dp/dl(m)',
     *       5x,'dln(p)/dl(m)'/
     *       ' ---------------------------------------',
     *       '----------------------------------------------')
      write(6,182) pp(1),dm(1),xlm(1),dl(1),xll(1)
 182  format(' Quadrature1',(5g15.5))
      write(6,183) pp(2),dm(2),xlm(2),dl(2),xll(2)
 183  format(' Quadrature2',(5g15.5))
      write(6,184) pp(3),dm(3),xlm(3),dl(3),xll(3)
 184  format(' Exact',6x,(5g15.5))
      write(6,187)
 187  format(' ---------------------------------------',
     *       '----------------------------------------------'/)
      rep(3) = 1.0d0
      e2 = hsec1990()
      ec = e2 - e1
      el(3) = ec
      e1 = e2

C vah HAVEXACT
 7000 continue
ccc---save exact results
c
      call matcopy(pp0,16,pp,16,1,16,1,1)
      call matcopy(dm0,16,dm,16,1,16,1,1)
      call matcopy(xlm0,16,xlm,16,1,16,1,1)
      call matcopy(dl0,16,dl,16,1,16,1,1)
      call matcopy(xll0,16,xll,16,1,16,1,1)
      call matcopy(el0,16,el,16,1,16,1,1)
      call matcopy(rep0,16,rep,16,1,16,1,1)
c
ccc---write "true" values, using exact/quadrature average
c
      do 185 i = 1,16
         aa(i,1) = pp(i)
         aa(i,2) = dm(i)
         aa(i,3) = xlm(i)
         aa(i,4) = dl(i)
         aa(i,5) = xll(i)
         aa(i,6) = rep(i)
         aa(i,7) = el(i)
 185  continue
c
ccc---check each row of aa, if any element is -999.0d0, then delete
c     the row.
c
      call packr(aa,16,3,7,-999.0d0,icount)
      if(icount.gt.0) then
         call meanc(aa,16,icount,7,res,7)
      else
         call fillv(res,7,-999.d0)
      endif
      write(6,188) (res(i),i=1,7)
 188  format(' Will use the following as TRUE values:'/
     *       ' ---------------------------------------',
     *       '----------------------------------------',
     *       '----------------------'/
     *       11x,'p',8x,'dp/dmu(m)',4x,'dln(p)/dmu(m)',5x,'dp/dl(m)',
     *       5x,'dln(p)/dl(m)',5x,'nrepet',5x,'time(sec)'/
     *       ' ---------------------------------------',
     *       '----------------------------------------',
     *       '----------------------'/
     *       (7g15.5)/
     *       ' ---------------------------------------',
     *       '----------------------------------------',
     *       '----------------------'/)
c
ccc   put exact values to the first results of each output files
c
      write(results,'(8(1x,g15.9))') (res(i),i=1,7),curseed
      irec=1
      call reswrite(cfilCFS ,results,irec)
      call reswrite(cfilNISE,results,irec)
      call reswrite(cfilNIST,results,irec)
      call reswrite(cfilKFS ,results,irec)
      call reswrite(cfilSDS ,results,irec)
      call reswrite(cfilGHK ,results,irec)
      call reswrite(cfilPCF ,results,irec)
      call reswrite(cfilDCS ,results,irec)
      call reswrite(cfilARSE,results,irec)
      call reswrite(cfilARSR,results,irec)
      call reswrite(cfilGSS ,results,irec)
      call reswrite(cfilAUS ,results,irec)
      call reswrite(cfilSUS ,results,irec)

      if (iendog .eq. 0) goto 9995
c                             ENDTIM
c
c***********************************************************************
c*    a simple timing test to set number of simulations endogenously   *
c***********************************************************************
c
      print *, 'BEGINTIM'
      print *,' '
c
ccc---CFS
      e1 = hsec1990()
      ir = irsimCFS
      call cfs(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,unBIG,parm,p,hu,hc)
      timeCFS = dmax1(1.0d0, hsec1990()-e1)
      print *, 'Time for ',ir,' simulations using CFS: ',timeCFS
      print *, ' '
c
      write(6,3000)
 3000 format(' ---------------------------------------',
     /       '-------------------------------------------------'/
     /       'Method',24x,'irsimBIG',7x,'# within',7x,'# equivalent',
     /       5x,'Number of'/
     /       45x,'time limit',11x,'to CFS',5x,'repetitions'/
     /       '---------------------------------------',
     /       '-------------------------------------------------')
c
ccc---NISE
      e1 = hsec1990()
      ir = irsimTRY
      call nise(mmaxim,irmaxim,
     /          m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimNSE = dmax1(10.0d0,out3)
      write(6,3001) irsimBIG,out1,out2,irsimNSE
 3001 format(' NIS_exponential:',15x,i5,2g20.5,i9)
c
ccc---NIST
      e1 = hsec1990()
      ir = irsimTRY
      call nist(mmaxim,irmaxim,
     /          m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimNST = dmax1(10.0d0,out3)
      write(6,3002) irsimBIG,out1,out2,irsimNST
 3002 format(' NIS_truncated_normal:',10x,i5,2g20.5,i9)
c
ccc---KFS
      e1 = hsec1990()
      ir = irsimTRY
      do 200 i = 1, m
         vec40(i) = b(i) - a(i)
 200  continue
      do 210 i = m+1, 2*m
         vec40(i) = dsqrt(d(i-m)**2 + xl(i-m)**2)
 210  continue
      val = 0.1d0 * fmin(vec40,2*m)
      call fillm(parm,irmaxim+1,2,val)
      call pkfs(mmaxim,irmaxim,
     /          m,xmu,w,wi,chw,a,b,ir,unBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimKFS = dmax1(10.0d0,out3)
      write(6,3003) irsimBIG,out1,out2,irsimKFS
 3003 format(' KFS:',27x,i5,2g20.5,i9)
c
ccc---SDS
      e1 = hsec1990()
      ir = irsimTRY
      call matcopy(aw,mmaxim,w,mmaxim,1,m,1,m)
c
c-----invv returns eigen values of w in tempmx1 in descending order,
c     so, sorts is used to make it in ascending order.
c
      call invv(aw,tempmxm,mmaxim,mmaxim,tempmx1,m,
     *          npr,nzr,nnr,rcond,*888)
      call sorts(tempmx1,m)
      val = (0.99d0 * tempmx1(1))**2
      call fillm(parm,irmaxim+1,2,val)
      do 220 j = 1,m
         do 221 i = 1,m
            if (i.ne.j) then
               tempmxm(i,j) = w(i,j)
            else
               tempmxm(i,j) = w(i,j) - val
            endif
 221     continue
 220  continue
      call chol(tempmxm,mmaxim,m,tempmxm,mmaxim)
      call sds(mmaxim,irmaxim,
     /         m,xmu,w,wi,tempmxm,a,b,ir,unBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimSDS = dmax1(10.0d0,out3)
      write(6,3004) irsimBIG,out1,out2,irsimSDS
 3004 format(' SDS:',27x,i5,2g20.5,i9)
c
ccc---GHK
      e1 = hsec1990()
      ir = irsimTRY
      call ghk(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimGHK = dmax1(10.0d0,out3)
      write(6,3005) irsimBIG,out1,out2,irsimGHK
 3005 format(' GHK:',27x,i5,2g20.5,i9)
c
ccc---PCF
      e1 = hsec1990()
      ir = irsimTRY
      parm(1,1) = detw
      call pcf(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,ussBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimPCF = dmax1(10.0d0,out3)
      write(6,3006) irsimBIG,out1,out2,irsimPCF
 3006 format(' PCF:',27x,i5,2g20.5,i9)
c
ccc---DCS
      e1 = hsec1990()
      ir = irsimTRY
      parm(1,1) = detw
      call dcs(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,usBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimDCS = dmax1(10.0d0,out3)
      write(6,3007) irsimBIG,out1,out2,irsimDCS
 3007 format(' DCS:',27x,i5,2g20.5,i9)
c
ccc---ARSE
      e1 = hsec1990()
      ir = irsimTRY
      issd = idnint(uuBIG(1,1)*1.0d5)
      do 230 i = 1,ir
         parm(i,1) = idnint( 1.0d5*randu(issd) )
 230  continue
      parm(ir+1,1) = ir
      call arse(mmaxim,irmaxim,
     /          m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimARE = dmax1(10.0d0,out3)
      write(6,3008) irsimBIG,out1,out2,irsimARE,dummy(1)
 3008 format(' ARS_exponential:',15x,i5,2g20.5,i9,
     *'    % not censored:',f10.7)
c
ccc---ARSR
      e1 = hsec1990()
      ir = irsimTRY
      issd = idnint(uuBIG(1,1)*1.0d5)
      do 231 i = 1,ir
         parm(i,1) = idnint( 1.0d5*randu(issd) )
 231  continue
      parm(ir+1,1) = ir
      call arsr(mmaxim,irmaxim,
     /          m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimARR = dmax1(10.0d0,out3)
      write(6,3009) irsimBIG,out1,out2,irsimARR,dummy(1)
 3009 format(' ARR_recursive_truncated_normal:',i5,2g20.5,i9,
     *'    % not censored:',f10.7)
c
ccc---GSS
      e1 = hsec1990()
      ir = irsimTRY
      gresampl = idnint( 5.0d0 + dble(float(ir)/20.0d0) )
cvah  added so as to make sure at least 10 gibbs resamplings are used:
      if(val.lt.10.d0) gresampl=10.d0
      call fillm(parm,irmaxim+1,2,gresampl)
      call gss(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimGSS = dmax1(10.0d0,out3)
      write(6,3010) irsimBIG,out1,out2,irsimGSS
 3010 format(' GSS:',27x,i5,2g20.5,i9)
c
ccc---AUS
      e1 = hsec1990()
      ir = irsimTRY
      do 240 i = 1,ir
         parm(i,1) = idnint( 1.0d5*randu(iseed) )
         parm(i,2) = 1.0d0
 240  continue
      parm(1,2) = idnint( 5.0d0 + dble(rsimBIG/20.0d0) )
      parm(2,2) = idnint( 10.0d0 + dble(rsimBIG/10.0d0) )
      parm(3,2) = idnint( 1.0d6*uuBIG(2,1) )
      call aus(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimAUS = dmax1(10.0d0,out3)
      write(6,3011) irsimBIG,out1,out2,irsimAUS,dummy(1)
 3011 format(' AUS:',27x,i5,2g20.5,i9,'    Max. # of CFS draws:',f10.1)
c
ccc---SUS
      e1 = hsec1990()
      ir = irsimTRY
      do 241 i = 1,ir
         parm(i,1) = idnint( 1.0d5*randu(iseed) )
         parm(i,2) = 1.0d0
 241  continue
      parm(1,2) = idnint( 5.0d0 + dble(rsimBIG/20.0d0) )
      parm(3,2) = idnint( 1.0d6*uuBIG(2,1) )
      call sus(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      etime = dmax1(1.0d0, hsec1990()-e1)
      out1 = hseclim / (etime/dble(float(ir)))
      out2 = dble(float(ir))*timecfs / etime
      out3 = dmin1(rsimBIG,out1,out2)
      irsimSUS = dmax1(10.0d0,out3)
      write(6,3012) irsimBIG,out1,out2,irsimSUS,dummy(1)
 3012 format(' SUS:',27x,i5,2g20.5,i9,'    Max. # of CFS draws:',f10.1/
     *       ' ---------------------------------------',
     *       '------------------------------------------------'/)
c
c
c
      write(6,3013)
 3013 format(' Number of Repetitions'/
     *       ' ------------------------------------------',
     *       '-------------------------------------------------------')
      write(6,3014) irsimCFS,irsimNSE,irsimNST,irsimKFS,irsimSDS,
     *              irsimGHK,irsimPCF,irsimDCS,irsimARE,irsimARR,
     *              irsimGSS,irsimAUS,irsimSUS
 3014 format(1x,
     *       'CFS       =',i5,'    NIS_E     =',i5,'    NIS_T_N   =',i5,
     *       '    KFS       =',i5,'    SDS       =',i5/1x,
     *       'GHK       =',i5,'    PCF       =',i5,'    DCS       =',i5,
     *       '    ARS_E     =',i5,'    ARS_R_T_N =',i5/1x,
     *       'GSS       =',i5,'    AUS       =',i5,'    SUS       =',i5/
     *    1x,'---------------------------------------------',
     *       '---------------------------------------------------')

ccc---insert (5)
c***********************************************************************
c*  a simple timing test to set number of simulations endogenously     *
c***********************************************************************
c
c     ENDTIM
 9995 continue
c

c     BEGINSIM
c
c-----asssign ID number to each of 13 methods
c
      methCFS  = 4
      methNISE = 5
      methNIST = 6
      methKFS  = 7
      methSDS  = 8
      methGHK  = 9
      methPCF  = 10
      methDCS  = 11
      methARSE = 12
      methARSR = 13
      methGSS  = 14
      methAUS  = 15
      methSUS  = 16
c
c
c-----initialize matrices to exact results
c
      call matcopy(pp,16,pp0,16,1,16,1,1)
      call matcopy(dm,16,dm0,16,1,16,1,1)
      call matcopy(xlm,16,xlm0,16,1,16,1,1)
      call matcopy(dl,16,dl0,16,1,16,1,1)
      call matcopy(xll,16,xll0,16,1,16,1,1)
      call matcopy(rep,16,rep0,16,1,16,1,1)
      call matcopy(el,16,el0,16,1,16,1,1)
c
c-----start of simulators
      print '(/a/a,i5,a,f13.1/a/)',
     *' ********************************',
     *' * Start of MC Repetition ',imcrepet,' *  Current seed:',curseed,
     *' ********************************'
      irec=imcrepet+1
c
ccc---CFS
      e1 = hsec1990()
      ir = irsimCFS
      call cfs(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,unBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
c-----notice that first three elements of the following vectors are
c     for the exact results
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5)','CFS:    P =',pploc,'     Time =',ec

c
      rep(methCFS) = ir
      el(methCFS) = ec
      pp(methCFS) = pploc
      dm(methCFS) = dmloc
      xlm(methCFS) = xlmloc
      dl(methCFS) = dlloc
      xll(methCFS) = xllloc
c
c-----time per replication
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methCFS),dm(methCFS),xlm(methCFS),dl(methCFS),
     /      xll(methCFS),rep(methCFS),el(methCFS),curseed
      call reswrite(cfilCFS,results,irec)
c
ccc---NISE
      e1 = hsec1990()
      ir = irsimNSE
      call nise(mmaxim,irmaxim,
     /          m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5)','NISE:   P =',pploc,'     Time =',ec
c
      rep(methNISE) = ir
      el(methNISE) = ec
      pp(methNISE) = pploc
      dm(methNISE) = dmloc
      xlm(methNISE) = xlmloc
      dl(methNISE) = dlloc
      xll(methNISE) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methNISE),dm(methNISE),xlm(methNISE),dl(methNISE),
     /      xll(methNISE),rep(methNISE),el(methNISE),curseed
      call reswrite(cfilNISE,results,irec)
c
ccc---NIST
      e1 = hsec1990()
      ir = irsimNST
      call nist(mmaxim,irmaxim,
     /          m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5)','NIST:   P =',pploc,'     Time =',ec
c
      rep(methNIST) = ir
      el(methNIST) = ec
      pp(methNIST) = pploc
      dm(methNIST) = dmloc
      xlm(methNIST) = xlmloc
      dl(methNIST) = dlloc
      xll(methNIST) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methNIST),dm(methNIST),xlm(methNIST),dl(methNIST),
     /      xll(methNIST),rep(methNIST),el(methNIST),curseed
      call reswrite(cfilNIST,results,irec)
c
ccc---KFS
      e1 = hsec1990()
      ir = irsimKFS
      do 250 i = 1, m
 250     vec40(i) = b(i) - a(i)
      continue
      do 251 i = m+1, 2*m
 251     vec40(i) = dsqrt(d(i-m)**2 + xl(i-m)**2)
      continue
      val = 0.1d0 * fmin(vec40,2*m)
      call fillm(parm,irmaxim+1,2,val)
      call pkfs(mmaxim,irmaxim,
     /          m,xmu,w,wi,chw,a,b,ir,unBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5,a,g13.5)','KFS:    P =',pploc,
     *      '     Time =',ec,'   Window =',parm(methKFS,1)
c
      rep(methKFS) = ir
      el(methKFS) = ec
      pp(methKFS) = pploc
      dm(methKFS) = dmloc
      xlm(methKFS) = xlmloc
      dl(methKFS) = dlloc
      xll(methKFS) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methKFS),dm(methKFS),xlm(methKFS),dl(methKFS),
     /      xll(methKFS),rep(methKFS),el(methKFS),curseed
      call reswrite(cfilKFS,results,irec)
c
ccc---SDS
      e1 = hsec1990()
      ir = irsimSDS
      call matcopy(aw,mmaxim,w,mmaxim,1,m,1,m)
      call invv(aw,tempmxm,mmaxim,mmaxim,tempmx1,m,
     *          npr,nzr,nnr,rcond,*888)
      call sorts(tempmx1,m)
      val = (0.99d0 * tempmx1(1))**2
      call fillm(parm,irmaxim+1,2,val)
      do 260 j = 1,m
         do 261 i = 1,m
            if (i.ne.j) then
               tempmxm(i,j) = w(i,j)
            else
               tempmxm(i,j) = w(i,j) - val
            endif
 261     continue
 260  continue
      call chol(tempmxm,mmaxim,m,tempmxm,mmaxim)
      call sds(mmaxim,irmaxim,
     /         m,xmu,w,wi,tempmxm,a,b,ir,unBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5,a,g13.5)','SDS:    P =',pploc,
     *      '     Time =',ec,'   Lambda =',parm(methsds,1)
      rep(methSDS) = ir
      el(methSDS) = ec
      pp(methSDS) = pploc
      dm(methSDS) = dmloc
      xlm(methSDS) = xlmloc
      dl(methSDS) = dlloc
      xll(methSDS) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methSDS),dm(methSDS),xlm(methSDS),dl(methSDS),
     /      xll(methSDS),rep(methSDS),el(methSDS),curseed
      call reswrite(cfilSDS,results,irec)
c
ccc---GHK
      e1 = hsec1990()
      ir = irsimGHK
      call ghk(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5)','KFS:    P =',pploc,'     Time =',ec
c
      rep(methGHK) = ir
      el(methGHK) = ec
      pp(methGHK) = pploc
      dm(methGHK) = dmloc
      xlm(methGHK) = xlmloc
      dl(methGHK) = dlloc
      xll(methGHK) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methGHK),dm(methGHK),xlm(methGHK),dl(methGHK),
     /      xll(methGHK),rep(methGHK),el(methGHK),curseed
      call reswrite(cfilGHK,results,irec)
c
ccc---PCF
      e1 = hsec1990()
      ir = irsimPCF
      parm(1,1) = detw
      call pcf(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,ussBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5)','GHK:    P =',pploc,'     Time =',ec
c
      rep(methPCF) = ir
      el(methPCF) = ec
      pp(methPCF) = pploc
      dm(methPCF) = dmloc
      xlm(methPCF) = xlmloc
      dl(methPCF) = dlloc
      xll(methPCF) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methPCF),dm(methPCF),xlm(methPCF),dl(methPCF),
     /      xll(methPCF),rep(methPCF),el(methPCF),curseed
      call reswrite(cfilPCF,results,irec)
c
ccc---DCS
      e1 = hsec1990()
      ir = irsimDCS
      parm(1,1) = detw
      call dcs(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,usBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5)','DCS:    P =',pploc,'     Time =',ec
c
      rep(methDCS) = ir
      el(methDCS) = ec
      pp(methDCS) = pploc
      dm(methDCS) = dmloc
      xlm(methDCS) = xlmloc
      dl(methDCS) = dlloc
      xll(methDCS) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methDCS),dm(methDCS),xlm(methDCS),dl(methDCS),
     /      xll(methDCS),rep(methDCS),el(methDCS),curseed
      call reswrite(cfilDCS,results,irec)
c
ccc---ARSE
      e1 = hsec1990()
      ir = irsimARE
      issd = idnint(uuBIG(1,1)*1.0d5)
      do 232 i = 1,ir
         parm(i,1) = idnint( 1.0d5*randu(issd) )
 232  continue
      parm(ir+1,1) = ir
      call arse(mmaxim,irmaxim,
     /          m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(4(a,g13.5))',
     *      ' ARSE:   P =',pploc,
     *      '     Time =',ec,'   Parm =',parm(methARSE,1),
     *      '  % not censored:',dummy(1)
c
      rep(methARSE) = ir
      el(methARSE) = ec
      pp(methARSE) = pploc
      dm(methARSE) = dmloc
      xlm(methARSE) = xlmloc
      dl(methARSE) = dlloc
      xll(methARSE) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methARSE),dm(methARSE),xlm(methARSE),dl(methARSE),
     /      xll(methARSE),rep(methARSE),el(methARSE),curseed
      call reswrite(cfilARSE,results,irec)
c
ccc---ARSR
      e1 = hsec1990()
      ir = irsimARR
      issd = idnint(uuBIG(1,1)*1.0d5)
      do 233 i = 1,ir
         parm(i,1) = idnint( 1.0d5*randu(issd) )
 233  continue
      parm(ir+1,1) = ir
      call arsr(mmaxim,irmaxim,
     /          m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(4(a,g13.5))',
     *      ' ARSR:   P =',pploc,
     *      '     Time =',ec,'   Parm =',parm(methARSR,1),
     *      '  % not censored:',dummy(1)
c
      rep(methARSR) = ir
      el(methARSR) = ec
      pp(methARSR) = pploc
      dm(methARSR) = dmloc
      xlm(methARSR) = xlmloc
      dl(methARSR) = dlloc
      xll(methARSR) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methARSR),dm(methARSR),xlm(methARSR),dl(methARSR),
     /      xll(methARSR),rep(methARSR),el(methARSR),curseed
      call reswrite(cfilARSR,results,irec)
c
ccc---GSS
      e1 = hsec1990()
      ir = irsimGSS
      if (gresampl.lt.0.d0) then
c     if < 0, this is not an endogenous timing run, so define gresampl:
         gresampl = idnint( 5.0d0 + dble(float(ir)/20.0d0) )
cvah     added so as to make sure at least 10 gibbs resamplings are used:
         if(val.lt.10.d0) gresampl=10.d0
      endif
      call fillm(parm,irmaxim+1,2,gresampl)
      call gss(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5,a,g13.5)','GSS:    P =',pploc,
     *      '     Time =',ec,'   Parm =',parm(methGSS,1)
c
      rep(methGSS) = ir
      el(methGSS) = ec
      pp(methGSS) = pploc
      dm(methGSS) = dmloc
      xlm(methGSS) = xlmloc
      dl(methGSS) = dlloc
      xll(methGSS) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methGSS),dm(methGSS),xlm(methGSS),dl(methGSS),
     /      xll(methGSS),rep(methGSS),el(methGSS),curseed
      call reswrite(cfilGSS,results,irec)
c
ccc---AUS
      e1 = hsec1990()
      ir = irsimAUS
      do 242 i = 1,ir
         parm(i,1) = idnint( 1.0d5*randu(iseed) )
         parm(i,2) = 1.0d0
 242  continue
      parm(1,2) = idnint( 5.0d0 + dble(rsimBIG/20.0d0) )
      parm(2,2) = idnint( 10.0d0 + dble(rsimBIG/10.0d0) )
      parm(3,2) = idnint( 1.0d6*uuBIG(2,1) )
      call aus(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5,a,g13.5,a,f10.1)','AUS:    P =',pploc,
     *      '     Time =',ec,'   Parm =',parm(methAUS,1),
     *      '  Max. # of CFS draws:',dummy(1)
c
      rep(methAUS) = ir
      el(methAUS) = ec
      pp(methAUS) = pploc
      dm(methAUS) = dmloc
      xlm(methAUS) = xlmloc
      dl(methAUS) = dlloc
      xll(methAUS) = xllloc
c
      ec = ec / dble(float(ir))
c
      write(results,'(8(1x,g15.9))')
     /      pp(methAUS),dm(methAUS),xlm(methAUS),dl(methAUS),
     /      xll(methAUS),rep(methAUS),el(methAUS),curseed
      call reswrite(cfilAUS,results,irec)
c
ccc---SUS
      e1 = hsec1990()
      ir = irsimSUS
      do 243 i = 1,ir
         parm(i,1) = idnint( 1.0d5*randu(iseed) )
         parm(i,2) = 1.0d0
 243  continue
      parm(1,2) = idnint( 5.0d0 + dble(rsimBIG/20.0d0) )
      parm(3,2) = idnint( 1.0d6*uuBIG(2,1) )
      call sus(mmaxim,irmaxim,
     /         m,xmu,w,wi,chw,a,b,ir,uuBIG,parm,p,hu,hc)
      call sim(mmaxim,m,xl,e1,p,hu,hc,ec,pploc,dmloc,xlmloc,
     /         dlloc,xllloc)
c
      if (ifullout.gt.0)
     *print '(1x,a,g13.5,a,g13.5,a,g13.5,a,f10.1)','SUS:    P =',pploc,
     *      '     Time =',ec,'   Parm =',parm(methSUS,1),
     *      '  Max. # of CFS draws:',dummy(1)
c
      rep(methSUS) = ir
      el(methSUS) = ec
      pp(methSUS) = pploc
      dm(methSUS) = dmloc
      xlm(methSUS) = xlmloc
      dl(methSUS) = dlloc
      xll(methSUS) = xllloc
      ec = ec / float(ir)
c
      write(results,'(8(1x,g15.9))')
     /      pp(methSUS),dm(methSUS),xlm(methSUS),dl(methSUS),
     /      xll(methSUS),rep(methSUS),el(methSUS),curseed
      call reswrite(cfilSUS,results,irec)
c
      if (ifullout.gt.0) then
c        print out full results
         write(6,1320)
 1320    format(1x/' Collected Results'/
     *       1x,'-----------------------------------------------',
     *          '------------------------------------------------'/
     *       1x,'           p        dp/dmu(m)   dln(p)/dmu(m) ',
     *          '  dp/dl(m)    dln(p)/dl(m)    nrepet    time(sec)'/
     *       1x,'-----------------------------------------------',
     *          '------------------------------------------------')
         print '(a,7g13.5)',' Quad.1:',pp(1),dm(1),xlm(1),dl(1),xll(1),
     *                        rep(1),el(1)
         print '(a,7g13.5)',' Quad.2:',pp(2),dm(2),xlm(2),dl(2),xll(2),
     *                        rep(2),el(2)
         print '(a,7g13.5)',' Exact: ',pp(3),dm(3),xlm(3),dl(3),xll(3),
     *                        rep(3),el(3)
         print '(a,7g13.5)',' CFS:   ',pp(methcfs),dm(methcfs),
     *                        xlm(methcfs),dl(methcfs),xll(methcfs),
     *                        rep(methcfs),el(methcfs)
         print '(a,7g13.5)',' NISE:  ',pp(methnise),dm(methnise),
     *                        xlm(methnise),dl(methnise),xll(methnise),
     *                        rep(methnise),el(methnise)
         print '(a,7g13.5)',' NIST:  ',pp(methnist),dm(methnist),
     *                        xlm(methnist),dl(methnist),xll(methnist),
     *                        rep(methnist),el(methnist)
         print '(a,7g13.5)',' KFS:   ',pp(methkfs),dm(methkfs),
     *                        xlm(methkfs),dl(methkfs),xll(methkfs),
     *                        rep(methkfs),el(methkfs)
         print '(a,7g13.5)',' SDS:   ',pp(methsds),dm(methsds),
     *                        xlm(methsds),dl(methsds),xll(methsds),
     *                        rep(methsds),el(methsds)
         print '(a,7g13.5)',' GHK:   ',pp(methghk),dm(methghk),
     *                        xlm(methghk),dl(methghk),xll(methghk),
     *                        rep(methghk),el(methghk)
         print '(a,7g13.5)',' PCF:   ',pp(methpcf),dm(methpcf),
     *                        xlm(methpcf),dl(methpcf),xll(methpcf),
     *                        rep(methpcf),el(methpcf)
         print '(a,7g13.5)',' DCS:   ',pp(methdcs),dm(methdcs),
     *                        xlm(methdcs),dl(methdcs),xll(methdcs),
     *                        rep(methdcs),el(methdcs)
         print '(a,7g13.5)',' ARSE:  ',pp(metharse),dm(metharse),
     *                        xlm(metharse),dl(metharse),xll(metharse),
     *                        rep(metharse),el(metharse)
         print '(a,7g13.5)',' ARSR:  ',pp(metharsr),dm(metharsr),
     *                        xlm(metharsr),dl(metharsr),xll(metharsr),
     *                        rep(metharsr),el(metharsr)
         print '(a,7g13.5)',' GSS:   ',pp(methgss),dm(methgss),
     *                        xlm(methgss),dl(methgss),xll(methgss),
     *                        rep(methgss),el(methgss)
         print '(a,7g13.5)',' AUS:   ',pp(methaus),dm(methaus),
     *                        xlm(methaus),dl(methaus),xll(methaus),
     *                        rep(methaus),el(methaus)
         print '(a,7g13.5)',' SUS:   ',pp(methsus),dm(methsus),
     *                        xlm(methsus),dl(methsus),xll(methsus),
     *                        rep(methsus),el(methsus)
         write(6,3021)
 3021    format(' -----------------------------------------------',
     *          '------------------------------------------------'/)
      endif
c
c
c     undocumented feature:
c     Will stop at the end of a Monte Carlo round, if stop.//resfil
c     exists and has MORE than 0 bytes -- careful!.
      call fsize('stop.'//resfil,itemp)
      if (itemp.gt.0) then
         print *,'---------'
         print *,'Will stop because file  stop.'//resfil//'  exists.'
         print *,'---------'
         goto 999
      endif
      if (imcrepet.lt.mcrepet) then
         imcrepet = imcrepet + 1
         goto 1000
c             DRAWRAND
      else
         goto 999
c             FINISH
      endif
 888  continue
      print *, 'error in inversion routine'
c     end of the program
 999  continue
c     FINISH
      end
