program simcov
      parameter (mtot=500, mlvl=100, mbch=25, mpar=10, mpts=100)
c
      parameter(mlim=max(mpar*mpar, mbch*mbch, mpar+mbch))
c
      implicit double precision (a-h, o-z)
      logical satt, debug
      character*80  rec
      character*30 flenme
      character*30 flenm2
      character*1  pch
      real err, bch, rnor
c
      dimension x(mlvl*mpar), xx(mtot*(mpar+mbch)), xpts(mpts*mpar),
     $          xtx(mpar*mpar), xtxi(mpar*mpar),
     $          u1(mtot*mpar), s1(mpar), v1(mpar*mpar),
     $          u2(mtot*(mpar+mbch)), s2(mpar+mbch),
     $          v2((mpar+mbch)*(mpar+mbch)),
     $          xn(mlvl*mbch), h(mlvl*mlvl),
     $          ip(mtot), iq(mtot), y(mtot), coef(mpar),
     $          eta0(mpts), eta1(mpts), tlm0(mpts), tlm1(mpts),
     $          xm(mpts), t(mpts),
     $          wk1(mtot*(mpar+mbch)), wk2(mlim), wk3(mtot+mpar+mbch),
     $          bch(mbch), err(mtot),
     $          cov(mpts), xmu(mpts)
      data one, zero/1.d0, 0.d0/, in1 /10/, in2/20/, debug/.FALSE./
c
c    -- Open the data file
      write (*,*) " Filename (without '.dat' extension) ?"
      read  (*,fmt='(a20)') flenme
      do 1 i=20, 1, -1
       if (flenme(i:i) .ne. ' ') then
         lstnbk = i
         go to 2
       end if
1     continue
2     continue
      iend = lstnbk +4
c
c     -- Open the data file, and the critical value file
c        (if a critical value file exists)
c        Either a .dat or a .DAT extension is allowed. 
c        Program quits if data file cannot be opened.
      lcflag = 1
      flenme(1:iend) = flenme(1:lstnbk) // '.dat'
      open (unit=in1, file=flenme(1:iend), status='old', iostat=lstat)
      if (lstat .ne. 0) then
         flenme(1:iend) = flenme(1:lstnbk) // '.DAT'
         open(unit=in1, file=flenme(1:iend), status='old', iostat=lstat)
         lcflag = 2
      end if
      if (lstat .ne. 0) then
         write (*,*) ' SIMCOV : Error opening data file; iostat = ',
     $        lstat
         stop
      end if
      if (lcflag .eq. 1) then
         flenm2(1:iend) = flenme(1:lstnbk) // '.crt'
      else
         flenm2(1:iend) = flenme(1:lstnbk) // '.CRT'
      end if
      open (unit=in2, file=flenm2(1:iend), status='old',iostat=lstat)
c
c     -- Print out header, and indicate whether Satterthwaite
c        approximation is being used.
      write (*,*)
      write (*,*)
     $' SIMCOV : One-Sided Random-Effect Regression Tolerance Limits'
      write (*,*) ' (Version 1.0, April 1995) '
      write (*,*)
      if (lstat .eq. 0) then
         write (*,*)
     $          ' *** Simulated pivot critical values from file ',
     $           flenm2(1:iend),' will be used.'
         satt = .FALSE.
      else
         write (*,*)
     $          ' *** Simulated pivot critical value file ',
     $                 flenm2(1:iend),' not found.'
         write (*,*) '    Satterthwaite approximation will be used.'
         satt = .TRUE.
      end if
      irec = 0
c
c    -- Read in parameters
c         ntot  -- Total number of observations
c         nlvl  -- Number of levels of fixed effect
c         nbch  -- Number of batches
c         npar  -- Number of regression parameters
c         npts  -- Number of points on regression surface at which
c                  the tolerance limit is to be evaluated
c         cont  -- Probability content for tolerance limit
c                  (=.99 for A-Basis; .90 for B-Basis)
c         conf  -- Confidence for lower tolerance limit
c                  (=.95 for both A- and B-Basis)
c
3     continue
         read (in1,fmt='(a80)') rec
         irec = irec +1
         if (rec(1:1) .eq. '#') go to 3
         read (unit=rec,fmt=*,iostat=lstat) ntot, nlvl, nbch, npar,
     $                                  npts, cont, conf
         if (lstat .ne. 0) go to 99
c
c     -- Make sure that the problem size doesn't exceed the
c        dataspace of the program
c
         ier = 0
         if (ntot .gt. mtot) then
            write (*,*) ' SIMCOV : Error -- ntot exceeds mtot '
            ier = 1
         end if
         if (nlvl .gt. mlvl) then
            write (*,*) ' SIMCOV : Error -- nlvl exceeds mlvl '
            ier = 1
         end if
         if (nbch .gt. mbch) then
            write (*,*) ' SIMCOV : Error -- nbch exceeds mbch '
            ier = 1
         end if
         if (npar .gt. mpar) then
            write (*,*) ' SIMCOV : Error -- npar exceeds mpar '
            ier = 1
         end if
         if (npts .gt. mpts) then
            write (*,*) ' SIMCOV : Error -- npts exceeds mpts  '
            ier = 1
         end if
         if (ier .ne. 0) then
            write (*,*)
            write (*,*) ' The values in the first non-comment '
            write (*,*) ' line of your input file indicate that '
            write (*,*) ' performing this analysis may exceed'
            write (*,*) ' the array bounds of this program. To '
            write (*,*) ' increase these bounds, modify the '
            write (*,*) ' PARAMETER and DIMENSION statements in '
            write (*,*) ' the main program source file, then '
            write (*,*) ' create a new executable program '
            write (*,*)
            stop
         end if
c
c     -- Now read in the levels of the fixed effects.  All
c          two dimensional arrays used in this program are
c          packed into one-dimensional arrays by columns.
      do 10 i=1, nlvl
 4    continue
         read (in1,fmt='(a80)') rec
         irec = irec +1
         if (rec(1:1) .eq. '#') go to 4
         read (unit=rec,fmt=*, iostat=lstat)
     $        (x((j-1)*nlvl+i), j=1, npar)
         if (lstat .ne. 0) go to 99
 10   continue
c
c     -- Next read in the fixed level indicators ip(i)
c         and the batch indicators iq(i).
      do 20 i=1, ntot
 5    continue
         read (in1,fmt='(a80)') rec
         irec = irec +1
         if (rec(1:1) .eq. '#') go to 5
         read (unit=rec,fmt=*,iostat=lstat) ip(i), iq(i)
         if (lstat .ne. 0) go to 99
 20   continue
c
c     -- Read in the points on the regression surface at which
c          the lower tolerance limits are to be calculated
      do 30 i=1, npts
 6    continue
         read (in1,fmt='(a)') rec
         irec = irec +1
         if (rec(1:1) .eq. '#') go to 6
         read (unit=rec,fmt=*,iostat=lstat)
     $        (xpts((j-1)*npts+i), j=1, npar)
         if (lstat .ne. 0) go to 99
 30   continue
c
c     -- Calculated the tolerance limits.  The first call is
c          to a routine which does all of the calculations
c          which do not depend on y.  The succeeding calls
c          performs the remaining calculations, which do
c          depend on y, and returns the result in the
c          vector 'xm' of point estimates, and the
c          vector 't' of tolerance limits.
      call regini (nlvl, npar, ntot, nbch, npts, x, xpts, ip, iq,
     $             cont, conf, xx, xtx, xtxi, xn, h, u1,
     $             s1, v1, u2, s2, v2, tlm0, tlm1, eta0, eta1,
     $             satt, in2, wk1, wk2, wk3)
c
c     -- Get input from user for simulation
      write (*,*)
      write (*,*) ' Number of simulation replicates ?'
      read  (*,*)   nsim
      write (*,*) ' Integer seed ?'
      read  (*,*)   iseed
      write (*,*) ' Number of values for intraclass correlation ?'
      read  (*,*)   nrho
      drho = dfloat(1)/(nrho -1)
 7    write (*,*) '1= use same random numbers for each rho'
      write (*,*) '0= use different random numbers for each rho ?'
      read  (*,*) itype
      if (itype .ne. 0 .and. itype .ne. 1) go to 7
      write (*,*) 
      write (*,*)
      write (*,*) '       rho       confidence'
      write (*,*)
c
c     -- Initialize the random number generator
      z = rnor (iseed)
c
c     -- Loop over intraclass correlation values
      do 70 irho = 1, nrho
      rho = (irho -1)*drho
c
c     -- If itype=1, reset seed for each rho
      if (itype .eq. 1) z = rnor (iseed)
      pct = -ppnd16(cont, ifault)
      sdb = sqrt(rho)
      sdw = sqrt(1-rho)
      if (nbch .eq. 1) then
          sdb = 0
          sdw = 1
      end if
      do 49 idx=1, npts
         cov(idx) = 0
         xmu(idx) = 0
 49   continue
c
c     -- Loop over simulation replicates
      do 50 idx=1, nsim
c
c     -- Generate pseudo-random dataset
      do 511 i=1, nbch
        bch (i) = rnor (0)
 511  continue
      do 521 i=1, ntot
        err(i) = rnor(0)
 521  continue
      do 51 i=1, ntot
         y(i) = bch(iq(i)) *sdb +err(i)*sdw
 51   continue
c
c     -- Calculate the tolerance limit for this dataset
      call regdat (npar, ntot, nbch, npts, xpts, y, coef, u1, s1,
     $             v1, u2, tlm0, tlm1, eta0, eta1, wk1, xm, t)
 3000 format (11f15.6)
c
      do 40 i=1, npts
       if (t(i) .lt. pct) cov(i) = cov(i) +1
       if (debug) then
         if (t(i) .gt. pct) then
            pch = '+'
        else
           pch = ' '
        end if
        write (1,1000) xm(i), t(i), tlm1(i), eta1(i)
1000    format (4f12.6)
        end if
 40   continue
 50   continue
      do 60 i=1, npts
         write (*,2000) rho, cov(i)/nsim
 2000    format(3f12.4)
 60   continue
 70   continue
      go to 100
 99   continue
c
c     -- Error in input file
      write (*,*) ' SIMCOV : Error reading input file '
      write (*,*) ' record number = ',irec
      write (*,*) ' io-status     = ',lstat
      write (*,*) ' bad record    = ',rec
 100  continue
      stop
      end