C....*...1.........2.........3.........4.........5.........6.........7.*.......8
C     REGR8    7/25/73
C
C     PURPOSE
C     MULTIPLE REGRESSION
C
C     USAGE
C     CALL REGR8(Y,X,B,E,N,K,G,VAR,IER)
C
C     ARGUMENTS
C     Y   - AN N BY 1 INPUT VECTOR
C           ELEMENTS OF Y ARE REAL*8
C     X   - AN N BY K INPUT MATRIX STORED COLUMNWISE (STORAGE MODE OF 0)
C           ELEMENTS OF X ARE REAL*8
C     B   - A  K BY 1 VECTOR OF REGRESSION COEFFICIENTS
C           ELEMENTS OF B ARE REAL*8
C     E   - AN N BY 1 VECTOR OF RESIDUALS
C           ELEMENTS OF E ARE REAL*8
C     N   - NUMBER OF ROWS IN X, LENGTH OF Y AND E
C           INTEGER
C     K   - NUMBER OF COLUMNS IN X, LENGTH OF B, NUMBER OF ROWS AND
C           COLUMNS IN G.  N MUST BE .GE. K.
C           INTEGER
C     G   - A K BY K MATRIX CONTAINING THE ESTIMATED VARIANCE-COVARIANCE
C           MATRIX OF THE REGRESSION COEFFICIENTS.  STORED COLUMNWISE
C           (STORAGE MODE OF 0)
C           ELEMENTS OF G ARE REAL*8
C     VAR - ESTIMATED VARIANCE WITH N-K DEGREES FREEDOM (SSE/(N-K)).
C           REAL*8
C     IER - ERROR PARAMETER CODED AS FOLLOWS:
C           IER=0  NO ERROR.
C           IER>0  X IS NOT FULL RANK (RANK OF X IS K-IER).
C           INTEGER
C
C     REMARKS
C     ON RETURN:
C     B=INVERSE(X'X)X'Y
C     E=Y-XB
C     VAR=SUM(E(I)**2)/(N-K)
C     G=VAR*INVERSE(X'X)
C     THE USAGE
C     CALL REGR8(Y,X,B,Y,N,K,G,VAR,IER)
C     IS PERMISSIBLE.  Y WILL CONTAIN E ON RETURN.
C
      SUBROUTINE REGR8(Y,X,B,E,N,K,G,VAR,IER)
      implicit real*8 (a-h,o-z)
      save
      REAL*8 Y(1),X(1),B(1),E(1),VAR
      REAL*8 G(1),GDUMMY(36)
      IF(K.LE.6) GO TO 10
      CALL Z2REGR8(Y,X,B,E,N,K,G,VAR,IER)
      RETURN
10    CALL Z2REGR8(Y,X,B,E,N,K,GDUMMY,VAR,IER)
      KK=K*K
      DO 20 I=1,KK
20    G(I)=GDUMMY(I)
      RETURN
      END
      SUBROUTINE Z2REGR8(Y,X,B,E,N,K,G,VAR,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 Y(1),X(1),B(1),E(1),VAR
      REAL*8 G(1)
C     COMPUTE AND STORE THE ARRAY |X'X X'Y| IN G WITH STORAGE MODE 1.
C                                 |Y'X Y'Y|
      IBXY=(K*K+K)/2
      IAYY=IBXY+K+1
      IBS=IAYY
      KLESS1=K-1
      KPLUS1=K+1
      DO 10 I=1,IAYY
10    G(I)=0.D0
      SUM=0.D0
      DO 20 L=1,N
      YL=Y(L)
      SUM=SUM+YL*YL
      DO 20 I=1,K
      LI=N*(I-1)+L
      XLI=X(LI)
      IA=IBXY+I
      G(IA)=G(IA)+XLI*YL
      DO 20 J=I,K
      IJ=(J*J-J)/2+I
      LJ=N*(J-1)+L
      XLJ=X(LJ)
20    G(IJ)=G(IJ)+XLI*XLJ
      G(IAYY)=SUM
C     CONVERT G TO A CORRELATION MATRIX.
      DO 30 I=1,KPLUS1
      IASI=IBS+I
      II=(I*I+I)/2
      SI=DSQRT(G(II))
      IF(SI.LT.1.D-50) SI=1.D0
30    G(IASI)=1.D0/SI
      DO 35 I=1,KPLUS1
      IASI=IBS+I
      DO 35 J=I,KPLUS1
      IJ=(J*J-J)/2+I
      IASJ=IBS+J
35    G(IJ)=G(IJ)*G(IASI)*G(IASJ)
      G(IBS+K+1)=SI
C     SWEEP G.
      TOL=1.D-13
      IER=0
      DO 49 L=1,K
      LL=(L*L+L)/2
      D=G(LL)
      IF(D.GT.TOL) GO TO 44
      DO 41 J=L,KPLUS1
      LJ=(J*J-J)/2+L
41    G(LJ)=0.D0
      IF(L.EQ.1) GO TO 43
      LLESS1=L-1
      DO 42 I=1,LLESS1
      IL=LL-I
42    G(IL)=0.D0
43    IER=IER+1
      GO TO 49
44    D=1.D0/D
      DO 45 I=1,KPLUS1
      DO 45 J=I,KPLUS1
      IF((I.EQ.L).OR.(J.EQ.L)) GO TO 45
      IJ=(J*J-J)/2+I
      IF(I.LT.L) GIL=G((L*L-L)/2+I)
      IF(I.GT.L) GIL=G((I*I-I)/2+L)
      IF(L.LT.J) GLJ=G((J*J-J)/2+L)
      IF(L.GT.J) GLJ=-G((L*L-L)/2+J)
      G(IJ)=G(IJ)-GIL*GLJ*D
45    CONTINUE
      DO 46 J=L,KPLUS1
      LJ=(J*J-J)/2+L
46    G(LJ)=G(LJ)*D
      IF(L.EQ.1) GO TO 48
      LLESS1=L-1
      DO 47 I=1,LLESS1
      IL=LL-I
47    G(IL)=-G(IL)*D
48    G(LL)=D
49    CONTINUE
C     RESCALE G TO A SWEPT SUM OF SQUARES AND CROSS PRODUCTS MATRIX.
      DO 50 I=1,KPLUS1
      IASI=IBS+I
      DO 50 J=I,KPLUS1
      IJ=(J*J-J)/2+I
      IASJ=IBS+J
50    G(IJ)=G(IJ)*G(IASI)*G(IASJ)
C     COMPUTE E, VAR
      SSE=0.D0
      DO 60 I=1,N
      SUM=0.D0
      DO 55 J=1,K
      IJ=N*(J-1)+I
      XIJ=X(IJ)
      IABJ=IBXY+J
55    SUM=SUM+XIJ*G(IABJ)
      EI=Y(I)-SUM
      SSE=SSE+EI*EI
60    E(I)=EI
      IF (N.GT.K) THEN
        V=SSE/DFLOAT(N-K)
      ELSE
        V=SSE
      END IF
      VAR=V
C     ASSIGN COEFICIENTS TO B.
      DO 65 I=1,K
      IABI=IBXY+I
65    B(I)=G(IABI)
C     CONVERT G VARIANCE-COVARIANCE MATRIX - STORAGE MODE 0
      DO 70 J=1,K
      JJ=K-J+1
      DO 70 I=1,JJ
      II=JJ-I+1
      L0=K*(JJ-1)+II
      L1=II+(JJ*JJ-JJ)/2
70    G(L0)=G(L1)*V
      DO 80 J=1,K
      DO 80 I=1,J
      LU=K*(J-1)+I
      LL=K*(I-1)+J
80    G(LL)=G(LU)
      RETURN
      END