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