C
C SUBROUTINE 'MFSD' OF THE IBM SCIENTIFIC SUBROUTINE PACKAGE.
C
C NOTE:  TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS
C        (1) HAVE BEEN CHANGED TO (*).
C
C     ..................................................................
C
C        SUBROUTINE MFSD
C
C        PURPOSE
C           FACTOR A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX
C
C        USAGE
C           CALL MFSD(A,N,EPS,IER)
C
C        DESCRIPTION OF PARAMETERS
C           A      - UPPER TRIANGULAR PART OF THE GIVEN SYMMETRIC
C                    POSITIVE DEFINITE N BY N COEFFICIENT MATRIX.
C                    ON RETURN A CONTAINS THE RESULTANT UPPER
C                    TRIANGULAR MATRIX.
C           N      - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX.
C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
C                    TOLERANCE FOR TEST ON LOSS OF 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           THE PRODUCT OF RETURNED DIAGONAL TERMS IS EQUAL TO THE
C           SQUARE-ROOT OF THE DETERMINANT OF THE GIVEN MATRIX.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           SOLUTION IS DONE USING THE SQUARE-ROOT METHOD OF CHOLESKY.
C           THE GIVEN MATRIX IS REPRESENTED AS PRODUCT OF TWO TRIANGULAR
C           MATRICES, WHERE THE LEFT HAND FACTOR IS THE TRANSPOSE OF
C           THE RETURNED RIGHT HAND FACTOR.
C
C     ..................................................................
C
      SUBROUTINE MFSD(A,N,EPS,IER)
C
C
      DIMENSION A(*)
      DOUBLE PRECISION DPIV,DSUM
C
C        TEST ON WRONG INPUT PARAMETER N
      IF(N-1) 12,1,1
    1 IER=0
C
C        INITIALIZE DIAGONAL-LOOP
      KPIV=0
      DO 11 K=1,N
      KPIV=KPIV+K
      IND=KPIV
      LEND=K-1
C
C        CALCULATE TOLERANCE
      TOL=ABS(EPS*A(KPIV))
C
C        START FACTORIZATION-LOOP OVER K-TH ROW
      DO 11 I=K,N
      DSUM=0.D0
      IF(LEND) 2,4,2
C
C        START INNER LOOP
    2 DO 3 L=1,LEND
      LANF=KPIV-L
      LIND=IND-L
    3 DSUM=DSUM+DBLE(A(LANF)*A(LIND))
C        END OF INNER LOOP
C
C        TRANSFORM ELEMENT A(IND)
    4 DSUM=DBLE(A(IND))-DSUM
      IF(I-K) 10,5,10
C
C        TEST FOR NEGATIVE PIVOT ELEMENT AND FOR LOSS OF SIGNIFICANCE
    5 IF(SNGL(DSUM)-TOL) 6,6,9
    6 IF(DSUM) 12,12,7
    7 IF(IER) 8,8,9
    8 IER=K-1
C
C        COMPUTE PIVOT ELEMENT
    9 DPIV=DSQRT(DSUM)
      A(KPIV)=DPIV
      DPIV=1.D0/DPIV
      GO TO 11
C
C        CALCULATE TERMS IN ROW
   10 A(IND)=DSUM*DPIV
   11 IND=IND+I
C
C        END OF DIAGONAL-LOOP
      RETURN
   12 IER=-1
      RETURN
      END
C
C=======================================================================
C