c
c   glm.f
c
c   This program can be the basis for a linear models
c   package.  Currently it creates no plots, checks no
c   residuals, and deals only with fixed effects models.
c   On the other hand, it can handle unbalanced data, nested
c   effects, analysis of covariance, and arbitrary interactions.
c   (The X matrix must be provided in a file by the user).
c   It also performs power and sample size calculations.
c
c
c   Steve Verrill
c
c
c   10/9/84
c   2/8/85
c   2/16/85
c   12/5/90 (and earlier) - modified to run on a Sun
c   5/91 - more modifications
c   7/99 - some cleanup
c

c   The parameter and dimension statements in the main, ccheck, 
c   hyptest, estim, and power subroutines must be in accord.

      parameter (nrdx=400,numpar=75,nrdc=numpar,numest=numpar)

      implicit double precision (a-h,o-z)

      dimension x(nrdx,numpar),y(nrdx)
      dimension c(nrdc,numest)
      dimension est(numest),cli(numest),cui(numest),
     x          clj(numest),cuj(numest)
      dimension cb1(5),sigma(5),cb(10)

      character*80 dataf,results

      common /hec2/ resss,iresdof
      common /ataic/ atai(numest,numest)
      common /alphac/ alpha(3)
      common /nrowc/ nrow
      common /threshmc/ threshm

c
c   Determine machine precision
c

      xmult = .99d0
      delta = 1.0d0
 1    delta = xmult*delta
      if ((1.0d0 + delta) .gt. 1.0d0) go to 1

      delta = delta/xmult

c      threshm = sqrt(delta)
      threshm = delta*10.0d0

      write(6,2)
 2    format(/,1x,'What name do you want for your results file?',/)

      read(5,3) results
 3    format(a80)

      open(unit=9,file=results,status='unknown')

      write(6,4)
      write(9,4)
 4    format(/,1x,'What is the name of the data file?',/)

      read(5,3) dataf
      write(9,3) dataf

      open(unit=3,file=dataf,status='old')
      open(unit=4,file='runf',status='unknown')
      open(unit=7,file='runfc',status='unknown')
      open(unit=8,file='csandas',status='unknown')

      iend = 0

c   The input that is used to produce the current run
c   is written to runfc.
c   If, later, you want to reproduce the run with
c   a different data file,
c   copy runfc to runf and answer yes to the
c   following question:

      write(6,10)
      write(9,10)
 10   format(/,1x,'Do you want to read the run control cards',
     x' from the file runf?',/,1x,'0 - no   1 - yes',/)

      read(5,*) irunf
      write(9,*) irunf
      if (irunf .eq. 1) go to 5000

      write(6,1030)
      write(9,1030)
 1030 format(/,1x,'How many rows are there in the X matrix',
     x' (400 or fewer)?',/)

      read(5,*) nrow
      write(9,*) nrow
      write(7,*) nrow

      write(6,1040)
      write(9,1040)
 1040 format(/,1x,'How many columns are there in the X matrix',
     x' (75 or fewer)?',/)

      read(5,*) ncol
      write(9,*) ncol
      write(7,*) ncol

c
c   Read the X matrix
c

      do 1100 i = 1,nrow

         read(3,*) (x(i,j),j = 1,ncol)
         write(9,*) (x(i,j),j = 1,ncol)

 1100 continue


      write(6,1152)
      write(9,1152)
 1152 format(/,1x,'Do you want to evaluate the design of an experiment',
     x', or',/,1x,'do you want to analyze the results from an ',
     x'experiment?',//,1x,'Design - 0   Analysis - 1',/)

      read(5,*) idfl
      write(9,*) idfl
      write(7,*) idfl
      if (idfl .eq. 0) go to 3004

c
c   Read the y vector
c

      do 1160 i = 1,nrow

         read(3,*) y(i)
         write(9,*) y(i)

 1160 continue

c
c   Hypothesis testing and Estimation
c

 1170 write(6,1200)
      write(9,1200)
 1200 format(/,1x,"Do you wish to test the hypothesis",
     x" c1'b = c2'b = ... = cq'b = 0 ",/,1x,"for some q and some c's?",
     x"     0 - no   1 - yes",/)

      read(5,*) ihypf
      write(9,*) ihypf
      write(7,*) ihypf

      if (ihypf .eq. 0) then

         if (iend .eq. 1) stop
         
         iend = 1
         go to 2000

      endif

      write(6,1210)
      write(9,1210)
 1210 format(/,1x,'What is q?',/)

      read(5,*) nq
      write(9,*) nq
      write(7,*) nq

      write(6,1220)
      write(9,1220)
 1220 format(/,1x,'What are the c primes?',/)

      write(8,1225)
 1225 format(//,1x,"The c's are (as rows)",/)

      do 1300 iii = 1,nq

         read(5,*) (c(k,iii), k = 1,ncol)
         write(9,*) (c(k,iii), k = 1,ncol)
         write(7,*) (c(k,iii), k = 1,ncol)
         write(8,*) (c(k,iii), k = 1,ncol)

 1300 continue

      call ccheck(c,ncol,nq,notind,nrdc,numest)

      iend = 0

      if (notind .eq. 1) go to 1170

c
c   Perform the hypothesis tests
c

      call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
     x             notestf,0)

      iend = 0

      if (notestf .eq. 1) go to 1170
      
c
c   Write the results of the hypothesis tests
c

      write(6,1402) hypss
      write(9,1402) hypss
 1402 format(///,1x,'The hypothesis sum of squares is ',2x,e12.4,/)

      write(6,1410) ihypdof
      write(9,1410) ihypdof
 1410 format(/,1x,'The hypothesis dof is',2x,i3,/)

      write(6,1420) resss
      write(9,1420) resss
 1420 format(//,1x,'The residual sum of squares is ',2x,e12.4,/)

      write(6,1430) iresdof
      write(9,1430) iresdof
 1430 format(/,1x,'The residual dof is',2x,i3,/)

      xnum = ihypdof
      xden = iresdof
      f = (hypss/xnum)/(resss/xden)
      write(6,1440) f
      write(9,1440) f
 1440 format(//,1x,'The f statistic is',2x,e12.4,/)

c
c   cdff is a routine from DCDFLIB.  It performs calculations
c   associated with the F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

      call cdff(1,p,q,f,xnum,xden,istatus,bound)

      if (istatus .ne. 0) then

         write(6,1442) istatus
 1442    format(/,1x,'Error 1442: The istatus value from an',
     x   ' F routine was ',i2,'.',/)

         stop

      endif

      pvalue = q

      write(6,1445) pvalue
      write(9,1445) pvalue
 1445 format(/,1x,'The p value is',2x,e12.4,/)


c
c   Estimation
c

 2000 write(6,2200)
      write(9,2200)
 2200 format(/,1x,"Do you wish to estimate linear combinations",
     x" of the parameters c1'b, c2'b, ... , cq'b",/,1x,
     x"for some q and for some c's?     0 - no   1 - yes",/)

      read(5,*) iestf
      write(9,*) iestf
      write(7,*) iestf

      if (iestf .eq. 0) then

         if (iend .eq. 1) stop

         iend = 1
         go to 1170

      endif

      isame = 0

      if (ihypf .eq. 0) go to 2209

      write(6,2205)
      write(9,2205)
 2205 format(/,1x,"Are the c's the same as those used in the",
     x" hypothesis testing?   0 - no   1 - yes",/)

      read(5,*) isame
      write(9,*) isame
      write(7,*) isame

      if (isame .eq. 1) go to 2500

 2209 write(6,2210)
      write(9,2210)
 2210 format(/,1x,'What is q?',/)

      read(5,*) nq
      write(9,*) nq
      write(7,*) nq

      write(6,2220)
      write(9,2220)
 2220 format(/,1x,'What are the c primes?',/)

      write(8,2225)
 2225 format(//,1x,"The c's are (as rows)",/)

      do 2300 iii = 1,nq

         read(5,*) (c(k,iii), k = 1,ncol)
         write(9,*) (c(k,iii), k = 1,ncol)
         write(7,*) (c(k,iii), k = 1,ncol)
         write(8,*) (c(k,iii), k = 1,ncol)

 2300 continue

 2500 write(6,2510)
      write(9,2510)
 2510 format(/,1x,'What are the individual and the joint alphas?',/,
     x1x,'(confidence level = (1 - alpha)*100 %)',/)

      read(5,*) alphai,alphaj
      write(9,*) alphai,alphaj
      write(7,*) alphai,alphaj

c
c   Do the estimations
c

      call estim(x,y,nrdx,nrow,ncol,c,nrdc,nq,alphai,alphaj,
     x           est,cli,cui,clj,cuj,notestf)

      iend = 0

      if (notestf .eq. 1) go to 2000

c
c   Write the results of the estimations
c

      conli = (1.0 - alphai)*100.0
      conlj = (1.0 - alphaj)*100.0

      write(6,2590) conli
      write(9,2590) conli
 2590 format(///,1x,'The confidence level for the individual',
     x' confidence intervals is ',f5.1,' percent',/)

      write(6,2595) conlj
      write(9,2595) conlj
 2595 format(/,1x,'The confidence level for the joint',
     x' confidence intervals is ',f5.1,' percent',/)

      write(6,2600)
      write(9,2600)
 2600 format(/,1x,'The estimates, individual confidence intervals,',
     x' and joint confidence intervals are',/)

      write(6,2602)
      write(9,2602)
 2602 format(/,12x,'estimate',14x,'individual',21x,'joint',/)

      do 2620 i = 1,nq

         write(6,2610) i,est(i),cli(i),cui(i),clj(i),cuj(i)
         write(9,2610) i,est(i),cli(i),cui(i),clj(i),cuj(i)
 2610    format(/,1x,"c",i2,"'b",4x,e12.4,4x,e12.4,1x,e12.4,
     x   4x,e12.4,1x,e12.4)

 2620 continue

      go to 1170

c
c   Experimental design material --- estimability and power
c

 3004 alpha(1) = .10
      alpha(2) = .05
      alpha(3) = .01

 3005 write(6,3010)
      write(9,3010)
 3010 format(/,1x,"Do you wish to check whether some c'bs are estimable"
     x,"?   no - 0   yes - 1",/)

      read(5,*) iestable
      write(9,*) iestable
      write(7,*) iestable

      if (iestable .eq. 0) then

         if (iend .eq. 1) stop

         iend = 1
         go to 4000

      endif

      write(6,3020)
      write(9,3020)
 3020 format(/,1x,'How many c primes are there?',/)

      read(5,*) nq
      write(9,*) nq
      write(7,*) nq

      write(6,3030)
      write(9,3030)
 3030 format(/,1x,'What are the c primes?',/)

      do 3050 iii = 1,nq

         read(5,*) (c(k,iii), k = 1,ncol)
         write(9,*) (c(k,iii), k = 1,ncol)
         write(7,*) (c(k,iii), k = 1,ncol)

 3050 continue

      call ccheck(c,ncol,nq,notind,nrdc,numest)

      iend = 0

      if (notind .eq. 1) go to 3005

c
c   Check the estimability of the c's
c

      call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
     x             notestf,1)

      iend = 0

      if (notestf .eq. 1) go to 3005

c
c   Perform a power calculation
c

 4000 write(6,4010)
      write(9,4010)
 4010 format(/,1x,'Do you want to perform a power calculation?',
     x'   no - 0   yes - 1',/)

      read(5,*) ipowfl
      write(9,*) ipowfl
      write(7,*) ipowfl

      if (ipowfl .eq. 0) then

         if (iend .eq. 1) stop

         iend = 1
         go to 3005

      endif

      ialready = 0

      if (iestable .eq. 0) go to 4040

      write(6,4030)
      write(9,4030)
 4030 format(/,1x,'Are the c primes the same as those for which',
     x' you just performed estimability calculations?',/,1x,
     x'no - 0   yes - 1',/)

      read(5,*) ialready
      write(9,*) ialready
      write(7,*) ialready

      if (ialready .eq. 1) go to 4100

 4040 write(6,4050)
      write(9,4050)
 4050 format(/,1x,"What is the q in the hypothesis c1'b=...=cq'b",
     x"=0 ?",/)

      read(5,*) nq
      write(9,*) nq
      write(7,*) nq

      write(6,4060)
      write(9,4060)
 4060 format(/,1x,'What are the c primes?',/)

      do 4070 iii = 1,nq

         read(5,*) (c(k,iii), k = 1,ncol)
         write(9,*) (c(k,iii), k = 1,ncol)
         write(7,*) (c(k,iii), k = 1,ncol)

 4070 continue

      call ccheck(c,ncol,nq,notind,nrdc,numest)

      iend = 0

      if (notind .eq. 1) go to 4000

 4100 if (nq .gt. 1) go to 4109

      write(6,4101)
      write(9,4101)
 4101 format(/,1x,"How many c'b values (5 or fewer) do you want ",
     x"in the power table?",/)

      read(5,*) numcb
      write(9,*) numcb
      write(7,*) numcb

      write(6,4102)
      write(9,4102)
 4102 format(/,1x,"What are the c'b values?",/)

      read(5,*) (cb1(i),i = 1,numcb)
      write(9,*) (cb1(i),i = 1,numcb)
      write(7,*) (cb1(i),i = 1,numcb)

      go to 4125

 4109 write(6,4110)
      write(9,4110)
 4110 format(/,1x,"What are the values of the c'bs under the ",
     x"alternative hypothesis?",/)

      read(5,*) (cb(i),i = 1,nq)
      write(9,*) (cb(i),i = 1,nq)
      write(7,*) (cb(i),i = 1,nq)

 4125 write(6,4130)
      write(9,4130)
 4130 format(/,1x,'How many values (5 or fewer) of sigma ',
     x'do you want in the power table?',/)

      read(5,*) numsig
      write(9,*) numsig
      write(7,*) numsig

      write(6,4140)
      write(9,4140)
 4140 format(/,1x,'What are the sigma values?',/)

      read(5,*) (sigma(i),i = 1,numsig)
      write(9,*) (sigma(i),i = 1,numsig)
      write(7,*) (sigma(i),i = 1,numsig)

      write(6,4160)
      write(9,4160)
 4160 format(/,1x,'What power do you wish to achieve?',/)

      read(5,*) powerl
      write(9,*) powerl
      write(7,*) powerl

      if (ialready .eq. 1) go to 4170

      call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
     x             notestf,2)

      iend = 0

      if (notestf .eq. 1) go to 4000

 4170 call power(nq,cb1,numcb,sigma,numsig,cb,powerl,iresdof)

      iend = 0

      go to 3005

c
c   End of experimental design material
c

c
c   Reading the run control cards from a file
c

 5000 write(6,5030)
      write(9,5030)
 5030 format(/,1x,'How many rows are there in the X matrix',
     x' (400 or fewer)?',/)

      read(4,*) nrow
      write(6,*) nrow
      write(9,*) nrow

      write(6,5040)
      write(9,5040)
 5040 format(/,1x,'How many columns are there in the X matrix',
     x' (75 or fewer)?',/)

      read(4,*) ncol
      write(6,*) ncol
      write(9,*) ncol

c
c   Read the X matrix
c

      do 5100 i = 1,nrow

         read(3,*) (x(i,j),j = 1,ncol)
         write(9,*) (x(i,j),j = 1,ncol)

 5100 continue

      iend = 0

      write(6,5152)
      write(9,5152)
 5152 format(/,1x,'Do you want to evaluate the design of an experiment',
     x', or',/,1x,'do you want to analyze the results from an ',
     x'experiment?',//,1x,'Design - 0   Analysis - 1',/)

      read(4,*) idfl
      write(6,*) idfl
      write(9,*) idfl
      if (idfl .eq. 0) go to 7004

c
c   Read the y vector
c

      do 5160 i = 1,nrow

         read(3,*) y(i)
         write(9,*) y(i)

 5160 continue

c
c   Hypothesis testing and Estimation
c

 5170 write(6,5200)
      write(9,5200)
 5200 format(/,1x,"Do you wish to test the hypothesis",
     x" c1'b = c2'b = ... = cq'b = 0 ",/,1x,"for some q and some c's?",
     x"     0 - no   1 - yes",/)

      read(4,*) ihypf
      write(6,*) ihypf
      write(9,*) ihypf

      if (ihypf .eq. 0) then

         if (iend .eq. 1) stop

         iend = 1
         go to 6000

      endif

      write(6,5210)
      write(9,5210)
 5210 format(/,1x,'What is q?',/)

      read(4,*) nq
      write(6,*) nq
      write(9,*) nq

      write(6,5220)
      write(9,5220)
 5220 format(/,1x,'What are the c primes?',/)

      write(8,5225)
 5225 format(//,1x,"The c's are (as rows)",/)

      do 5300 iii = 1,nq

         read(4,*) (c(k,iii), k = 1,ncol)
         write(6,*) (c(k,iii), k = 1,ncol)
         write(9,*) (c(k,iii), k = 1,ncol)
         write(8,*) (c(k,iii), k = 1,ncol)

 5300 continue

      call ccheck(c,ncol,nq,notind,nrdc,numest)

      iend = 0

      if (notind .eq. 1) go to 5170

c
c   Perform the hypothesis tests
c

      call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
     x             notestf,0)

      iend = 0

      if (notestf .eq. 1) go to 5170

c
c   Write the results of the hypothesis tests
c

      write(6,1402) hypss
      write(9,1402) hypss


      write(6,1410) ihypdof
      write(9,1410) ihypdof


      write(6,1420) resss
      write(9,1420) resss


      write(6,1430) iresdof
      write(9,1430) iresdof


      xnum = ihypdof
      xden = iresdof
      f = (hypss/xnum)/(resss/xden)
      write(6,1440) f
      write(9,1440) f


c
c   cdff is a routine from DCDFLIB.  It performs calculations
c   associated with the F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

      call cdff(1,p,q,f,xnum,xden,istatus,bound)

      if (istatus .ne. 0) then

         write(6,5442) istatus
 5442    format(/,1x,'Error 5442: The istatus value from an',
     x   ' F routine was ',i2,'.',/)

         stop

      endif

      pvalue = q

      write(6,1445) pvalue
      write(9,1445) pvalue



c
c   Estimation
c

 6000 write(6,6200)
      write(9,6200)
 6200 format(/,1x,"Do you wish to estimate linear combinations",
     x" of the parameters c1'b, c2'b, ... , cq'b",/,1x,
     x"for some q and for some c's?     0 - no   1 - yes",/)

      read(4,*) iestf
      write(6,*) iestf
      write(9,*) iestf

      if (iestf .eq. 0) then

         if (iend .eq. 1) stop

         iend = 1
         go to 5170

      endif

      isame = 0

      if (ihypf .eq. 0) go to 6209

      write(6,6205)
      write(9,6205)
 6205 format(/,1x,"Are the c's the same as those used in the",
     x" hypothesis testing?   0 - no   1 - yes",/)

      read(4,*) isame
      write(6,*) isame
      write(9,*) isame

      if (isame .eq. 1) go to 6500

 6209 write(6,6210)
      write(9,6210)
 6210 format(/,1x,'What is q?',/)

      read(4,*) nq
      write(6,*) nq
      write(9,*) nq

      write(6,6220)
      write(9,6220)
 6220 format(/,1x,'What are the c primes?',/)

      write(8,6225)
 6225 format(//,1x,"The c's are (as rows)",/)

      do 6300 iii = 1,nq

         read(4,*) (c(k,iii), k = 1,ncol)
         write(6,*) (c(k,iii), k = 1,ncol)
         write(9,*) (c(k,iii), k = 1,ncol)
         write(8,*) (c(k,iii), k = 1,ncol)

 6300 continue

 6500 write(6,6510)
      write(9,6510)
 6510 format(/,1x,'What are the individual and the joint alphas?',/,
     x1x,'(confidence level = (1 - alpha)*100 %)',/)

      read(4,*) alphai,alphaj
      write(6,*) alphai,alphaj
      write(9,*) alphai,alphaj

c
c   Do the estimations
c

      call estim(x,y,nrdx,nrow,ncol,c,nrdc,nq,alphai,alphaj,
     x           est,cli,cui,clj,cuj,notestf)

      iend = 0

      if (notestf .eq. 1) go to 6000

c
c   Write the results of the estimations
c

      conli = (1.0 - alphai)*100.0
      conlj = (1.0 - alphaj)*100.0

      write(6,6590) conli
      write(9,6590) conli
 6590 format(///,1x,'The confidence level for the individual',
     x' confidence intervals is ',f5.1,' percent',/)

      write(6,6595) conlj
      write(9,6595) conlj
 6595 format(/,1x,'The confidence level for the joint',
     x' confidence intervals is ',f5.1,' percent',/)

      write(6,6600)
      write(9,6600)
 6600 format(/,1x,'The estimates, individual confidence intervals,',
     x' and joint confidence intervals are',/)

      write(6,6602)
      write(9,6602)
 6602 format(/,12x,'estimate',14x,'individual',21x,'joint',/)

      do 6620 i = 1,nq

         write(6,6610) i,est(i),cli(i),cui(i),clj(i),cuj(i)
         write(9,6610) i,est(i),cli(i),cui(i),clj(i),cuj(i)
 6610    format(/,1x,"c",i2,"'b",4x,e12.4,4x,e12.4,1x,e12.4,
     x   4x,e12.4,1x,e12.4)

 6620 continue

      go to 5170

c
c   Experimental design material --- estimability and power
c

 7004 alpha(1) = .10
      alpha(2) = .05
      alpha(3) = .01

 7005 write(6,7010)
      write(9,7010)
 7010 format(/,1x,"Do you wish to check whether some c'bs are estimable"
     x,"?   no - 0   yes - 1",/)

      read(4,*) iestable
      write(6,*) iestable
      write(9,*) iestable

      if (iestable .eq. 0) then

         if (iend .eq. 1) stop

         iend = 1
         go to 8000

      endif

      write(6,7020)
      write(9,7020)
 7020 format(/,1x,'How many c primes are there?',/)

      read(4,*) nq
      write(6,*) nq
      write(9,*) nq

      write(6,7030)
      write(9,7030)
 7030 format(/,1x,'What are the c primes?',/)

      do 7050 iii = 1,nq

         read(4,*) (c(k,iii), k = 1,ncol)
         write(6,*) (c(k,iii), k = 1,ncol)
         write(9,*) (c(k,iii), k = 1,ncol)

 7050 continue

      call ccheck(c,ncol,nq,notind,nrdc,numest)

      iend = 0

      if (notind .eq. 1) go to 7005

c
c   Check the estimability of the c's
c

      call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
     x             notestf,1)

      iend = 0

      if (notestf .eq. 1) go to 7005

c
c   Perform a power calculation
c

 8000 write(6,8010)
      write(9,8010)
 8010 format(/,1x,'Do you want to perform a power calculation?',
     x'   no - 0   yes - 1',/)

      read(4,*) ipowfl
      write(6,*) ipowfl
      write(9,*) ipowfl

      if (ipowfl .eq. 0) then

         if (iend .eq. 1) stop

         iend = 1
         go to 7005

      endif

      ialready = 0

      if (iestable .eq. 0) go to 8040

      write(6,8030)
      write(9,8030)
 8030 format(/,1x,'Are the c primes the same as those for which',
     x' you just performed estimability calculations?',/,1x,
     x'no - 0   yes - 1',/)

      read(4,*) ialready
      write(6,*) ialready
      write(9,*) ialready

      if (ialready .eq. 1) go to 8100

 8040 write(6,8050)
      write(9,8050)
 8050 format(/,1x,"What is the q in the hypothesis c1'b=...=cq'b",
     x"=0 ?",/)

      read(4,*) nq
      write(6,*) nq
      write(9,*) nq

      write(6,8060)
      write(9,8060)
 8060 format(/,1x,'What are the c primes?',/)

      do 8070 iii = 1,nq

         read(4,*) (c(k,iii), k = 1,ncol)
         write(6,*) (c(k,iii), k = 1,ncol)
         write(9,*) (c(k,iii), k = 1,ncol)

 8070    continue

      call ccheck(c,ncol,nq,notind,nrdc,numest)

      iend = 0

      if (notind .eq. 1) go to 8000

 8100 if (nq .gt. 1) go to 8109

      write(6,8101)
      write(9,8101)
 8101 format(/,1x,"How many c'b values (5 or fewer) do you want ",
     x"in the power table?",/)

      read(4,*) numcb
      write(6,*) numcb
      write(9,*) numcb

      write(6,8102)
      write(9,8102)
 8102 format(/,1x,"What are the c'b values?",/)

      read(4,*) (cb1(i),i = 1,numcb)
      write(6,*) (cb1(i),i = 1,numcb)
      write(9,*) (cb1(i),i = 1,numcb)

      go to 8125

 8109 write(6,8110)
      write(9,8110)
 8110 format(/,1x,"What are the values of the c'bs under the ",
     x"alternative hypothesis?",/)

      read(4,*) (cb(i),i = 1,nq)
      write(6,*) (cb(i),i = 1,nq)
      write(9,*) (cb(i),i = 1,nq)

 8125 write(6,8130)
      write(9,8130)
 8130 format(/,1x,'How many values (5 or fewer) of sigma ',
     x'do you want in the power table?',/)

      read(4,*) numsig
      write(6,*) numsig
      write(9,*) numsig

      write(6,8140)
      write(9,8140)
 8140 format(/,1x,'What are the sigma values?',/)

      read(4,*) (sigma(i),i = 1,numsig)
      write(6,*) (sigma(i),i = 1,numsig)
      write(9,*) (sigma(i),i = 1,numsig)

      write(6,8160)
      write(9,8160)
 8160 format(/,1x,'What power do you wish to achieve?',/)

      read(4,*) powerl
      write(6,*) powerl
      write(9,*) powerl

      if (ialready .eq. 1) go to 8170

      call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
     x             notestf,2)

      iend = 0

      if (notestf .eq. 1) go to 8000

 8170 call power(nq,cb1,numcb,sigma,numsig,cb,powerl,iresdof)

      iend = 0

      go to 7005

c
c   End of experimental design material
c

      end




c
c   Matrix addition
c

      subroutine matadd(a,b,c,nrow,ncol,nrdim)

      implicit double precision (a-h,o-z)

      dimension a(nrdim,1),b(nrdim,1),c(nrdim,1)

      do 20 i = 1,nrow

         do 10 j = 1,ncol

            c(i,j) = a(i,j) + b(i,j)

 10      continue

 20   continue

      return

      end




c
c   Matrix subtraction
c

      subroutine matsub(a,b,c,nrow,ncol,nrdim)

      implicit double precision (a-h,o-z)

      dimension a(nrdim,1),b(nrdim,1),c(nrdim,1)

      do 20 i = 1,nrow

         do 10 j = 1,ncol

            c(i,j) = a(i,j) - b(i,j)

 10      continue

 20   continue

      return

      end




c
c   Matrix multiplication
c

      subroutine matmult(a,nra,nca,b,nrb,ncb,c,nrda,nrdb)

      implicit double precision (a-h,o-z)

      dimension a(nrda,1),b(nrdb,1),c(nrda,1)

      if (nca .ne. nrb) then

         write(6,10)
         write(9,10)
10       format(1x,'Matrix multiplication error:',/,
     x   1x,'The a column and b row dimensions do not match.')
         stop

      endif

      do 50 j = 1,ncb

         do 40 i = 1,nra

            c(i,j) = 0.0

            do 30 k = 1,nca

               c(i,j) = c(i,j) + a(i,k)*b(k,j)

 30         continue

 40      continue

 50   continue

      return

      end




c
c   Matrix transposition
c

      subroutine mattran(a,nra,nca,at,nrda,ncda)

      implicit double precision (a-h,o-z)

      dimension a(nrda,1),at(ncda,1)

      do 20 i = 1,nra

         do 10 j = 1,nca

            at(j,i) = a(i,j)

 10      continue

 20   continue

      return

      end




      function dot(x,y,n)

      implicit double precision (a-h,o-z)

      dimension x(n),y(n)

      dot = 0.0

      do 10 i = 1,n

         dot = dot + x(i)*y(i)

10    continue

      return

      end




      function vnorm(x,n)

      implicit double precision (a-h,o-z)

      dimension x(n)
      external dot

      vnorm = sqrt(dot(x,x,n))

      return

      end




      subroutine vector(c,nrow,j,v,nrdc)

      implicit double precision (a-h,o-z)

      dimension c(nrdc,1),v(nrow)

      do 10 i = 1,nrow

         v(i) = c(i,j)

10    continue

      return

      end




      subroutine ccheck(c,ncol,nq,notind,nrdc,numest)

c   This subroutine checks whether the columns of c are
c   linearly independent

c   The parameter and dimension statements in the main, ccheck, 
c   hyptest, estim, and power subroutines must be in accord.

      parameter (numpar=75,numestd=numpar)

      implicit double precision (a-h,o-z)

      dimension c(nrdc,numestd),uc(numpar,numestd),
     xvc(numestd,numestd),
     xxlamc(numestd),e(numestd)
      dimension ctemp(numpar,numestd)
      dimension work(numpar)

      common /threshmc/ threshm

      notind = 0

      do 5 i = 1,ncol

         do 4 j = 1,nq

            ctemp(i,j) = c(i,j)

 4       continue

 5    continue

c
c   ssvdc is the LINPACK singular value decomposition routine
c

      call dsvdc(ctemp,numpar,ncol,nq,xlamc,e,uc,numpar,
     xvc,numestd,work,21,info)

      xnorm = 0.0

      do 10 i = 1,nq

         xnorm = xnorm + xlamc(i)**2

 10   continue

      xnorm = sqrt(xnorm)
      thresh = xnorm*threshm

      izerof = 0

      do 20 i = 1,nq

         if (xlamc(i) .gt. thresh) go to 20
         ifirst = i
         izerof = 1
         go to 30

 20   continue

 30   if (izerof .eq. 0) return
      irank = ifirst - 1

      write(6,40)
      write(9,40)
 40   format(/,1x,'The columns of c = (c1 ... cq) appear',
     x' to be linearly dependent.',/)

      write(6,50) irank
      write(9,50) irank
 50   format(/,1x,'The rank of c appears to be ',i2,/)

      write(6,60)
      write(9,60)
 60   format(/,1x,'The set of all d such that (c1 ... cq)d = 0',
     x' (approximately)',/,1x,'is the space spanned by the ',
     x'following vectors',//)

      do 80 i = 1,nq

         write(6,70) (vc(i,j),j = ifirst,nq)
         write(9,70) (vc(i,j),j = ifirst,nq)
 70      format(1x,10(e10.2,2x))

 80   continue

      notind = 1

      return

      end




      subroutine hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
     x             notestf,idesign)

c
c   If idesign equals 0,
c   this subroutine checks whether the c's are estimable and
c   if they are it tests the hypothesis c1'b = ... = cq'b = 0.
c   If idesign equals 1 (estimability calculations requested) or
c   idesign equals 2 (power calculations requested),
c   this subroutine checks the estimability of the c's and provides for
c   power calculations.
c

c   The parameter and dimension statements in the main, ccheck, 
c   hyptest, estim, and power subroutines must be in accord.

      parameter (nrowd=400,numpar=75,numest=numpar)

      implicit double precision (a-h,o-z)

      dimension x(nrdx,ncol),y(nrdx),c(nrdc,ncol),e(numpar),work(nrowd)
      dimension f(numpar,numpar),ft(numpar,numpar),fft(numpar,numpar)
      dimension vtemp(nrowd),temp1(nrowd),temp2(nrowd)
      dimension g(nrowd,numpar),h(nrowd,numpar)
      dimension ua(nrowd,numest),xlama(numest),va(numest,numest)
      dimension notest(numest)
      dimension at(numest,nrowd),ata(numest,numest)
      dimension det(2)

      common /hec/ ux(nrowd,numpar),vx(numpar,numpar),xlamx(numpar),
     x             a(nrowd,numest)
      common /hec2/ resss,iresdof
      common /ataic/ atai(numest,numest)
      common /ixrankc/ ixrank
      common /threshmc/ threshm

      external vnorm

c
c   notestf is set to 1 if some of the c's are not estimable
c

      notestf = 0

      do 5 i = 1,nrow

         do 4 j = 1,ncol

            g(i,j) = x(i,j)

 4       continue

 5    continue

c
c   dsvdc is the LINPACK singular value decomposition routine
c

      call dsvdc(g,nrdx,nrow,ncol,xlamx,e,ux,nrdx,vx,numpar,
     xwork,21,info)

c
c   If info is not equal to 0 print an error message and stop.
c

      if (info .ne. 0) then

         write(6,7) info
 7       format(1x,'ERROR hyptest 7:  The info value from a singular',
     x   ' value decomposition was ',i3,'.',/)
         stop

      endif

c
c   Obtain the threshold that is used to determine whether singular values
c   are "essentially zero"
c

      xnorm = 0.0

      do 10 i = 1,ncol

         xnorm = xnorm + xlamx(i)**2

 10   continue

      xnorm = sqrt(xnorm)
      thresh = xnorm*threshm

c
c   Determine the rank of X
c

      do 120 i = 1,ncol

         if (xlamx(i) .gt. thresh) go to 120

         ixrank = i - 1
         go to 130

 120  continue

      ixrank = ncol

c
c   Determine whether the c's are estimable (there exists an a
c   such that X'a = c)
c

 130  iresdof = nrow - ixrank

      do 150 i = 1,ncol

         do 140 j = 1,ixrank

            f(i,j) = vx(i,j)

 140     continue

 150  continue

c
c   Obtain the transpose of f
c

      do 22 i = 1,ncol

         do 12 j = 1,ixrank

            ft(j,i) = f(i,j)

 12      continue

 22   continue

c
c   Multiply f by ft to obtain fft
c

      do 52 j = 1,ncol

         do 42 i = 1,ncol

            fft(i,j) = 0.0

            do 32 k = 1,ixrank

               fft(i,j) = fft(i,j) + f(i,k)*ft(k,j)

 32         continue

 42      continue

 52   continue

c
c   fft is the projection operator onto the space spanned by
c   the columns of V (or the columns of X')
c

      mm = 0

      do 200 i = 1,nq

         call vector(c,ncol,i,vtemp,nrdc)
         thresh = vnorm(vtemp,ncol)*threshm
         call matmult(fft,ncol,ncol,vtemp,ncol,1,temp1,numpar,nrowd)
         call matsub(vtemp,temp1,temp2,ncol,1,nrowd)
         if (vnorm(temp2,ncol) .lt. thresh) go to 200

         mm = mm + 1
         notest(mm) = i

 200  continue

      if (mm .eq. 0) go to 300

      write(6,210) (notest(i),i = 1,mm)
      write(9,210) (notest(i),i = 1,mm)
 210  format(///,1x,"The following c's do not appear to be estimable",
     1//,(1x,10(i3,2x)),/)

      notestf = 1

      return

c
c   Calculate the ai's that satisfy    X'ai = ci
c

 300  if (idesign .ne. 1) go to 305

      write(6,302)
      write(9,302)
 302  format(///,1x,"The c's all appear to be estimable",/)

 305  do 320 i = 1,nrow

         do 310 j = 1,ixrank

            g(i,j) = ux(i,j)/xlamx(j)

 310     continue

 320  continue

c
c   H is the matrix that transforms c's into a's
c   (H = u(nxr)*diaginv(rxr)*vt(rxp))
c

      call matmult(g,nrow,ixrank,ft,ixrank,ncol,h,nrowd,numpar)
      call matmult(h,nrow,ncol,c,ncol,nq,a,nrowd,nrdc)
      if (idesign .eq. 0) go to 323

      call mattran(a,nrow,nq,at,nrowd,numest)
      call matmult(at,nq,nrow,a,nrow,nq,ata,numest,nrowd)

c
c   Taken together, dpofa and dpodi yield an inverse of the
c   ata matrix.  dpofa makes use of the Cholesky algorithm 
c   to factor the positive definite matrix ata.
c   dpofa and dpodi are LINPACK routines.  
c   See the LINPACK manual
c   for details.  (Dongarra, Moler, Bunch, and Stewart,
c   LINPACK Users' Guide, SIAM, 1979)
c

      call dpofa(ata,numest,nq,info)

      if (info .ne. 0) then

         write(6,9999) info
 9999    format(/,1x,'ERROR 9999:  The info value from a matrix',
     x   ' inversion was',1x,i3,'.',/)
         stop

      endif

c
c   Calculate the inverse of ata
c

      job = 1

      call dpodi(ata,numest,nq,det,job)

      do 322 i = 1,nq

         do 321 j = i,nq

            atai(i,j) = ata(i,j)
            atai(j,i) = ata(i,j)

 321     continue

 322  continue

      return


 323  do 325 i = 1,nrow

         do 324 j = 1,nq

            g(i,j) = a(i,j)

 324     continue

 325  continue

c
c   Print the a vectors
c

      write(8,326)
 326  format(//,1x,"The a's are (as columns)",/)

      do 328 i = 1,nrow

         write(8,327) (a(i,j),j = 1,nq)
 327     format(10(1x,e11.3))

 328  continue

c
c   Perform a singular value decomposition of A
c   in order to obtain an orthonormal basis of the
c   linear span of A.
c   dsvdc is a LINPACK singular value decomposition routine.
c

      call dsvdc(g,nrowd,nrow,nq,xlama,e,ua,nrowd,va,numest,
     xwork,20,infoa)

      if (infoa .ne. 0) then

         write(6,329) infoa
 329     format(/,1x,'ERROR 329: The info value from a singular',
     x   ' value decomposition was ',i3,'.',/)
         stop

      endif

c
c   obtain the hypothesis sum of squares
c
      hypss = 0.0

      do 340 j = 1,nq

         temp0 = 0.0

         do 330 i = 1,nrow

            temp0 = temp0 + ua(i,j)*y(i)

 330     continue

         hypss = hypss + temp0**2

 340  continue


      ihypdof = nq

c
c   obtain the residual sum of squares
c

      xss = 0.0

      do 440 j = 1,ixrank

         temp0 = 0.0

         do 430 i = 1,nrow

            temp0 = temp0 + ux(i,j)*y(i)

 430     continue

         xss = xss + temp0**2

 440  continue

      totss = 0.0

      do 450 i = 1,nrow

         totss = totss + y(i)**2

 450  continue

      resss = totss - xss

      return

      end




      subroutine estim(x,y,nrdx,nrow,ncol,c,nrdc,nc,alphai,alphaj,
     x           est,cli,cui,clj,cuj,notestf)

c
c   This subroutine checks whether the c's are estimable and
c   if they are it obtains confidence intervals.
c

c   The parameter and dimension statements in the main, ccheck, 
c   hyptest, estim, and power subroutines must be in accord.

      parameter (nrowd=400,numpar=75,numest=numpar)

      implicit double precision (a-h,o-z)

      dimension x(nrdx,ncol),y(nrdx),c(nrdc,numest),e(numpar),
     xwork(nrowd)
      dimension f(numpar,numpar),ft(numpar,numpar),fft(numpar,numpar)
      dimension vtemp(nrowd),temp1(nrowd),temp2(nrowd)
      dimension g(nrowd,numpar),h(nrowd,numpar)
      dimension est(numest),cli(numest),cui(numest),clj(numest),
     xcuj(numest)
      dimension notest(numest)

c
c   Notice the restrictions on the number of c's
c

      dimension xlamc(numpar),uc(numpar,numpar),vc(numpar,numpar)

      common /hec/ ux(nrowd,numpar),vx(numpar,numpar),xlamx(numpar),
     x             a(nrowd,numest)
      common /hec2/ resss,iresdof
      common /ixrankc/ ixrank
      common /threshmc/ threshm

      external vnorm

      notestf = 0

      do 5 i = 1,nrow

         do 4 j = 1,ncol

            g(i,j) = x(i,j)

 4       continue

 5    continue

c
c   ssvdc is the LINPACK singular value decomposition routine
c

      call dsvdc(g,nrowd,nrow,ncol,xlamx,e,ux,nrowd,vx,numpar,work,
     x21,info)

c
c   If info is not equal to 0 print an error message and stop.
c

      if (info .ne. 0) then

         write(6,7) info
 7       format(1x,'ERROR estim 7:  The info value from a singular',
     x   ' value decomposition was ',i3,'.',/)
         stop

      endif

c
c   Obtain the threshold that is used to determine whether singular values
c   are "essentially zero"
c

      xnorm = 0.0

      do 10 i = 1,ncol

         xnorm = xnorm + xlamx(i)**2

 10   continue

      xnorm = sqrt(xnorm)
      thresh = xnorm*threshm

c
c   Determine the rank of X
c

      do 120 i = 1,ncol

         if (xlamx(i) .gt. thresh) go to 120

         ixrank = i - 1
         go to 1000

 120  continue

      ixrank = ncol

c
c   Determine the number of linearly independent c's
c

 1000 if (ncol .gt. nc) then

         do 1003 i = 1,nrowd

            do 1002 j = 1,numpar

               g(i,j) = 0.0

 1002       continue

 1003    continue

c
c   C is made up of columns of c's
c

         do 1005 i = 1,ncol

            do 1004 j = 1,nc

               g(i,j) = c(i,j)

 1004       continue

 1005    continue

c
c   ssvdc is the LINPACK singular value decomposition routine
c

         call dsvdc(g,nrowd,ncol,nc,xlamc,e,uc,numpar,vc,numpar,work,
     x   21,info)

c
c   Obtain the threshold that is used to determine whether singular values
c   are "essentially zero"
c
         xnorm = 0.0

         do 1010 i = 1,nc

            xnorm = xnorm + xlamc(i)**2

 1010    continue

         xnorm = sqrt(xnorm)
         thresh = xnorm*threshm

c
c   Determine the rank of C
c

         do 1120 i = 1,nc

            if (xlamc(i) .gt. thresh) go to 1120

            icrank = i - 1
            go to 130

 1120    continue

         icrank = nc

      else

c
c   nc >= ncol
c

         do 2003 i = 1,nrowd

            do 2002 j = 1,numpar

               g(i,j) = 0.0

 2002       continue

 2003    continue

         do 2005 i = 1,ncol

            do 2004 j = 1,nc

               g(j,i) = c(i,j)

 2004       continue

 2005    continue

c
c   ssvdc is the LINPACK singular value decomposition routine
c

         call dsvdc(g,nrowd,nc,ncol,xlamc,e,uc,numpar,vc,numpar,work,
     x   21,info)

c
c   Obtain the threshold that is used to determine whether singular values
c   are "essentially zero"
c

         xnorm = 0.0

         do 2010 i = 1,ncol

            xnorm = xnorm + xlamc(i)**2

 2010    continue

         xnorm = sqrt(xnorm)
         thresh = xnorm*threshm

c
c   Determine the rank of C
c

         do 2120 i = 1,ncol

            if (xlamc(i) .gt. thresh) go to 2120

            icrank = i - 1
            go to 130

 2120    continue

         icrank = ncol

      endif

c
c   Determine whether the c's are estimable (there exists an a
c   such that X'a = c)
c

 130  do 150 i = 1,ncol

         do 140 j = 1,ixrank

            f(i,j) = vx(i,j)

 140     continue

 150  continue

c
c   Obtain the transpose of f
c

      do 22 i = 1,ncol

         do 12 j = 1,ixrank

            ft(j,i) = f(i,j)

 12      continue

 22   continue

c
c   Multiply f by ft to obtain fft
c

      do 52 j = 1,ncol

         do 42 i = 1,ncol

            fft(i,j) = 0.0

            do 32 k = 1,ixrank

               fft(i,j) = fft(i,j) + f(i,k)*ft(k,j)

 32         continue

 42      continue

 52   continue

c
c   fft is the projection operator onto the space spanned by
c   the columns of V (or the columns of X')
c

      mm = 0

      do 200 i = 1,nc

         call vector(c,ncol,i,vtemp,nrdc)
         thresh = vnorm(vtemp,ncol)*threshm
         call matmult(fft,ncol,ncol,vtemp,ncol,1,temp1,numpar,nrowd)
         call matsub(vtemp,temp1,temp2,ncol,1,nrowd)
         if (vnorm(temp2,ncol) .lt. thresh) go to 200

         mm = mm + 1
         notest(mm) = i

 200  continue

      if (mm .eq. 0) go to 300

      write(6,210) (notest(i),i = 1,mm)
      write(9,210) (notest(i),i = 1,mm)
 210  format(///,1x,"The following c's do not appear to be estimable",
     x//,(1x,10(i2,2x)),/)

      notestf = 1

      return

c
c   Calculate the ai's that satisfy  X'ai = ci
c

 300  do 320 i = 1,nrow

         do 310 j = 1,ixrank

            g(i,j) = ux(i,j)/xlamx(j)

 310     continue

 320  continue

c
c   H is the matrix that transforms c's into a's
c   (h = u(nxr)*diaginv(rxr)*vt(rxp))
c

      call matmult(g,nrow,ixrank,ft,ixrank,ncol,h,nrowd,numpar)
      call matmult(h,nrow,ncol,c,ncol,nc,a,nrowd,nrdc)

c
c   Print the a vectors
c

      write(8,326)
 326  format(//,1x,"The a's are (as columns)",/)

      do 328 i = 1,nrow

         write(8,327) (a(i,j),j = 1,nc)
 327     format(10(1x,e11.3))

 328  continue

c
c   Obtain the residual vector
c

      xss = 0.0

      do 440 j = 1,ixrank

         temp0 = 0.0

         do 430 i = 1,nrow

            temp0 = temp0 + ux(i,j)*y(i)

 430     continue

         xss = xss + temp0**2

 440  continue

      totss = 0.0

      do 450 i = 1,nrow

         totss = totss + y(i)**2

 450  continue

      resss = totss - xss

      iresdof = nrow - ixrank

c
c   Calculate the individual and the joint confidence intervals
c

 500  s = sqrt(resss/iresdof)

      xxx = icrank
      yyy = iresdof

c
c   cdft is a routine from DCDFLIB.  It performs calculations
c   associated with the t distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

      p = 1.0 - alphai/2.0
      q = 1.0d0 - p

      call cdft(2,p,q,t,yyy,istatus,bound)

      if (istatus .ne. 0) then

         write(6,502) istatus
 502     format(/,1x,'ERROR 502:  The istatus value from a T inverse',
     x   ' routine was ',i2,'.',/)

         stop

      endif

c
c   cdff is a routine from DCDFLIB.  It performs calculations
c   associated with the F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

      p = 1.0 - alphaj
      q = alphaj

      call cdff(2,p,q,fvalue,xxx,yyy,istatus,bound)

      if (istatus .ne. 0) then

         write(6,504) istatus
 504     format(/,1x,'Error 504: The istatus value from an',
     x   ' F inverse',
     x   ' routine was ',i2,'.',/)

         stop

      endif

      delta1 = s*t
      delta2 = s*sqrt(icrank*fvalue)

      do 600 i = 1,nc

         aa = 0.0
         xxx = 0.0

         do 550 j = 1,nrow

            aa = aa + a(j,i)*a(j,i)
            xxx = xxx + a(j,i)*y(j)

 550     continue

         aa = sqrt(aa)
         est(i) = xxx

         delta = aa*delta1
         cli(i) = xxx - delta
         cui(i) = xxx + delta

         delta = aa*delta2
         clj(i) = xxx - delta
         cuj(i) = xxx + delta

 600  continue

      return

      end




      subroutine power(nq,cb1,numcb,sigma,numsig,cb,powerl,iresdof)

c
c   This subroutine performs power calculations.
c

c   The parameter and dimension statements in the main, ccheck, 
c   hyptest, estim, and power subroutines must be in accord.


      parameter (numest=75)

      implicit double precision (a-h,o-z)

      dimension cb1(5),sigma(5),cb(numest)
      dimension powera(5,5,3)
      dimension numrep(5,5,3)

      dimension halfwidt(3,5)

      common /alphac/ alpha(3)
      common /ataic/ atai(numest,numest)
      common /nrowc/ nrow
      common /ixrankc/ ixrank

      if (nq .gt. 1) go to 1000

      d = sqrt(atai(1,1))

      do 200 ii = 1,3

         f = iresdof

c
c   cdft is a routine from DCDFLIB.  It performs calculations
c   associated with the t distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

         p = 1.0 - alpha(ii)/2.0
         call cdft(2,p,1.0d0-p,xbase,f,istatus,bound)

         if (istatus .ne. 0) then

            write(6,10) istatus
 10         format(/,1x,'ERROR 10:  The istatus value from a T inverse',
     x      ' routine was ',i2,'.',/)

            stop

         endif

         do 100 i = 1,numsig

            do 90 j = 1,numcb

               xncp = (cb1(j)*d)/sigma(i)

c
c   cdffnc is a routine from DCDFLIB.  It performs calculations
c   associated with the non-central F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

               f = iresdof

               call cdffnc(1,p,q,xbase*xbase,1.0d0,f,xncp*xncp,
     x                     istatus,bound)

               if (istatus .ne. 0) then

                  write(6,1025) istatus

                  stop

               endif

               powera(i,j,ii) = 1.0d0 - p

c
c   Find the number of replications needed to achieve the
c   desired power.  ("1 replication" means the original design,
c   "k replications" means k copies of the original design.)
c

               numrep(i,j,ii) = 1

               if (powerl .le. powera(i,j,ii)) go to 90

               numl = 1
               numh = 50

               xncpt = xncp*sqrt(50.0)
               idof = 50*nrow - ixrank
               f = idof

c
c   cdft is a routine from DCDFLIB.  It performs calculations
c   associated with the t distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

               p = 1.0 - alpha(ii)/2.0
               call cdft(2,p,1.0d0-p,x,f,istatus,bound)

               if (istatus .ne. 0) then

                  write(6,10) istatus

                  stop

               endif

c
c   cdffnc is a routine from DCDFLIB.  It performs calculations
c   associated with the non-central F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

               call cdffnc(1,p,q,x*x,1.0d0,f,xncpt*xncpt,
     x                     istatus,bound)

               if (istatus .ne. 0) then

                  write(6,1025) istatus

                  stop

               endif

               pt = 1.0d0 - p

               if (pt .ge. powerl) go to 50

c
c   numrep is set to 100 so that ** is printed in the power table.
c   ** indicates that the necessary number of replications is
c   greater than 50.
c

               numrep(i,j,ii) = 100
               go to 90

 50            idel = numh - numl

               if (idel .lt. 2) then

                  numrep(i,j,ii) = numh
                  go to 90

               endif
                  
               numt = numl + idel/2
               xnumt = numt
               xncpt = xncp*sqrt(xnumt)
               idof = numt*nrow - ixrank
               f = idof

c
c   cdft is a routine from DCDFLIB.  It performs calculations
c   associated with the t distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

               p = 1.0 - alpha(ii)/2.0
               call cdft(2,p,1.0d0-p,x,f,istatus,bound)

               if (istatus .ne. 0) then

                  write(6,10) istatus

                  stop

               endif

c
c   cdffnc is a routine from DCDFLIB.  It performs calculations
c   associated with the non-central F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

               call cdffnc(1,p,q,x*x,1.0d0,f,xncpt*xncpt,
     x                     istatus,bound)

               if (istatus .ne. 0) then

                  write(6,1025) istatus

                  stop

               endif


               pt = 1.0d0 - p


               if (pt .ge. powerl) go to 70

               numl = numt
               go to 50

 70            numh = numt
               go to 50

 90         continue

c
c   Calculate approximate half-widths of confidence intervals
c

            halfwidt(ii,i) = xbase*sigma(i)/d

 100     continue

 200  continue

      write(6,205)
      write(9,205)
 205  format(///,1x,"Power tables",//,1x,"(Here q = 1, and we consider",
     x" a range of c'b values and a range of sigma values.)",/)

      do 300 ii = 1,3

         write(6,310) alpha(ii)
         write(9,310) alpha(ii)
 310     format(///,1x,'The significance level is  ',f4.2,//)

         write(6,320)
         write(9,320)
 320     format(/,1x,"The c'b values are on the top.  The sigma values",
     x   " are along the left side.",//)

         write(6,330) (cb1(i),i = 1,numcb)
         write(9,330) (cb1(i),i = 1,numcb)
 330     format(/,18x,5(e12.4,3x))

         write(6,332)
         write(9,332)
 332     format(/)

         do 340 i = 1,numsig

            write(6,335) sigma(i),(powera(i,j,ii),j = 1,numcb)
            write(9,335) sigma(i),(powera(i,j,ii),j = 1,numcb)
 335        format(1x,e12.4,5x,5(e12.4,3x))

 340     continue

 300  continue


      write(6,355)
      write(9,355)
 355  format(///,1x,"Halfwidth tables",//,1x,"(Here q = 1,",
     x" and we consider",
     x" a range of sigma values.)",/)

      do 390 ii = 1,3

         write(6,360) alpha(ii)
         write(9,360) alpha(ii)
 360     format(///,1x,'The significance level is  ',f4.2,//)

         write(6,365)
         write(9,365)
 365     format(/,1x,"The sigma values",
     x   " are along the left side.",//)

         do 375 i = 1,numsig

            write(6,370) sigma(i),halfwidt(ii,i)
            write(9,370) sigma(i),halfwidt(ii,i)
 370        format(1x,e12.4,5x,5(e12.4,3x))

 375     continue

 390  continue


      write(6,405)
      write(9,405)
 405  format(///,1x,'Replication tables',//,
     x1x,'These tables give us the number of replications of the ',
     x'current design',/,1x,'that are needed to achieve the desired',
     x' power.',/)

      write(6,410) powerl
      write(9,410) powerl
 410  format(/,1x,'The desired power is  ',f4.2,//)

      do 500 ii = 1,3

         write(6,420) alpha(ii)
         write(9,420) alpha(ii)
 420     format(///,1x,'The significance level is  ',f4.2,//)

         write(6,430)
         write(9,430)
 430     format(/,1x,"The c'b values are on the top.  the sigma values",
     x   " are along the left side.",//)

         write(6,440) (cb1(i),i = 1,numcb)
         write(9,440) (cb1(i),i = 1,numcb)
 440     format(/,18x,5(e12.4,3x))

         write(6,442)
         write(9,442)
 442     format(/)

         do 450 i = 1,numsig

            write(6,445) sigma(i),(numrep(i,j,ii),j = 1,numcb)
            write(9,445) sigma(i),(numrep(i,j,ii),j = 1,numcb)
 445        format(1x,e12.4,5x,5(5x,i2,8x))

 450     continue

 500  continue

      return



 1000 do 1200 ii = 1,3

         p = 1.0 - alpha(ii)
         d1 = nq
         d2 = iresdof

c
c   cdff is a routine from DCDFLIB.  It performs calculations
c   associated with the F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

         call cdff(2,p,1.0d0-p,xbase,d1,d2,istatus,bound)

         if (istatus .ne. 0) then

            write(6,1005) istatus
 1005       format(/,1x,'Error 1005: The istatus value from an',
     x      ' F inverse',
     x      ' routine was ',i2,'.',/)

            stop

         endif

         do 1100 i = 1,numsig

            xncp = 0.0

            do 1020 j = 1,nq

               do 1010 k = 1,nq

                  xncp = xncp + atai(j,k)*cb(j)*cb(k)

 1010          continue

 1020       continue

            xncp = xncp/(sigma(i)**2)

c
c   call the non-central f routine and obtain powera(i,1,ii)
c

c
c   cdffnc is a routine from DCDFLIB.  It performs calculations
c   associated with the non-central F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

            d2 = iresdof

            call cdffnc(1,p,q,xbase,d1,d2,xncp,istatus,bound)

            powera(i,1,ii) = 1.0d0 - p
   
            if (istatus .ne. 0) then

               write(6,1025) istatus
 1025          format(/,1x,'Error 1025: The istatus value from',
     x         ' a non-central F cdf',
     x         ' routine was ',i2,'.',/)

               stop

            endif

c
c   find the number of replications needed to achieve the
c   desired power ("1 replication" means the original design,
c   "k replications" means k copies of the original design)
c
            numrep(i,1,ii) = 1

            if (powerl .le. powera(i,1,ii)) go to 1100

            numl = 1
            numh = 50

            xncpt = xncp*50.0d0
            idof = 50*nrow - ixrank
            d2 = idof

c
c   cdff is a routine from DCDFLIB.  It performs calculations
c   associated with the F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

            p = 1.0 - alpha(ii)

            call cdff(2,p,1.0d0-p,x,d1,d2,istatus,bound)

            if (istatus .ne. 0) then

               write(6,1035) istatus
 1035          format(/,1x,'Error 1035: The istatus value from an',
     x         ' F inverse',
     x         ' routine was ',i2,'.',/)

               stop

            endif

c
c   cdffnc is a routine from DCDFLIB.  It performs calculations
c   associated with the non-central F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

            call cdffnc(1,p,q,x,d1,d2,xncpt,istatus,bound)

            pt = 1.0d0 - p
   
            if (istatus .ne. 0) then

               write(6,1040) istatus
 1040          format(/,1x,'Error 1040: The istatus value from',
     x         ' a non-central F cdf',
     x         ' routine was ',i2,'.',/)

               stop

            endif

c
            if (pt .ge. powerl) go to 1050
c
c   numrep is set to 100 so that ** is printed in the power table.
c   ** indicates that the necessary number of replications is
c   greater than 50.
c
            numrep(i,1,ii) = 100
            go to 1100

 1050       idel = numh - numl
            if (idel .lt. 2) then
               numrep(i,1,ii) = numh
               go to 1100
            endif

 1060       numt = numl + idel/2
            xnumt = numt
            xncpt = xncp*xnumt
            idof = numt*nrow - ixrank
            d2 = idof

c
c   cdff is a routine from DCDFLIB.  It performs calculations
c   associated with the F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

            p = 1.0 - alpha(ii)

            call cdff(2,p,1.0d0-p,x,d1,d2,istatus,bound)

            if (istatus .ne. 0) then

               write(6,1062) istatus
 1062          format(/,1x,'Error 1062: The istatus value from an',
     x         ' F inverse',
     x         ' routine was ',i2,'.',/)

               stop

            endif

c
c   cdffnc is a routine from DCDFLIB.  It performs calculations
c   associated with the non-central F distribution.   DCDFLIB is a
c   public domain library of "routines for cumulative distribution
c   functions, their inverses, and their parameters."  It was produced
c   by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and 
c   Kathy Russell of the Department of Biomathematics, 
c   M.D. Anderson Cancer Center, The University of Texas.  
c   The library can be found at
c
c   http://odin.mdacc.tmc.edu/anonftp/source.html
c

            call cdffnc(1,p,q,x,d1,d2,xncpt,istatus,bound)

            pt = 1.0d0 - p
   
            if (istatus .ne. 0) then

               write(6,1064) istatus
 1064          format(/,1x,'Error 1064: The istatus value from',
     x         ' a non-central F cdf',
     x         ' routine was ',i2,'.',/)

               stop

            endif


            if (pt .ge. powerl) go to 1070
            numl = numt
            go to 1050

 1070       numh = numt
            go to 1050

 1100    continue

 1200 continue

      write(6,1205)
      write(9,1205)
 1205 format(///,1x,'Power tables',//)

      write(6,1210) (cb(i),i = 1,nq)
      write(9,1210) (cb(i),i = 1,nq)

 1210 format(/,1x,"The c1'b, ... , cq'b values are",//,1x,
     x5(e12.4,2x),/,1x,5(e12.4,2x),/)

      do 1300 ii = 1,3

         write(6,1310) alpha(ii)
         write(9,1310) alpha(ii)
 1310    format(///,1x,'The significance level is  ',f4.2,//)

         write(6,1320)
         write(9,1320)
 1320    format(/,1x,'The sigma and associated power values are',//)

         do 1340 i = 1,numsig

            write(6,1335) sigma(i),(powera(i,1,ii))
            write(9,1335) sigma(i),(powera(i,1,ii))
 1335       format(1x,e12.4,5x,e12.4)

 1340    continue

 1300 continue

      write(6,1405)
      write(9,1405)
 1405 format(///,1x,'Replication tables',//,
     x1x,'These tables give us the number of replications of the ',
     x'current design',/,1x,'that are needed to achieve the desired',
     x' power.',/)

      write(6,1410) powerl
      write(9,1410) powerl
 1410 format(/,1x,'The desired power is  ',f4.2,//)

      do 1500 ii = 1,3

         write(6,1420) alpha(ii)
         write(9,1420) alpha(ii)
 1420    format(///,1x,'The significance level is  ',f4.2,//)

         write(6,1430)
         write(9,1430)
 1430    format(/,1x,'The sigma and associated',
     x   ' replication values are',//)

         do 1450 i = 1,numsig

            write(6,1445) sigma(i),numrep(i,1,ii)
            write(9,1445) sigma(i),numrep(i,1,ii)
 1445       format(1x,e12.4,10x,i2)

 1450    continue

 1500 continue

      return

      end