C
C     ..................................................................
C
C        SUBROUTINE APFS
C
C        PURPOSE
C           PERFORM SYMMETRIC FACTORIZATION OF THE MATRIX OF THE NORMAL
C           EQUATIONS FOLLOWED BY CALCULATION OF THE LEAST SQUARES FIT
C           OPTIONALLY
C
C        USAGE
C           CALL APFS(WORK,IP,IRES,IOP,EPS,ETA,IER)
C
C        DESCRIPTION OF PARAMETERS
C           WORK  - GIVEN SYMMETRIC COEFFICIENT MATRIX, STORED
C                   COMPRESSED, I.E UPPER TRIANGULAR PART COLUMNWISE.
C                   THE GIVEN RIGHT HAND SIDE OCCUPIES THE NEXT IP
C                   LOCATIONS IN WORK. THE VERY LAST COMPONENT OF WORK
C                   CONTAINS THE SQUARE SUM OF FUNCTION VALUES E0
C                   THIS SCHEME OF STORAGE ALLOCATION IS PRODUCED E.G.
C                   BY SUBROUTINE APLL.
C                   THE GIVEN MATRIX IS FACTORED IN THE FORM
C                   TRANSPOSE(T)*T AND THE GIVEN RIGHT HAND SIDE IS
C                   DIVIDED BY TRANSPOSE(T).
C                   THE UPPER TRIANGULAR FACTOR T IS RETURNED IN WORK IF
C                   IOP EQUALS ZERO.
C                   IN CASE OF NONZERO IOP THE CALCULATED SOLUTIONS ARE
C                   STORED IN THE COLUMNS OF TRIANGULAR ARRAY WORK OF
C                   CORRESPONDING DIMENSION AND E0  IS REPLACED BY THE
C                   SQUARE SUM OF THE ERRORS FOR FIT OF DIMENSION IRES.
C                   THE TOTAL DIMENSION OF WORK IS (IP+1)*(IP+2)/2
C           IP    - NUMBER OF FUNDAMENTAL FUNCTIONS USED FOR LEAST
C                   SQUARES FIT
C           IRES  - DIMENSION OF CALCULATED LEAST SQUARES FIT.
C                   LET N1, N2, DENOTE THE FOLLOWING NUMBERS
C                   N1 = MAXIMAL DIMENSION FOR WHICH NO LOSS OF
C                        SIGNIFICANCE WAS INDICATED DURING FACTORIZATION
C                   N2 = SMALLEST DIMENSION FOR WHICH THE SQUARE SUM OF
C                        THE ERRORS DOES NOT EXCEED TEST=ABS(ETA*FSQ)
C                   THEN IRES=MINO(IP,N1) IF IOP IS NONNEGATIVE
C                   AND  IRES=MINO(IP,N1,N2) IF IOP IS NEGATIVE
C           IOP   - INPUT PARAMETER FOR SELECTION OF OPERATION
C                   IOP = 0 MEANS TRIANGULAR FACTORIZATION, DIVISION OF
C                           THE RIGHT HAND SIDE BY TRANSPOSE(T) AND
C                           CALCULATION OF THE SQUARE SUM OF ERRORS IS
C                           PERFORMED ONLY
C                   IOP = +1 OR -1 MEANS THE SOLUTION OF DIMENSION IRES
C                           IS CALCULATED ADDITIONALLY
C                   IOP = +2 OR -2 MEANS ALL SOLUTIONS FOR DIMENSION ONE
C                           UP TO IRES ARE CALCULATED ADDITIONALLY
C           EPS   - RELATIVE TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
C                   A SENSIBLE VALUE IS BETWEEN 1.E-3 AND 1.E-6
C           ETA   - RELATIVE TOLERANCE FOR TOLERATED SQUARE SUM OF
C                   ERRORS. A REALISTIC VALUE IS BETWEEN 1.E0 AND 1.E-6
C           IER   - RESULTANT ERROR PARAMETER
C                   IER =-1 MEANS NONPOSITIVE IP
C                   IER = 0 MEANS NO LOSS OF SIGNIFICANCE DETECTED
C                           AND SPECIFIED TOLERANCE OF ERRORS REACHED
C                   IER = 1 MEANS LOSS OF SIGNIFICANCE DETECTED OR
C                           SPECIFIED TOLERANCE OF ERRORS NOT REACHED
C
C        REMARKS
C           THE ABSOLUTE TOLERANCE USED INTERNALLY FOR TEST ON LOSS OF
C           SIGNIFICANCE IS TOL=ABS(EPS*WORK(1)).
C           THE ABSOLUTE TOLERANCE USED INTERNALLY FOR THE SQUARE SUM OF
C           ERRORS IS ABS(ETA*FSQ).
C           IOP GREATER THAN 2 HAS THE SAME EFFECT AS IOP = 2.
C           IOP LESS THAN -2 HAS THE SAME EFFECT AS IOP =-2.
C           IRES = 0 MEANS THE ABSOLUTE VALUE OF EPS IS NOT LESS THAN
C           ONE AND/OR WORK(1) IS NOT POSITIVE AND/OR IP IS NOT POSITIVE
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           CALCULATION OF THE LEAST SQUARES FITS IS DONE USING
C           CHOLESKYS SQUARE ROOT METHOD FOR SYMMETRIC FACTORIZATION.
C           THE INCORPORATED TEST ON LOSS OF SIGNIFICANCE MEANS EACH
C           RADICAND MUST BE GREATER THAN THE INTERNAL ABSOLUTE
C           TOLERANCE TOL=ABS(EPS*WORK(1)).
C           IN CASE OF LOSS OF SIGNIFICANCE IN THE ABOVE SENSE ONLY A
C           SUBSYSTEM OF THE NORMAL EQUATIONS IS SOLVED.
C           IN CASE OF NEGATIVE IOP THE TRIANGULAR FACTORIZATION IS
C           TERMINATED PREMATURELY EITHER IF THE SQUARE SUM OF THE
C           ERRORS DOES NOT EXCEED ETA*FSQ OR IF THERE IS INDICATION
C           FOR LOSS OF SIGNIFICANCE
C
C     ..................................................................
C
      SUBROUTINE APFS(WORK,IP,IRES,IOP,EPS,ETA,IER)
C
C
C        DIMENSIONED DUMMY VARIABLES
      DIMENSION WORK(1)
      IRES=0
C
C        TEST OF SPECIFIED DIMENSION
      IF(IP)1,1,2
C
C        ERROR RETURN IN CASE OF ILLEGAL DIMENSION
    1 IER=-1
      RETURN
C
C        INITIALIZE FACTORIZATION PROCESS
    2 IPIV=0
      IPP1=IP+1
      IER=1
      ITE=IP*IPP1/2
      IEND=ITE+IPP1
      TOL=ABS(EPS*WORK(1))
      TEST=ABS(ETA*WORK(IEND))
C
C        START LOOP OVER ALL ROWS OF WORK
      DO 11 I=1,IP
      IPIV=IPIV+I
      JA=IPIV-IRES
      JE=IPIV-1
C
C        FORM SCALAR PRODUCT NEEDED TO MODIFY CURRENT ROW ELEMENTS
      JK=IPIV
      DO 9 K=I,IPP1
      SUM=0.
      IF(IRES)5,5,3
    3 JK=JK-IRES
      DO 4 J=JA,JE
      SUM=SUM+WORK(J)*WORK(JK)
    4 JK=JK+1
    5 IF(JK-IPIV)6,6,8
C
C        TEST FOR LOSS OF SIGNIFICANCE
    6 SUM=WORK(IPIV)-SUM
      IF(SUM-TOL)12,12,7
    7 SUM=SQRT(SUM)
      WORK(IPIV)=SUM
      PIV=1./SUM
      GOTO 9
C
C        UPDATE OFF-DIAGONAL TERMS
    8 SUM=(WORK(JK)-SUM)*PIV
      WORK(JK)=SUM
    9 JK=JK+K
C
C        UPDATE SQUARE SUM OF ERRORS
      WORK(IEND)=WORK(IEND)-SUM*SUM
C
C        RECORD ADDRESS OF LAST PIVOT ELEMENT
      IRES=IRES+1
      IADR=IPIV
C
C        TEST FOR TOLERABLE ERROR IF SPECIFIED
      IF(IOP)10,11,11
   10 IF(WORK(IEND)-TEST)13,13,11
   11 CONTINUE
      IF(IOP)12,22,12
C
C        PERFORM BACK SUBSTITUTION IF SPECIFIED
   12 IF(IOP)14,23,14
   13 IER=0
   14 IPIV=IRES
   15 IF(IPIV)23,23,16
   16 SUM=0.
      JA=ITE+IPIV
      JJ=IADR
      JK=IADR
      K=IPIV
      DO 19 I=1,IPIV
      WORK(JK)=(WORK(JA)-SUM)/WORK(JJ)
      IF(K-1)20,20,17
   17 JE=JJ-1
      SUM=0.
      DO 18 J=K,IPIV
      SUM=SUM+WORK(JK)*WORK(JE)
      JK=JK+1
   18 JE=JE+J
      JK=JE-IPIV
      JA=JA-1
      JJ=JJ-K
   19 K=K-1
   20 IF(IOP/2)21,23,21
   21 IADR=IADR-IPIV
      IPIV=IPIV-1
      GOTO 15
C
C        NORMAL RETURN
   22 IER=0
   23 RETURN
      END
