C
C     ..................................................................
C
C        SUBROUTINE GDATA
C
C        PURPOSE
C           GENERATE INDEPENDENT VARIABLES UP TO THE M-TH POWER (THE
C           HIGHEST DEGREE POLYNOMIAL SPECIFIED) AND COMPUTE MEANS,
C           STANDARD DEVIATIONS, AND CORRELATION COEFFICIENTS.  THIS
C           SUBROUTINE IS NORMALLY CALLED BEFORE SUBROUTINES ORDER,
C           MINV AND MULTR IN THE PERFORMANCE OF A POLYNOMIAL
C           REGRESSION.
C
C        USAGE
C           CALL GDATA (N,M,X,XBAR,STD,D,SUMSQ)
C
C        DESCRIPTION OF PARAMETERS
C           N     - NUMBER OF OBSERVATIONS.
C           M     - THE HIGHEST DEGREE POLYNOMIAL TO BE FITTED.
C           X     - INPUT MATRIX (N BY M+1) .  WHEN THE SUBROUTINE IS
C                   CALLED, DATA FOR THE INDEPENDENT VARIABLE ARE
C                   STORED IN THE FIRST COLUMN OF MATRIX X, AND DATA FOR
C                   THE DEPENDENT VARIABLE ARE STORED IN THE LAST
C                   COLUMN OF THE MATRIX.  UPON RETURNING TO THE
C                   CALLING ROUTINE, GENERATED POWERS OF THE INDEPENDENT
C                   VARIABLE ARE STORED IN COLUMNS 2 THROUGH M.
C           XBAR  - OUTPUT VECTOR OF LENGTH M+1 CONTAINING MEANS OF
C                   INDEPENDENT AND DEPENDENT VARIABLES.
C           STD   - OUTPUT VECTOR OF LENGTH M+1 CONTAINING STANDARD
C                   DEVIATIONS OF INDEPENDENT AND DEPENDENT VARIABLES.
C           D     - OUTPUT MATRIX (ONLY UPPER TRIANGULAR PORTION OF THE
C                   SYMMETRIC MATRIX OF M+1 BY M+1) CONTAINING CORRELA-
C                   TION COEFFICIENTS.  (STORAGE MODE OF 1)
C           SUMSQ - OUTPUT VECTOR OF LENGTH M+1 CONTAINING SUMS OF
C                   PRODUCTS OF DEVIATIONS FROM MEANS  OF INDEPENDENT
C                   AND DEPENDENT VARIABLES.
C
C        REMARKS
C           N MUST BE GREATER THAN M+1.
C           IF M IS EQUAL TO 5 OR GREATER, SINGLE PRECISION MAY NOT BE
C           SUFFICIENT TO GIVE SATISFACTORY COMPUTATIONAL RESULTS.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           REFER TO B. OSTLE, 'STATISTICS IN RESEARCH', THE IOWA STATE
C           COLLEGE PRESS, 1954, CHAPTER 6.
C
C     ..................................................................
C
      SUBROUTINE GDATA (N,M,X,XBAR,STD,D,SUMSQ)
      DIMENSION X(1),XBAR(1),STD(1),D(1),SUMSQ(1)
C
C        ...............................................................
C
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C        STATEMENT WHICH FOLLOWS.
C
C     DOUBLE PRECISION X,XBAR,STD,D,SUMSQ,T1,T2
C
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C        ROUTINE.
C
C        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
C        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  SQRT AND ABS IN
C        STATEMENT 180 MUST BE CHANGED TO DSQRT AND DABS.
C
C        ...............................................................
C
C     GENERATE INDEPENDENT VARIABLES
C
      IF(M-1) 105, 105, 90
   90 L1=0
      DO 100 I=2,M
      L1=L1+N
      DO 100 J=1,N
      L=L1+J
      K=L-N
  100 X(L)=X(K)*X(J)
C
C     CALCULATE MEANS
C
  105 MM=M+1
      DF=N
      L=0
      DO 115 I=1,MM
      XBAR(I)=0.0
      DO 110 J=1,N
      L=L+1
  110 XBAR(I)=XBAR(I)+X(L)
  115 XBAR(I)=XBAR(I)/DF
C
      DO 130 I=1,MM
  130 STD(I)=0.0
C
C     CALCULATE SUMS OF CROSS-PRODUCTS OF DEVIATIONS
C
      L=((MM+1)*MM)/2
      DO 150 I=1,L
  150 D(I)=0.0
      DO 170 K=1,N
      L=0
      DO 170 J=1,MM
      L2=N*(J-1)+K
      T2=X(L2)-XBAR(J)
      STD(J)=STD(J)+T2
      DO 170 I=1,J
      L1=N*(I-1)+K
      T1=X(L1)-XBAR(I)
      L=L+1
  170 D(L)=D(L)+T1*T2
      L=0
      DO 175 J=1,MM
      DO 175 I=1,J
      L=L+1
  175 D(L)=D(L)-STD(I)*STD(J)/DF
      L=0
      DO 180 I=1,MM
      L=L+I
      SUMSQ(I)=D(L)
  180 STD(I)= SQRT( ABS(D(L)))
C
C     CALCULATE CORRELATION COEFFICIENTS
C
      L=0
      DO 190 J=1,MM
      DO 190 I=1,J
      L=L+1
  190 D(L)=D(L)/(STD(I)*STD(J))
C
C     CALCULATE STANDARD DEVIATIONS
C
      DF=SQRT(DF-1.0)
      DO 200 I=1,MM
  200 STD(I)=STD(I)/DF
      RETURN
      END