C.hr DSWEEP
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DSWEEP   10/4/72
C
C     PURPOSE
C     INVERT A POSITIVE DEFINITE MATRIX IN PLACE BY A DIAGONAL SWEEP.
C
C     USAGE
C     CALL DSWEEP(A,N,EPS,IER)
C
C     ARGUMENTS
C     A   - SYMMETRIC POSITIVE DEFINITE N BY N MATRIX STORED COLUMNWISE
C           (STORAGE MODE OF 0).  ON RETURN CONTAINS THE INVERSE OF A
C           STORED COLUMNWISE.
C           REAL*8
C     N   - NUMBER OF ROWS AND COLUMNS OF A.
C           INTEGER
C     EPS - INPUT CONSTANT USED AS A RELATIVE TOLERANCE IN TESTING FOR
C           DEGENERATE RANK.  A REASONABLE VALUE FOR EPS IS 1.D-13.
C           REAL*8
C     IER - ERROR PARAMETER CODED AS FOLLOWS
C           IER=0  NO ERROR, RANK OF A IS N.
C           IER.GT.0  A IS SINGULAR, RANK OF A IS N-IER.
C           INTEGER
C
C     REMARK
C     IF IER.GT.0 THEN DSWEEP RETURNS A G-INVERSE.
C
C     REFERENCE
C     SCHATZOFF, M. ET AL.  EFFICIENT CALCULATION OF ALL POSSIBLE REG-
C     RESSIONS.  TECHNOMETRICS, 10. 769-79 (NOVEMBER 1968)
C
C     PROGRAMMER
C     DR. A. RONALD GALLANT
C     DEPARTMANT OF STATISTICS
C     NORTH CAROLINA STATE UNIVERSITY
C     RALEIGH, NORTH CAROLINA  27695-8203
C
C
C
      SUBROUTINE DSWEEP(A,N,EPS,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 A(N*N)
      TOL=0.D0
      DO 5 I=1,N
      II=-N+N*I+I
      TEST=A(II)
5     IF(TEST.GT.TOL) TOL=TEST
      TOL=TOL*EPS
      IER=0
      DO 50 K=1,N
      KK=-N+N*K+K
      AKK=A(KK)
      IF(AKK.GT.TOL) GO TO 20
      DO 10 J=K,N
      KJ=-N+N*J+K
10    A(KJ)=0.D0
      IF(K.EQ.1) GO TO 16
      KLESS1=K-1
      DO 15 I=1,KLESS1
      IK=KK-I
15    A(IK)=0.D0
16    IER=IER+1
      GO TO 50
20    D=1.D0/AKK
      DO 25 I=1,N
      DO 25 J=I,N
      IF((I.EQ.K).OR.(J.EQ.K)) GO TO 25
      IJ=N*(J-1)+I
      IF(I.LT.K) AIK=A(-N+N*K+I)
      IF(I.GT.K) AIK=A(-N+N*I+K)
      IF(K.LT.J) AKJ=A(-N+N*J+K)
      IF(K.GT.J) AKJ=-A(-N+N*K+J)
      A(IJ)=A(IJ)-AIK*AKJ*D
25    CONTINUE
      DO 30 J=K,N
      KJ=-N+N*J+K
30    A(KJ)=A(KJ)*D
      IF(K.EQ.1) GO TO 36
      KLESS1=K-1
      DO 35 I=1,KLESS1
      IK=KK-I
35    A(IK)=-A(IK)*D
36    A(KK)=D
50    CONTINUE
      DO 55 I=1,N
      DO 55 J=I,N
      IF(I.EQ.J) GO TO 55
      IJ=-N+N*J+I
      JI=-N+N*I+J
      A(JI)=A(IJ)
55    CONTINUE
      RETURN
      END