C
C     ..................................................................
C
C        SUBROUTINE MFSS
C
C        PURPOSE
C           GIVEN A SYMMETRIC POSITIVE SEMI DEFINITE MATRIX , MFSS WILL
C           (1) DETERMINE THE RANK AND LINEARLY INDEPENDENT ROWS AND
C               COLUMNS
C           (2) FACTOR A SYMMETRIC SUBMATRIX OF MAXIMAL RANK
C           (3) EXPRESS NONBASIC ROWS IN TERMS OF BASIC ONES,
C               EXPRESS NONBASIC COLUMNS IN TERMS OF BASIC ONES
C               EXPRESS BASIC VARIABLES IN TERMS OF FREE ONES
C           SUBROUTINE MFSS MAY BE USED AS A PREPARATORY STEP FOR THE
C           CALCULATION OF THE LEAST SQUARES SOLUTION OF MINIMAL
C           LENGTH OF A SYSTEM OF LINEAR EQUATIONS WITH SYMMETRIC
C           POSITIVE SEMI-DEFINITE COEFFICIENT MATRIX
C
C        USAGE
C           CALL MFSS(A,N,EPS,IRANK,TRAC)
C
C        DESCRIPTION OF PARAMETERS
C           A     - UPPER TRIANGULAR PART OF GIVEN SYMMETRIC SEMI-
C                   DEFINITE MATRIX STORED COLUMNWISE IN COMPRESSED FORM
C                   ON RETURN A CONTAINS THE MATRIX T AND, IF IRANK IS
C                   LESS THAN N, THE MATRICES U AND TU
C           N     - DIMENSION OF GIVEN MATRIX A
C           EPS   - TESTVALUE FOR ZERO AFFECTED BY ROUND-OFF NOISE
C           IRANK - RESULTANT VARIABLE, CONTAINING THE RANK OF GIVEN
C                   MATRIX A IF A IS SEMI-DEFINITE
C                   IRANK = 0 MEANS A HAS NO POSITIVE DIAGONAL ELEMENT
C                             AND/OR EPS IS NOT ABSOLUTELY LESS THAN ONE
C                   IRANK =-1 MEANS DIMENSION N IS NOT POSITIVE
C                   IRANK =-2 MEANS COMPLETE FAILURE, POSSIBLY DUE TO
C                             INADEQUATE RELATIVE TOLERANCE EPS
C           TRAC  - VECTOR OF DIMENSION N CONTAINING THE
C                   SOURCE INDEX OF THE I-TH PIVOT ROW IN ITS I-TH
C                   LOCATION, THIS MEANS THAT TRAC CONTAINS THE
C                   PRODUCT REPRESENTATION OF THE PERMUTATION WHICH
C                   IS APPLIED TO ROWS AND COLUMNS OF A IN TERMS OF
C                   TRANSPOSITIONS
C
C        REMARKS
C           EPS MUST BE ABSOLUTELY LESS THAN ONE. A SENSIBLE VALUE IS
C           SOMEWHERE IN BETWEEN 10**(-4) AND 10**(-6)
C           THE ABSOLUTE VALUE OF INPUT PARAMETER EPS IS USED AS
C           RELATIVE TOLERANCE.
C           IN ORDER TO PRESERVE SYMMETRY ONLY PIVOTING ALONG THE
C           DIAGONAL IS BUILT IN.
C           ALL PIVOTELEMENTS MUST BE GREATER THAN THE ABSOLUTE VALUE
C           OF EPS TIMES ORIGINAL DIAGONAL ELEMENT
C           OTHERWISE THEY ARE TREATED AS IF THEY WERE ZERO
C           MATRIX A REMAINS UNCHANGED IF THE RESULTANT VALUE IRANK
C           EQUALS ZERO
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           THE SQUARE ROOT METHOD WITH DIAGONAL PIVOTING IS USED FOR
C           CALCULATION OF THE RIGHT HAND TRIANGULAR FACTOR.
C           IN CASE OF AN ONLY SEMI-DEFINITE MATRIX THE SUBROUTINE
C           RETURNS THE IRANK X IRANK UPPER TRIANGULAR FACTOR T OF A
C           SUBMATRIX OF MAXIMAL RANK, THE IRANK X (N-IRANK) MATRIX U
C           AND THE (N-IRANK) X (N-IRANK) UPPER TRIANGULAR TU SUCH
C           THAT TRANSPOSE(TU)*TU=I+TRANSPOSE(U)*U
C
C     ..................................................................
C
      SUBROUTINE MFSS(A,N,EPS,IRANK,TRAC)
C
C
C        DIMENSIONED DUMMY VARIABLES
      DIMENSION A(1),TRAC(1)
      DOUBLE PRECISION SUM
C
C        TEST OF SPECIFIED DIMENSION
      IF(N)36,36,1
C
C        INITIALIZE TRIANGULAR FACTORIZATION
    1 IRANK=0
      ISUB=0
      KPIV=0
      J=0
      PIV=0.
C
C        SEARCH FIRST PIVOT ELEMENT
      DO 3 K=1,N
      J=J+K
      TRAC(K)=A(J)
      IF(A(J)-PIV)3,3,2
    2 PIV=A(J)
      KSUB=J
      KPIV=K
    3 CONTINUE
C
C        START LOOP OVER ALL ROWS OF A
      DO 32 I=1,N
      ISUB=ISUB+I
      IM1=I-1
    4 KMI=KPIV-I
      IF(KMI)35,9,5
C
C        PERFORM PARTIAL COLUMN INTERCHANGE
    5 JI=KSUB-KMI
      IDC=JI-ISUB
      JJ=ISUB-IM1
      DO 6 K=JJ,ISUB
      KK=K+IDC
      HOLD=A(K)
      A(K)=A(KK)
    6 A(KK)=HOLD
C
C        PERFORM PARTIAL ROW INTERCHANGE
      KK=KSUB
      DO 7 K=KPIV,N
      II=KK-KMI
      HOLD=A(KK)
      A(KK)=A(II)
      A(II)=HOLD
    7 KK=KK+K
C
C        PERFORM REMAINING INTERCHANGE
      JJ=KPIV-1
      II=ISUB
      DO 8 K=I,JJ
      HOLD=A(II)
      A(II)=A(JI)
      A(JI)=HOLD
      II=II+K
    8 JI=JI+1
    9 IF(IRANK)22,10,10
C
C        RECORD INTERCHANGE IN TRANSPOSITION VECTOR
   10 TRAC(KPIV)=TRAC(I)
      TRAC(I)=KPIV
C
C        MODIFY CURRENT PIVOT ROW
      KK=IM1-IRANK
      KMI=ISUB-KK
      PIV=0.
      IDC=IRANK+1
      JI=ISUB-1
      JK=KMI
      JJ=ISUB-I
      DO 19 K=I,N
      SUM=0.D0
C
C        BUILD UP SCALAR PRODUCT IF NECESSARY
      IF(KK)13,13,11
   11 DO 12 J=KMI,JI
      SUM=SUM-A(J)*A(JK)
   12 JK=JK+1
   13 JJ=JJ+K
      IF(K-I)14,14,16
   14 SUM=A(ISUB)+SUM
C
C        TEST RADICAND FOR LOSS OF SIGNIFICANCE
      IF(SUM-ABS(A(ISUB)*EPS))20,20,15
   15 A(ISUB)=DSQRT(SUM)
      KPIV=I+1
      GOTO 19
   16 SUM=(A(JK)+SUM)/A(ISUB)
      A(JK)=SUM
C
C        SEARCH FOR NEXT PIVOT ROW
      IF(A(JJ))19,19,17
   17 TRAC(K)=TRAC(K)-SUM*SUM
      HOLD=TRAC(K)/A(JJ)
      IF(PIV-HOLD)18,19,19
   18 PIV=HOLD
      KPIV=K
      KSUB=JJ
   19 JK=JJ+IDC
      GOTO 32
C
C        CALCULATE MATRIX OF DEPENDENCIES U
   20 IF(IRANK)21,21,37
   21 IRANK=-1
      GOTO 4
   22 IRANK=IM1
      II=ISUB-IRANK
      JI=II
      DO 26 K=1,IRANK
      JI=JI-1
      JK=ISUB-1
      JJ=K-1
      DO 26 J=I,N
      IDC=IRANK
      SUM=0.D0
      KMI=JI
      KK=JK
      IF(JJ)25,25,23
   23 DO 24 L=1,JJ
      IDC=IDC-1
      SUM=SUM-A(KMI)*A(KK)
      KMI=KMI-IDC
   24 KK=KK-1
   25 A(KK)=(SUM+A(KK))/A(KMI)
   26 JK=JK+J
C
C        CALCULATE I+TRANSPOSE(U)*U
      JJ=ISUB-I
      PIV=0.
      KK=ISUB-1
      DO 31 K=I,N
      JJ=JJ+K
      IDC=0
      DO 28 J=K,N
      SUM=0.D0
      KMI=JJ+IDC
      DO 27 L=II,KK
      JK=L+IDC
   27 SUM=SUM+A(L)*A(JK)
      A(KMI)=SUM
   28 IDC=IDC+J
      A(JJ)=A(JJ)+1.D0
      TRAC(K)=A(JJ)
C
C        SEARCH NEXT DIAGONAL ELEMENT
      IF(PIV-A(JJ))29,30,30
   29 KPIV=K
      KSUB=JJ
      PIV=A(JJ)
   30 II=II+K
      KK=KK+K
   31 CONTINUE
      GOTO 4
   32 CONTINUE
   33 IF(IRANK)35,34,35
   34 IRANK=N
   35 RETURN
C
C        ERROR RETURNS
C
C        RETURN IN CASE OF ILLEGAL DIMENSION
   36 IRANK=-1
      RETURN
C
C        INSTABLE FACTORIZATION OF I+TRANSPOSE(U)*U
   37 IRANK=-2
      RETURN
      END
