SUBROUTINE POLY(Y,X,W,N,IDEG,IWRITE,B,SDB,S,DF,PRED,RES)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES A LEAST SQUARES
C              POLYNOMIAL FIT (OF DEGREE = IDEG) OF THE
C              RESPONSE VARIABLE DATA IN THE SINGLE PRECISION
C              VECTOR Y AS A FUNCTION OF THE INDEPENDENT
C              VARIABLE DATA IN THE SINGLE PRECISION
C              VECTOR X.
C     INPUT  ARGUMENTS--Y      = SINGLE PRECISION VECTOR OF 
C                                RESPONSE DATA (THAT IS, THE
C                                DEPENDENT VARIABLE).
C                     --X      = SINGLE PRECISION VECTOR OF 
C                                THE INDEPENDENT VARIABLE.
C                     --W      = THE SINGLE PRECISION VECTOR
C                                OF WEIGHTS FOR THE RESPONSE
C                                VARIABLE.
C                     --N      = THE INTEGER VALUE OF THE SAMPLE SIZE.
C                     --IDEG   = THE INTEGER VALUE OF THE DESIRED
C                                DEGREE OF THE POLYNOMIAL
C                                TO BE FIT.
C                     --IWRITE = THE INTEGER VALUE WHICH IF ZERO WILL 
C                                RESULT IN NO PRINTED OUTPUT, AND IF
C                                NON-ZERO (E.G., 1) WILL RESULT IN
C                                SOME LIMITED PRINTED OUTPUT
C                                (COEFFICIENTS, STANDARD DEVIATIONS OF
C                                COEFFICIENTS, RESIDUAL STANDARD DEVIATION).
C     OUTPUT ARGUMENTS--B      = THE SINGLE PRECISION VECTOR OF
C                                ESTIMATED REGRESSION COEFFICIENTS.
C                     --SDB    = THE SINGLE PRECISION VECTOR OF
C                                ESTIMATED STANDARD DEVIATIONS OF THE 
C                                ESTIMATED REGRESSION COEFFICIENTS.
C                     --S      = THE ESTIMATED RESIDUAL STANDARD
C                                DEVIATION.
C                     --DF     = THE DEGREES OF FREEDOM
C                                ASSOCIATED WITH THE RESIDUAL
C                                STANDARD DEVIATION =
C                                NUMBER OF OBSERVATIONS MINUS
C                                NUMBER OF PARAMETERS =
C                                N - (IDEG + 1).
C                     --PRED   = THE SINGLE PRECISION VECTOR OF
C                                PREDICTED VALUES FROM THE
C                                LEAST SQUARES FIT.
C                     --RES    = THE SINGLE PRECISION VECTOR OF
C                                RESIDUALS FROM THE LEAST SQUARES FIT.
C     SUBROUTINES NEEDED--DECOMP, INVXWX, DOT, FCDF.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING LABORATORY
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE:  301-921-2315
C     ORIGINAL VERSION--MARCH     1974
C     UPDATED--OCTOBER   1974 
C     UPDATED--MARCH     1975 
C     UPDATED--MAY       1975 
C     UPDATED--JULY      1975 
C     UPDATED--SEPTEMBER 1975 
C     UPDATED--NOVEMBER  1975.
C     UPDATED--FEBRUARY  1976.
C     UPDATED--JUNE      1976.
C     UPDATED--OCTOBER   1976.
C     UPDATED--MAY       1977.
C     UPDATED--JUNE      1977.
C
C---------------------------------------------------------------------
C
      DIMENSION Y(1),X(1),W(1),B(1),SDB(1),PRED(1),RES(1)
      DIMENSION B2(50)
      DIMENSION F(3000),WRES(3000),G(50),H(50)
      DIMENSION Q(10000),R(2500),D(50),IPIVOT(50) 
      COMMON /BLOCK2/ WS(15000)
      COMMON /BLOCK3/ DUM1(3000),DUM2(3000)
      EQUIVALENCE (Q(1),WS(1))
      EQUIVALENCE (R(1),WS(10001))
      EQUIVALENCE (D(1),WS(12501))
      EQUIVALENCE (IPIVOT(1),WS(12551)) 
C
      IPR=6
      AN=N
      K=IDEG+1
      AK=K
      NK=N*K
      KK=K*K
      NMAX=3000
      KMAX=50
      NKMAX=10000
C
C-----START POINT-----------------------------------------------------
C
C     WRITE OUT THE TITLE
C
      IF(IWRITE.EQ.0)GOTO150
      WRITE(IPR,999)
      WRITE(IPR,999)
      WRITE(IPR,101)
      WRITE(IPR,102)N
      WRITE(IPR,103)IDEG
C
C     PRE-SET THE OUTPUT VARIABLES AND VECTORS
C     TO A LARGE VALUE
C     IN CASE A PREMATURE EXIT OCCURS DUE TO NUMERICAL INSTABILITY.
C
  150 VALUE=(10.0**10)+1000.0 
      S=VALUE
      IF(K.LE.0)GOTO155
      DO160I=1,K
      B(I)=VALUE
      SDB(I)=VALUE
  160 CONTINUE
  155 IF(N.LE.0)GOTO175
      DO170I=1,N
      RES(I)=VALUE
  170 CONTINUE
C
C     CHECK THE INPUT ARGUMENTS N AND K 
C
  175 IF(N.LE.0.OR.N.GT.NMAX)WRITE(IPR,105)N,NMAX 
      IF(N.LE.0.OR.N.GT.NMAX)RETURN
      IF(K.LE.0.OR.K.GT.KMAX)WRITE(IPR,110)K,KMAX 
      IF(K.LE.0.OR.K.GT.KMAX)RETURN
      IF(K.GT.N)WRITE(IPR,115)K,N
      IF(K.GT.N)RETURN
C
C     INSPECT THE WEIGHT VECTOR W--IF ALL ELEMENTS ARE IDENTICAL,
C     THEN RESET ALL ELEMENTS TO 1.0.  THIS AVOIDS THE
C     PROBLEM OF AN UNDEFINED EMPTY WEIGHT VECTOR W WHEN
C     IN FACT AN EQUAL WEIGHTING SCHEME IS DESIRED.
C
      IWFLAG=0
      WHOLD=W(1)
      DO600I=1,N
      IF(W(I).EQ.WHOLD)GOTO600
      GOTO850
  600 CONTINUE
      IWFLAG=1
  850 IF(IWFLAG.EQ.0.AND.IWRITE.NE.0)WRITE(IPR,851)
      IF(IWFLAG.EQ.1.AND.IWRITE.NE.0)WRITE(IPR,852)
C
C     COMPUTE THE ORIGINAL FORM FOR THE Q MATRIX
C     WHICH WILL BE IDENTICAL TO THE DATA MATRIX X
C     OF INDEPENDENT VARIABLES IF THE WEIGHTS
C     SPECIFIED ARE ALL EQUAL.  NOTE THAT THE
C     DATA MATRIX X IS NEVER COMPUTED AS SUCH.
C     NOTE THAT THE DATA MATRIX X IS NOT TO BE
C     CONFUSED WITH THE SINGLE INDEPENDENT
C     VARIABLE VECTOR X.
C     THE Q MATRIX WILL BE CHANGED IN THE DECOMP SUBROUTINE.
C
      DO100J=1,K
      IF(J.EQ.1)GOTO250
      IF(J.EQ.2)GOTO350
      JM1=J-1
      DO200I=1,N
      IQARG1=(I-1)*K+J
      IQARG2=(I-1)*K+JM1
      Q(IQARG1)=Q(IQARG2)*X(I)
  200 CONTINUE
      GOTO100
  250 DO300I=1,N
      IQARG=(I-1)*K+1
      Q(IQARG)=1.0
  300 CONTINUE
      GOTO100
  350 DO400I=1,N
      IQARG=(I-1)*K+2
      Q(IQARG)=X(I) 
  400 CONTINUE
  100 CONTINUE
      IF(IWFLAG.EQ.1)GOTO1440 
      DO1300I=1,N
      DO1400J=1,K
      IQARG=(I-1)*K+J
      Q(IQARG)=Q(IQARG)*SQRT(W(I))
 1400 CONTINUE
 1300 CONTINUE
C
C     COMPUTE ETA AND TOL (FOR THE UNIVAC 1108, ETA = 2**-27)
C     WHICH WILL BE USED IN THE DECOMP SUBROUTINE 
C
 1440 ETA=1.0
 1450 ETA=0.5*ETA
      ETAP1=ETA+1.0 
      IF(ETAP1.GT.1.0)GOTO1450
      TOL=ETA*AK
      NM5=N-5
      NKM5=NK-5
CCCCC WRITE(IPR,1505)(Y(I),I=1,6)
CCCCC WRITE(IPR,1505)(Y(I),I=NM5,N)
CCCCC WRITE(IPR,1505)(X(I),I=1,6)
CCCCC WRITE(IPR,1505)(X(I),I=NM5,N)
CCCCC WRITE(IPR,1505)(W(I),I=1,6)
CCCCC WRITE(IPR,1505)(W(I),I=NM5,N)
CCCCC WRITE(IPR,1505)(Q(IQARG),IQARG=1,6)
CCCCC WRITE(IPR,1505)(Q(IQARG),IQARG=NKM5,NK)
CCCCC WRITE(IPR,1505)(Q(IQARG),IQARG=1,NK)
C
      CALL DECOMP(N,K,ETA,TOL,IRANK,INSING)
C
CCCCC WRITE(IPR,1505)(Q(IQARG),IQARG=1,NK)
CCCCC WRITE(IPR,1505)(R(IRARG),IRARG=1,KK)
CCCCC WRITE(IPR,1505)(D(J),J=1,K)
      IF(INSING.EQ.1)GOTO1550 
      WRITE(IPR,1510)IRANK,K
      RETURN
 1550 KP1=K+1
C
C     *************************************************************** 
C
C     THE PURPOSE OF THIS NEXT SEGMENT (BETWEEN THE STARRED LINES)
C     IS TO SOLVE FOR THE DESIRED REGRESSION COEFFICIENTS.
C     ITERATIVE REFINEMENT IS USED.
C     A SECOND OUTPUT FROM THIS SEGMENT IS THE RESIDUALS FROM THE
C     FINAL FIT.
C     A THIRD OUTPUT FROM THIS SEGMENT IS AN INDICATION (IN THE
C     VARIABLE ICONV) AS TO WHETHER THE ITERATIVE REFINEMENT
C     CONVERGED OR NOT.
C     X--USED IN THIS SEGMENT 
C     Q--USED IN THIS SEGMENT 
C     R--USED IN THIS SEGMENT 
C     D--USED IN THIS SEGMENT 
C     IPIVOT--USED IN THIS SEGMENT
C
C
      ETA2=ETA*ETA
      B2(KP1)=-1.0
      DO 50010 I=1,N
      F(I)=Y(I)
      WRES(I)=0.0
      RES(I)=0.0
50010 CONTINUE
      DO 50020 J=1,K
      B2(J)=0.0
      G(J)=0.0
      H(J)=0.0
50020 CONTINUE
      M=0 
      AMDB2=0.0
      AMDR2=0.0
C
C     BEGIN THE M-TH ITERATION STEP IN THE ITERATIVE REFINEMENT
C
50030 IF (M.LT.2) GO TO 50040 
      IF (((64.*AMDB2.LT.AMDB1).AND.(AMDB2.GT.ETA2*AMB)).OR.((64.*AMDR2.
     1LT.AMDR1).AND.(AMDR2.GT.ETA2*AMR))) GO TO 50040
      GO TO 50250
50040 AMDB1=AMDB2
      AMDR1=AMDR2
      AMDB2=0.0
      AMDR2=0.0
      IF (M.EQ.0) GO TO 50100 
C
C     BEGIN FORMING NEW RESIDUALS
C
      IF(IWFLAG.EQ.0)GOTO50044
      DO50045I=1,N
      WRES(I)=WRES(I)+F(I)
      RES(I)=RES(I)+F(I)
50045 CONTINUE
      GOTO50055
50044 DO 50050 I=1,N
      WRES(I)=WRES(I)+F(I)*SQRT(W(I))
      RES(I)=RES(I)+F(I)/SQRT(W(I))
50050 CONTINUE
50055 DO 50070 IS=1,IRANK
      J=IPIVOT(IS)
      B2(J)=B2(J)+G(IS)
      DO 50060 L=1,N
      IF(J.GE.2)GOTO50065
      DUM1(L)=1.0
      GOTO50060
50065 DUM1(L)=X(L)**(J-1)
50060 CONTINUE
C
      CALL DOT(DUM1,WRES,1,N,0.0,G(IS)) 
50070 G(IS)=-G(IS)
C
      DO 50090 I=1,N
      E=RES(I)
      DO 50080 L=1,K
      IF(L.GE.2)GOTO50085
      DUM1(L)=1.0
      GOTO50080
50085 DUM1(L)=X(I)**(L-1)
50080 CONTINUE
      DUM1(KP1)=Y(I)
C
      CALL DOT(DUM1,B2,1,KP1,E,F(I))
      IF(IWFLAG.EQ.0)F(I)=-F(I)*SQRT(W(I))
      IF(IWFLAG.EQ.1)F(I)=-F(I)
50090 CONTINUE
C
C     END FORMING NEW RESIDUALS
C
50100 DO 50150 IS=1,IRANK
      J=IPIVOT(IS)
      IF (IS.NE.1) GO TO 50110
50110 ISM1=IS-1
      DO 50120 L=1,ISM1
      IRARG=(L-1)*K+IS
50120 DUM1(L)=R(IRARG)
      ANEGGI=-G(IS) 
C
      CALL DOT(DUM1,H,1,ISM1,ANEGGI,H(IS))
      H(IS)=-H(IS)
C
      E=-H(IS)
      DO 50130 L=1,N
      IQARG=(L-1)*K+J
50130 DUM1(L)=Q(IQARG)
C
      CALL DOT(DUM1,F,1,N,E,E)
      E=E/D(IS)
C
      G(IS)=E
      DO 50140 I=1,N
      IQARG=(I-1)*K+J
50140 F(I)=F(I)-E*Q(IQARG)
50150 CONTINUE
C
      DO 50210 IS=1,IRANK
      JS=IRANK+1-IS 
      JSP1=JS+1
      DO 50200 L=JSP1,IRANK
      IRARG=(JS-1)*K+L
50200 DUM1(L)=R(IRARG)
      ANEGGJ=-G(JS) 
C
      CALL DOT(DUM1,G,JSP1,IRANK,ANEGGJ,G(JS))
50210 G(JS)=-G(JS)
C
      DO 50220 IS=1,IRANK
50220 AMDB2=AMDB2+G(IS)*G(IS) 
      DO 50230 I=1,N
50230 AMDR2=AMDR2+F(I)*F(I)
      IF (M.NE.0) GO TO 50240 
      AMB=AMDB2
      AMR=AMDR2
50240 CONTINUE
CCCCC WRITE(IPR,50505)M
CCCCC WRITE(IPR,50506)(B2(I),I=1,K)
CCCCC DO5555I=1,N
CCCCC WRITE(IPR,50506)Y(I),X(I),RES(I),WRES(I),F(I),G(I),H(I)
C5555 CONTINUE
50505 FORMAT(I8)
50506 FORMAT(1H ,8F10.5)
C
C     END THE M-TH ITERATION STEP IN THE ITERATIVE REFINEMENT
C
      M=M+1
      GO TO 50030
50250 IF ((AMDR2.GT.4.*ETA2*AMR).AND.(AMDB2.GT.4.*ETA2*AMB)) GO TO 50260
      ICONV=1
      GO TO 1700
50260 ICONV=0
C
C     *************************************************************** 
C
 1700 IF(ICONV.EQ.1)GOTO1750
      WRITE(IPR,1520)
      RETURN
C
C     ADJUST THE R MATRIX
C     WHICH WILL BE USED IN THE INVERT (INVRR) SUBROUTINE.
C
1750  DO1800J=1,K
      IRARG=(J-1)*K+J
      R(IRARG)=SQRT(D(J))
      IF(J.EQ.K)GOTO1800
      JP1=J+1
      DO1900L=JP1,K 
      IRARG1=(J-1)*K+L
      IRARG2=(J-1)*K+J
      R(IRARG1)=R(IRARG1)*R(IRARG2)
 1900 CONTINUE
 1800 CONTINUE
CCCCC WRITE(IPR,1505)(R(IRARG),IRARG=1,KK)
C
      CALL INVXWX(N,K)
C
CCCCC WRITE(IPR,1505)(R(IRARG),IRARG=1,KK)
C
C     COMPUTE STATISTICAL CALCULATIONS AND THEN WRITE OUT COEFFICIENTS
C     AND STANDARD DEVIATIONS OF COEFFICIENTS ALONG WITH THE
C     RESIDUAL STANDARD DEVIATION.
C
      DO3040I=1,K
      B(I)=B2(I)
 3040 CONTINUE
      DO3070I=1,N
      IF(RES(I).EQ.0.0)GOTO3070
      GOTO3050
 3070 CONTINUE
      WRITE(IPR,3080)
 3080 FORMAT(' ',10X,'NOTE THAT AN EXACT FIT HAS BEEN OBTAINED')
 3050 SUM=0.0
      IF(IWFLAG.EQ.0)GOTO3060 
      DO3055I=1,N
      SUM=SUM+RES(I)*RES(I)
 3055 CONTINUE
      RESSS=SUM
      GOTO3105
 3060 DO3100I=1,N
      SUM=SUM+RES(I)*RES(I)*W(I)
 3100 CONTINUE
 3105 IF(K.EQ.N.AND.IWRITE.NE.0)WRITE(IPR,3110)K
      IF(K.EQ.N)GOTO3250
 3110 FORMAT(1H ,10X,72HNOTE THAT THE NUMBER OF COEFFICIENTS K
     1  = THE SAMPLE SIZE N = ,I8)
      RESDF=AN-AK
      DF=RESDF
      IRESDF=AN-AK+0.5
      RESMS=RESSS/RESDF
      S=SQRT(RESMS) 
 3250 CONTINUE
CCCCC WRITE(IPR,3205)
      DO3300I=1,N
      PRED(I)=Y(I)-RES(I)
CCCCC WRITE(IPR,3305)Y(I),PRED(I),RES(I)
 3300 CONTINUE
      IF(K.EQ.N)GOTO3450
      GOTO3480
 3450 IF(IWRITE.NE.0)WRITE(IPR,3455)
 3455 FORMAT(1H ,21H        J        B(J))
      DO3460I=1,K
      IF(IWRITE.NE.0)WRITE(IPR,3405)J,B(J)
 3460 CONTINUE
      RETURN
 3480 IF(IWRITE.NE.0)WRITE(IPR,3310)S
      IF(IWRITE.NE.0)WRITE(IPR,3311)IRESDF
      IF(IWRITE.NE.0)WRITE(IPR,3312)
      IF(IWRITE.NE.0)WRITE(IPR,3315)
      DO3400J=1,K
      IRARG=(J-1)*K+J
      SDB(J)=S*SQRT(R(IRARG)) 
      T=B(J)/SDB(J) 
      IF(IWRITE.NE.0)WRITE(IPR,3405)J,B(J),SDB(J),T
 3400 CONTINUE
C
C     COMPUTE THE COVARIANCE AND CORRELATION MATRIX OF THE COEFFICIENTS
C
      DO3500I=1,K
      DO3600J=1,K
      IRARG=(I-1)*K+J
      R(IRARG)=R(IRARG)*S*S
 3600 CONTINUE
 3500 CONTINUE
      DO2200I=1,K
      DO2300J=1,K
      IF(I.EQ.J)GOTO2300
      IRARG1=(I-1)*K+J
      IRARG2=(I-1)*K+I
      IRARG3=(J-1)*K+J
      R(IRARG1)=R(IRARG1)/SQRT(R(IRARG2)*R(IRARG3))
 2300 CONTINUE
 2200 CONTINUE
      DO2400J=1,K
      IRARG=(J-1)*K+J
      R(IRARG)=1.0
 2400 CONTINUE
CCCCC WRITE(IPR,2405)
      DO2500I=1,K
      IRAMIN=(I-1)*K+1
      IRAMAX=I*K
CCCCC WRITE(IPR,2505)(R(IRARG),IRARG=IRAMIN,IRAMAX)
 2500 CONTINUE
C
C     CHECK FOR REPLICATION.
C     IF SO, COMPUTE A POOLED STANDARD DEVIATION. 
C     THEN COMPUTE A LACK OF FIT F TEST.
C
C     DETERMINE THE NUMBER OF DISTINCT SUBSETS
C
      NUMSET=0
      DO4200I=1,N
      IF(NUMSET.EQ.0)GOTO4350 
      DO4300J=1,NUMSET
      IF(X(I).EQ.DUM1(J))GOTO4200
 4300 CONTINUE
 4350 NUMSET=NUMSET+1
      DUM1(NUMSET)=X(I)
 4200 CONTINUE
      IF(NUMSET.EQ.0)WRITE(IPR,4205)
      IF(NUMSET.EQ.0)RETURN
C
C     COPY OUT EACH SUBSET INTO THE DUM2 VECTOR
C     AND ANALYZE IT THEREIN
C
      IREPDF=0
      REPSS=0.0
      DO4600ISET=1,NUMSET
      NI=0
      DO4700I=1,N
      IF(X(I).EQ.DUM1(ISET))NI=NI+1
      IF(X(I).EQ.DUM1(ISET))DUM2(NI)=Y(I)
 4700 CONTINUE
      ANI=NI
      SUM=0.0
      DO5100I=1,NI
      SUM=SUM+DUM2(I)
 5100 CONTINUE
      YMEAN=SUM/ANI 
      SUM=0.0
      DO5200I=1,NI
      SUM=SUM+(DUM2(I)-YMEAN)**2
 5200 CONTINUE
      IREPDF=IREPDF+NI-1
      REPSS=REPSS+SUM
 4600 CONTINUE
      IF(IREPDF.LE.0)RETURN
      REPDF=IREPDF
      REPV=REPSS/REPDF
      REPSD=SQRT(REPV)
      IF(IWRITE.NE.0)WRITE(IPR,999)
      IF(IWRITE.NE.0)WRITE(IPR,4930)REPSD
      IF(IWRITE.NE.0)WRITE(IPR,4935)IREPDF
      IF(IWRITE.NE.0)WRITE(IPR,4933)NUMSET
      IFITDF=IRESDF-IREPDF
      IF(IFITDF.GE.1)GOTO5250 
      IF(IWRITE.NE.0)WRITE(IPR,4936)
      IF(IWRITE.NE.0)WRITE(IPR,4937)
      IF(IWRITE.NE.0)WRITE(IPR,4938)
      IF(IWRITE.NE.0)WRITE(IPR,4939)
      RETURN
 5250 FITDF=IFITDF
      FITSS=RESSS-REPSS
      FITMS=FITSS/FITDF
      FSTAT=FITMS/RESMS
      CALL FCDF(FSTAT,IFITDF,IREPDF,CDF)
      CDF2=100.0*CDF
      IF(IWRITE.NE.0)WRITE(IPR,4940)FSTAT,CDF2
      IF(IWRITE.NE.0)WRITE(IPR,4945)IFITDF,IREPDF 
C
  999 FORMAT(1H )
 1505 FORMAT(1H ,6E15.8)
 1510 FORMAT(1H ,51H*****ERROR--THE MATRIX IS SINGULAR--IT HAS IRANK = ,
     1I8,51H WHICH IS LESS THAN THE NUMBER OF COEFFICIENTS K = ,I8,6H **
     1***)
 1520 FORMAT(1H ,51H*****ERROR--THE ITERATIONS ARE NOT CONVERGING *****)
 1605 FORMAT(1H ,I8,5X,E15.8) 
 2005 FORMAT(1H ,8E15.8)
 2405 FORMAT(1H ,34HCORRELATION MATRIX OF COEFFICIENTS)
 2505 FORMAT(1H ,8E15.8)
 3205 FORMAT(1H ,44H      OBSERVED      PREDICTED      RESIDUALS)
 3305 FORMAT(1H ,8E15.8)
 3310 FORMAT(1H ,10X,30HRESIDUAL STANDARD DEVIATION = ,E15.8)
 3311 FORMAT(1H ,10X,30HRESIDUAL DEGREES OF FREEDOM = ,I8)
 3312 FORMAT(1H ,10X,13HCOEFFICIENTS:)
 3315 FORMAT(1H ,58H          J        B(J)      SD(B(J))        B(J)/SD
     1(B(J)))
 3405 FORMAT(1H ,3X,I8,8E15.8)
  101 FORMAT(1H ,28HLEAST SQUARES POLYNOMIAL FIT) 
  102 FORMAT(1H ,10X,16HSAMPLE SIZE N = ,I8)
  103 FORMAT(1H ,10X,9HDEGREE = ,I8)
  105 FORMAT(1H ,33H*****ERROR--THE SAMPLE SIZE N (= ,I8,40H) IS NON-POS
     1ITIVE OR LARGER THAN NMAX = ,I8,6H *****)
  110 FORMAT(1H ,52H*****ERROR--THE DESIRED NUMBER OF COEFFICIENTS K (=
     1,I8,40H) IS NON-POSITIVE OR LARGER THAN KMAX = ,I8,6H *****)
  115 FORMAT(1H ,52H*****ERROR--THE DESIRED NUMBER OF COEFFICIENTS K (=
     1,I8,38H) IS LARGER THAN THE SAMPLE SIZE N (= ,I8,7H) *****)
  851 FORMAT(1H ,10X,20HUNEQUAL WEIGHTS CASE)
  852 FORMAT(1H ,10X,18HEQUAL WEIGHTS CASE)
 4205 FORMAT(1H ,38HERROR IN POLY   SUBROUTINE--NUMSET = 0) 
 4930 FORMAT(1H ,44H          REPLICATION STANDARD DEVIATION  = ,D15.7)
 4935 FORMAT(1H ,44H          REPLICATION DEGREES OF FREEDOM  = ,I8)
 4933 FORMAT(1H ,44H          NUMBER OF DISTINCT SUBSETS      = ,I8)
 4936 FORMAT(1H ,51H          LACK OF FIT F TEST CANNOT BE DONE BECAUSE)
 4937 FORMAT(1H ,44H          HAVE ONLY 0 DEGREES OF FREEDOM IN ,        ,
     121HNUMERATOR OF F RATIO.)
 4938 FORMAT(1H ,49H          THIS HAPPENS WHEN NUMBER OF PARAMETERS ,
     16HFITTED )
 4939 FORMAT(1H ,45H          IS IDENTICAL TO NUMBER OF DISTINCT ,
     18HSUBSETS.)
 4940 FORMAT(1H ,32H          LACK OF FIT F RATIO = ,F10.4,7H = THE , 
     1F8.4,14H% POINT OF THE) 
 4945 FORMAT(1H ,30H          F DISTRIBUTION WITH ,I6,5H AND ,I6,
     119H DEGREES OF FREEDOM) 
C
      RETURN
      END