C
C     ..................................................................
C
C        SUBROUTINE DSINV
C
C        PURPOSE
C           INVERT A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX
C
C        USAGE
C           CALL DSINV(A,N,EPS,IER)
C
C        DESCRIPTION OF PARAMETERS
C           A      - DOUBLE PRECISION UPPER TRIANGULAR PART OF GIVEN
C                    SYMMETRIC POSITIVE DEFINITE N BY N COEFFICIENT
C                    MATRIX.
C                    ON RETURN A CONTAINS THE RESULTANT UPPER
C                    TRIANGULAR MATRIX IN DOUBLE PRECISION.
C           N      - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX.
C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED
C                    AS RELATIVE TOLERANCE FOR TEST ON LOSS OF
C                    SIGNIFICANCE.
C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
C                    IER=0  - NO ERROR
C                    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-
C                             TER N OR BECAUSE SOME RADICAND IS NON-
C                             POSITIVE (MATRIX A IS NOT POSITIVE
C                             DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI-
C                             FICANCE)
C                    IER=K  - WARNING WHICH INDICATES LOSS OF SIGNIFI-
C                             CANCE. THE RADICAND FORMED AT FACTORIZA-
C                             TION STEP K+1 WAS STILL POSITIVE BUT NO
C                             LONGER GREATER THAN ABS(EPS*A(K+1,K+1)).
C
C        REMARKS
C           THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE
C           STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS.
C           IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU-
C           LAR MATRIX IS STORED COLUMNWISE TOO.
C           THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL
C           CALCULATED RADICANDS ARE POSITIVE.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           DMFSD
C
C        METHOD
C           SOLUTION IS DONE USING FACTORIZATION BY SUBROUTINE DMFSD.
C
C     ..................................................................
C
      SUBROUTINE DSINV(A,N,EPS,IER)
C
C
      DIMENSION A(1)
      DOUBLE PRECISION A,DIN,WORK
C
C        FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE DMFSD
C        A = TRANSPOSE(T) * T
      CALL DMFSD(A,N,EPS,IER)
      IF(IER) 9,1,1
C
C        INVERT UPPER TRIANGULAR MATRIX T
C        PREPARE INVERSION-LOOP
    1 IPIV=N*(N+1)/2
      IND=IPIV
C
C        INITIALIZE INVERSION-LOOP
      DO 6 I=1,N
      DIN=1.D0/A(IPIV)
      A(IPIV)=DIN
      MIN=N
      KEND=I-1
      LANF=N-KEND
      IF(KEND) 5,5,2
    2 J=IND
C
C        INITIALIZE ROW-LOOP
      DO 4 K=1,KEND
      WORK=0.D0
      MIN=MIN-1
      LHOR=IPIV
      LVER=J
C
C        START INNER LOOP
      DO 3 L=LANF,MIN
      LVER=LVER+1
      LHOR=LHOR+L
    3 WORK=WORK+A(LVER)*A(LHOR)
C        END OF INNER LOOP
C
      A(J)=-WORK*DIN
    4 J=J-MIN
C        END OF ROW-LOOP
C
    5 IPIV=IPIV-MIN
    6 IND=IND-1
C        END OF INVERSION-LOOP
C
C        CALCULATE INVERSE(A) BY MEANS OF INVERSE(T)
C        INVERSE(A) = INVERSE(T) * TRANSPOSE(INVERSE(T))
C        INITIALIZE MULTIPLICATION-LOOP
      DO 8 I=1,N
      IPIV=IPIV+I
      J=IPIV
C
C        INITIALIZE ROW-LOOP
      DO 8 K=I,N
      WORK=0.D0
      LHOR=J
C
C        START INNER LOOP
      DO 7 L=K,N
      LVER=LHOR+K-I
      WORK=WORK+A(LHOR)*A(LVER)
    7 LHOR=LHOR+L
C        END OF INNER LOOP
C
      A(J)=WORK
    8 J=J+K
C        END OF ROW- AND MULTIPLICATION-LOOP
C
    9 RETURN
      END