      program  simpvt
c
c     Mark Vangel, April, 1995
c
c        Program `simpvt' simulates the critical values to be
c     used for tolerance limit factors in the limit of very 
c     large between-batch variability (equivalently, very small
c     within-batch variability).  These critical values should
c     be used instead of the Satterthwaite approximation when
c     the data are very unbalanced.
c
c        This program takes as input the same input file that
c     is used by the `recipe' tolerance limit program.  If
c     the input filename for `recipe' is `name.dat' then 
c     `simpvt' will produce an output file `name.crt' which
c     contains the simulated critical values.  This critical
c     value file will be read by the `recipe' program and used
c     instead of the Satterthwaite approximation.
c
      parameter (mtot=1000, mlvl=100, mbch=25, mpar=10, mpts=100,
     $           mpivsm=1000000)
      implicit double precision (a-h, o-z)
      character*80  rec
      character*30 flenme
      character*30 flenm2
      dimension x(mlvl*mpar), xx(mtot*(mpar+mbch)), xpts(mpts*mpar),
     $          u1(mtot*mpar), s1(mpar), v1(mpar*mpar),
     $          ip(mlvl), iq(mtot), wk1(mtot), wk2(mtot),
     $          vals(mpivsm), w(mpar)
      data in, iout/10, 20/
c
      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=in, file=flenme(1:iend), status='old', iostat=lstat)
      if (lstat .ne. 0) then
         flenme(1:iend) = flenme(1:lstnbk) // '.DAT'
         open(unit=in, file=flenme(1:iend), status='old', iostat=lstat)
         lcflag = 2
      end if
      if (lstat .ne. 0) then
         write (*,*) ' SIMPVT : 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=iout, file=flenm2(1:iend), status='new',iostat=lstat)
      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 (in,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 (*,*) ' SIMPVT : Error -- ntot exceeds mtot '
            ier = 1
         end if
         if (nlvl .gt. mlvl) then
            write (*,*) ' SIMPVT : Error -- nlvl exceeds mlvl '
            ier = 1
         end if
         if (nbch .gt. mbch) then
            write (*,*) ' SIMPVT : Error -- nbch exceeds mbch '
            ier = 1
         end if
         if (npar .gt. mpar) then
            write (*,*) ' SIMPVT : Error -- npar exceeds mpar '
            ier = 1
         end if
         if (npts .gt. mpts) then
            write (*,*) ' SIMPVT : 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 (in,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 data, along with the corresponding
c          fixed level indicators ip(i) and batch indicators iq(i).
      do 20 i=1, ntot
 5    continue
         read (in,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 (in,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      -- Build full data matrix from unique rows
      do 40 i=1, ntot
         do 50 j=1, npar
            xx((j-1)*ntot+i) = x((j-1)*nlvl+ip(i))
 50      continue
 40   continue
c
c    --  SVD of design matrix xx 
      ijob = 21
      ldu  = ntot
      ldv  = npar
      tol  = 1.d-7
      call dsvdc (xx , ntot, ntot, npar, s1, wk1, u1, ldu, v1, ldv,
     $            wk2, ijob, info)
c
c    -- Rank (xx)
      irk = 0
      do 60 i=1, npar
         if (abs(s1(i)) .lt. tol) go to 70
         irk = irk +1
 60   continue
 70   continue
c
c     -- Now, simulate the ratio quantiles for each of `npts'
c        vectors of independent variables
 7    continue
         write (*,*) ' Number of simulation replicates ?'
         read  (*,*)   nrep
         if (nrep .gt. mpivsm) then
            write (*,*) 
     $      " Number of replicates should not exceed 'mpivsm'"
            write (*,*) 
     $      "     (Either enter smaller value, or change the "
            write (*,*) 
     $      "      parameter 'mpivsm' and recompile the program)"
         go to 7
      end if
      write (*,*) ' Integer seed ?'
      read  (*,*)   iseed
      zcont = ppnd16(cont, ifault)
      write (*,*)
      write (*,*)
     $' SIMPVT : One-Sided Random-Effect Regression Tolerance Limits'
      write (*,*) ' (Version 1.0,  April 1995) '
      write (*,*)
      write (*,*) '     Simulated critical values '
      write (*,*) 
      write (*,*) ' Number of values     = ',npts
      write (*,*) ' Number of replicates = ',nrep
      write (*,*) ' Seed                 = ',iseed
      write (*,*) ' Input file           = ',flenme(1:iend)
      write (*,*) ' Output file          = ',flenm2(1:iend)
      write (*,*)
      write (iout,*) ' Data: ',flenme(1:iend),
     $              ' Num. Reps.: ',nrep,
     $              ' Seed: ',iseed
      do 100 i=1, npts
         z = rnor(iseed)
         call dcopy (npar, xpts(i), npts, w, 1)
         call simrat
     $     (u1,s1,v1,iq,w,nbch,ntot,npar,nrep,irk,zcont,conf,
     $      wk1,wk2,vals,quant)
         write (*,*) quant
         write (iout,*) quant
 100  continue
c
      go to 200
 99   continue
c
c     -- Error in input file
      write (*,*) ' SIMPVT : Error reading input file '
      write (*,*) ' record number = ',irec
      write (*,*) ' io-status     = ',lstat
      write (*,*) ' bad record    = ',rec
200   continue
      stop
      end







