C.hr DMFSD
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     SUBROUTINE DMFSD
C
C     PURPOSE
C     FACTOR A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX
C
C     USAGE
C     CALL DMFSD(A,N,EPS,IER)
C
C     ARGUMENTS
C     A   - UPPER TRIANGULAR PART OF GIVEN SYMMETRIC POSITIVE DEFINITE N
C           BY N MATRIX; ON RETURN, UPPER TRIANGULAR R, A=R'R.
C           REAL*8
C     N   - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX
C           INTEGER*4
C     EPS - INPUT CONSTANT WHICH IS USED AS RELATIVE TOLERANCE FOR TEST
C           ON LOSS OF SIGNIFICANCE, A REASONABLE VALUE IS 1.D-13.
C           REAL*8
C     IER - ERROR PARAMETER, CODED AS FOLLOWS
C           IER=0  - NO ERROR
C           IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAMETER N OR
C                    BECAUSE SOME RADICAND IS NON-POSITIVE (MATRIX A IS
C                    NOT POSITIVE DEFINITE, POSSIBLY DUE TO LOSS OF
C                    SIGNIFICANCE)
C           IER=K  - WARNING WHICH INDICATES LOSS OF SIGNIFICANCE.  THE
C                    RADICAND FORMED AT FACTORIZATION STEP K+1 WAS STILL
C                    POSITIVE BUT NO LONGER GREATER THAN
C                    ABS(EPS*A(K+1,K+1)
C
C     REMARK
C     THIS IS A COPY OF SSP ROUTINE DMFSD
C
      SUBROUTINE DMFSD(A,N,EPS,IER)
      save
      REAL*8 EPS,DPIV,DSUM
      REAL*8 A((N*N+N)/2)
      IF(N-1) 12,1,1
    1 IER=0
      KPIV=0
      DO 11 K=1,N
      KPIV=KPIV+K
      IND=KPIV
      LEND=K-1
      TOL=ABS(SNGL(EPS*A(KPIV)))
      DO 11 I=K,N
      DSUM=0.D0
      IF(LEND) 2,4,2
    2 DO 3 L=1,LEND
      LANF=KPIV-L
      LIND=IND-L
    3 DSUM=DSUM+A(LANF)*A(LIND)
    4 DSUM=A(IND)-DSUM
      IF(I-K) 10,5,10
    5 IF(SNGL(DSUM)-TOL) 6,6,9
    6 IF(DSUM) 12,12,7
    7 IF(IER) 8,8,9
    8 IER=K-1
    9 DPIV=DSQRT(DSUM)
      A(KPIV)=DPIV
      DPIV=1.D0/DPIV
      GO TO 11
   10 A(IND)=DSUM*DPIV
   11 IND=IND+I
      RETURN
   12 IER=-1
      RETURN
      END