      program recipe
c     RECIPE, a FORTRAN program for REGression Confidence Intervals
c        on PErcentiles (i.e., one-sided tolernace limits) for a
c        large class of mixed models.
c
c     Mark Vangel, NIST, Gaithersburg, MD 20899
c     Version 1.0, April 1995
c
c        RECIPE is a program for determining one-sided tolerance 
c     limits from mixed models involving two variance components.
c     The fixed part of the model can be arbitrary, and the random 
c     part need not be balanced in any way, except that it is 
c     necessary that each observation have the same variance, 
c     and that this variance equal the sum of the two components 
c     of variance (e.g., no interactions between fixed and random 
c     effects are permitted).  The theory is documented in
c
c     Vangel, M. G. (1995),``One-Sided $\beta$-Content 
c        Tolerance Limits for Mixed Models'', Proceedings
c        of the Section on Physical and Engineering Sciences,
c        American Statistical Association
c
c     Vangel, M. G. (1995), ``ANOVA Estimates for `Partially-
c        Balanced' Mixed Models'', draft manuscript.
c
c     Vangel, M. G. (1995), ``Design Allowables From Regression
c        Models Using Data From Several Batches'', Proceedings
c        of the 12th Symposium on Composite Materials: Testing
c        and Design, American Society for Testing and Materials,
c        STP ???, to appear.
c
c     The user's manual for this program is
c
c     Vangel, M. G. (1995), ``A User's manual for RECIPE'',
c        draft manuscript.
c      
c
c     Required subroutines:
c
c     LINPACK --  DSVDC, Singular value decomposition.
c     BLAS    --  Various level 1, 2, and 3 routines.  BLAS routines
c                 have been used as much as possible, for both
c                 readibility and efficiency.
c     SLATEC  --  DSORT, sort algorithm,
c                 RNOR, normal pseudo-random variables
c                   (and various low-level utility routines)
c     ZERO1   --  Root finder for noncentral-t inverse.  ACM algorithm
c                 which uses bisection and rational interpolation is
c                 used.
c     TNC     --  Lenth's algorithm for noncentral-t cdf, along with
c                 its required subroutines.
c
      parameter (mtot=100, mlvl=100, mbch=25, mpar=10, mpts=100)
c
c     MLIM   --   Set mlim = max(mpar*mpar, mbch*mbch, mpar+mbch)
c                 (Some PC compilers don't allow expressions in
c                  parameter statements; you might want to try the
c                  statement ...
c
c     parameter(mlim=max(mpar*mpar, mbch*mbch, mpar+mbch))
c
c             ... and see if your compiler accepts it)
c
      parameter(mlim=625)
c
c     Set the values in the parameter statement to suit the size of
c     your applications, or the memory constraints of your computer:
c
c     mtot   --  Maximum number of response values (Ys)
c     mlvl   --  Maximum number of fixed effect levels (distinct
c                sets of Xs)
c     mbch   --  Maximum number of batches (random effect levels)
c     mpar   --  Maximum number of fixed parameters (regression
c                coefficients)
c     mpts   --  Maximum number of points at which tolerance limit
c                is to be calculated
c
c
c     2017/09: Make the following tweaks (Alan Heckert):
c
c              1. Do not require the input file to have a ".dat"
c                 extension.
c              2. Replace the use of assigned GO TO as this is
c                 considered a deleted feature in current versions
c                 of the Fortran standard.
c
      implicit double precision (a-h, o-z)
      logical      satt
      character*80 rec
      character*80 flenme
      character*80 flenm2
c
      dimension x(1000), xx(3500), xpts(1000),
     $          xtx(100), xtxi(100),
     $          u1(1000), s1(10), v1(100),
     $          u2(3500), s2(35),
     $          v2(1225),
     $          xn(2500), h(10000),
     $          ip(100), iq(100), y(100), coef(10),
     $          eta0(100), eta1(100), tlm0(100), tlm1(100),
     $          xm(100), t(100),
     $          wk1(3500), wk2(625), wk3(135) 
c
c     -- If your compiler will allow it, PLEASE substitute the
c        following dimension statement for the above dimension
c        statement.  If you must use the above dimension statemnt
c        and you change any of the parameter values refer to the
c        following statement and CAREFULLY adjust the above
c        dimension statement accordingly.  I apologize for the
c        inelegance of this approach, but portability to
c        bush-league PC compilers makes it necessary.
c
c     dimension x(mlvl*mpar), xx(mtot*(mpar+mbch)), xpts(mpts*mpar),
c    $          xtx(mpar*mpar), xtxi(mpar*mpar),
c    $          u1(mtot*mpar), s1(mpar), v1(mpar*mpar),
c    $          u2(mtot*(mpar+mbch)), s2(mpar+mbch),
c    $          v2((mpar+mbch)*(mpar+mbch)),
c    $          xn(mlvl*mbch), h(mlvl*mlvl),
c    $          ip(mtot), iq(mtot), y(mtot), coef(mpar),
c    $          eta0(mpts), eta1(mpts), tlm0(mpts), tlm1(mpts),
c    $          xm(mpts), t(mpts),
c    $          wk1(mtot*(mpar+mbch)), wk2(mlim), wk3(mtot+mpar+mbch),
c    $          simvec(mpivsm), crit(mpts)
      data one, zero/1.d0, 0.d0/, in1 /10/, in2 /20/
c
c    -- Open the data file
      write (*,*) ' Filename (if no file name extension given, ',
     1            '.dat assumed) ?'
      read  (*,fmt='(a80)') flenme
      do 1 i=80, 1, -1
       if (flenme(i:i) .ne. ' ') then
         lstnbk = i
         go to 2
       end if
1     continue
2     continue
c
c     check if file name extension given
c
      iper=0
      do 11 i=lstnbk, 1, -1
       if (flenme(i:i) .eq. '.') then
         iper = i
         go to 12
       end if
11    continue
12    continue
c
      if (iper.gt.0) then
         open (unit=in1, file=flenme(1:lstnbk), iostat=lstat,
     $         status='old')
      else if (iper.eq.0) then
         iend = lstnbk +4
c
c        -- Open the data file, and the critical value file 
c           (if a critical value file exists)
         flenme((lstnbk+1):iend) = '.dat'
         open (unit=in1, file=flenme(1:iend), iostat=lstat,
     $         status='old')
         if (lstat .ne. 0) then
            flenme((lstnbk+1):iend) = '.DAT'
            open (unit=in1, file=flenme(1:iend), iostat=lstat,
     $            status='old')
         end if
      endif
c
      if (lstat .ne. 0) then
        write (*,*) 'RECIPE : Data file not found'
        stop
      end if
c
c     -- Now check for the critical value file
      satt = .TRUE.
      flenm2 = 'n'
      flenm2(1:lstnbk) = flenme(1:lstnbk)
      flenm2((lstnbk+1):iend) = '.crt'
      open (unit=in2, file=flenm2(1:iend), iostat=lstat,
     $    status='old')
      if (lstat .eq. 0) satt = .FALSE.
      if (lstat .ne. 0) then
        flenm2((lstnbk+1):iend) = '.CRT'
        open (unit=in2, file=flenm2(1:iend), iostat=lstat,
     $    status='old')
         if (lstat .eq. 0) satt = .FALSE.
      end if
      
c
c     -- Print out header, and indicate whether Satterthwaite
c        approximation is being used.
      write (*,*)
      write (*,*)
     $' RECIPE : One-Sided Random-Effect Regression Tolerance Limits'
      write (*,*) ' (Version 1.0, April 1995) '
      write (*,*)
      if (.NOT. satt) then
         write (*,*) 
     $          ' *** Simulated pivot critical values from file ',
     $           flenm2(1:iend),' will be used.'
      else
         write (*,*) 
     $          ' *** Simulated pivot critical value file ',
     $                 flenm2(1:iend),' not found.'
         write (*,*) '    Satterthwaite approximation will be used.'
      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 (*,*) ' RECIPE : Error -- ntot exceeds mtot '
            ier = 1
         end if
         if (nlvl .gt. mlvl) then
            write (*,*) ' RECIPE : Error -- nlvl exceeds mlvl '
            ier = 1
         end if
         if (nbch .gt. mbch) then
            write (*,*) ' RECIPE : Error -- nbch exceeds mbch '
            ier = 1
         end if
         if (npar .gt. mpar) then
            write (*,*) ' RECIPE : Error -- npar exceeds mpar '
            ier = 1
         end if
         if (npts .gt. mpts) then
            write (*,*) ' RECIPE : 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 data y(i), along with the corresponding
c          fixed level indicators ip(i) and 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), y(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 second call
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)
      call regdat (npar, ntot, nbch, npts, xpts, y, coef, u1, s1,
     $             v1, u2, tlm0, tlm1, eta0, eta1, wk1, xm, t)
c

      write (*,*)
     $  '       Probability    Confidence   Regression  Tolerance',
     $  ' Limit'
      do 40 i=1, npts
        write (*,1000) cont, conf, xm(i), t(i)
1000    format (2f15.2,2F15.6)
 40   continue
      go to 100
 99   continue
c
c     -- Error in input file
      write (*,*) ' RECIPE : Error reading input file '
      write (*,*) ' record number = ',irec
      write (*,*) ' io-status     = ',lstat
      write (*,*) ' bad record    = ',rec
 100  continue
      stop
      end
