C DAPF 10
C ..................................................................DAPF 20
C DAPF 30
C SUBROUTINE DAPFS DAPF 40
C DAPF 50
C PURPOSE DAPF 60
C PERFORM SYMMETRIC FACTORIZATION OF THE MATRIX OF THE NORMAL DAPF 70
C EQUATIONS FOLLOWED BY CALCULATION OF THE LEAST SQUARES FIT DAPF 80
C OPTIONALLY DAPF 90
C DAPF 100
C USAGE DAPF 110
C CALL DAPFS(WORK,IP,IRES,IOP,EPS,ETA,IER) DAPF 120
C DAPF 130
C DESCRIPTION OF PARAMETERS DAPF 140
C WORK - GIVEN SYMMETRIC COEFFICIENT MATRIX, STORED DAPF 150
C COMPRESSED, I.E UPPER TRIANGULAR PART COLUMNWISE. DAPF 160
C THE GIVEN RIGHT HAND SIDE OCCUPIES THE NEXT IP DAPF 170
C LOCATIONS IN WORK. THE VERY LAST COMPONENT OF WORK DAPF 180
C CONTAINS THE SQUARE SUM OF FUNCTION VALUES E0 DAPF 190
C THIS SCHEME OF STORAGE ALLOCATION IS PRODUCED E.G. DAPF 200
C BY SUBROUTINE APLL. DAPF 210
C THE GIVEN MATRIX IS FACTORED IN THE FORM DAPF 220
C TRANSPOSE(T)*T AND THE GIVEN RIGHT HAND SIDE IS DAPF 230
C DIVIDED BY TRANSPOSE(T). DAPF 240
C THE UPPER TRIANGULAR FACTOR T IS RETURNED IN WORK IFDAPF 250
C IOP EQUALS ZERO. DAPF 260
C IN CASE OF NONZERO IOP THE CALCULATED SOLUTIONS ARE DAPF 270
C STORED IN THE COLUMNS OF TRIANGULAR ARRAY WORK OF DAPF 280
C CORRESPONDING DIMENSION AND E0 IS REPLACED BY THE DAPF 290
C SQUARE SUM OF THE ERRORS FOR FIT OF DIMENSION IRES. DAPF 300
C THE TOTAL DIMENSION OF WORK IS (IP+1)*(IP+2)/2 DAPF 310
C WORK MUST BE OF DOUBLE PRECISION DAPF 320
C IP - NUMBER OF FUNDAMENTAL FUNCTIONS USED FOR LEAST DAPF 330
C SQUARES FIT DAPF 340
C IRES - DIMENSION OF CALCULATED LEAST SQUARES FIT. DAPF 350
C LET N1, N2, DENOTE THE FOLLOWING NUMBERS DAPF 360
C N1 = MAXIMAL DIMENSION FOR WHICH NO LOSS OF DAPF 370
C SIGNIFICANCE WAS INDICATED DURING FACTORIZATIONDAPF 380
C N2 = SMALLEST DIMENSION FOR WHICH THE SQUARE SUM OF DAPF 390
C THE ERRORS DOES NOT EXCEED TEST=ABS(ETA*FSQ) DAPF 400
C THEN IRES=MINO(IP,N1) IF IOP IS NONNEGATIVE DAPF 410
C AND IRES=MINO(IP,N1,N2) IF IOP IS NEGATIVE DAPF 420
C IOP - INPUT PARAMETER FOR SELECTION OF OPERATION DAPF 430
C IOP = 0 MEANS TRIANGULAR FACTORIZATION, DIVISION OF DAPF 440
C THE RIGHT HAND SIDE BY TRANSPOSE(T) AND DAPF 450
C CALCULATION OF THE SQUARE SUM OF ERRORS IS DAPF 460
C PERFORMED ONLY DAPF 470
C IOP = +1 OR -1 MEANS THE SOLUTION OF DIMENSION IRES DAPF 480
C IS CALCULATED ADDITIONALLY DAPF 490
C IOP = +2 OR -2 MEANS ALL SOLUTIONS FOR DIMENSION ONEDAPF 500
C UP TO IRES ARE CALCULATED ADDITIONALLY DAPF 510
C EPS - RELATIVE TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.DAPF 520
C A SENSIBLE VALUE IS BETWEEN 1.E-10 AND 1.E-15 DAPF 530
C ETA - RELATIVE TOLERANCE FOR TOLERATED SQUARE SUM OF DAPF 540
C ERRORS. A REALISTIC VALUE IS BETWEEN 1.E0 AND 1.E-15DAPF 550
C IER - RESULTANT ERROR PARAMETER DAPF 560
C IER =-1 MEANS NONPOSITIVE IP DAPF 570
C IER = 0 MEANS NO LOSS OF SIGNIFICANCE DETECTED DAPF 580
C AND SPECIFIED TOLERANCE OF ERRORS REACHED DAPF 590
C IER = 1 MEANS LOSS OF SIGNIFICANCE DETECTED OR DAPF 600
C SPECIFIED TOLERANCE OF ERRORS NOT REACHED DAPF 610
C DAPF 620
C REMARKS DAPF 630
C THE ABSOLUTE TOLERANCE USED INTERNALLY FOR TEST ON LOSS OF DAPF 640
C SIGNIFICANCE IS TOL=ABS(EPS*SNGL(WORK(1))). DAPF 650
C THE ABSOLUTE TOLERANCE USED INTERNALLY FOR THE SQUARE SUM OFDAPF 660
C ERRORS IS ABS(ETA*SNGL(FSQ)). DAPF 670
C IOP GREATER THAN 2 HAS THE SAME EFFECT AS IOP = 2. DAPF 680
C IOP LESS THAN -2 HAS THE SAME EFFECT AS IOP =-2. DAPF 690
C IRES = 0 MEANS THE ABSOLUTE VALUE OF EPS IS NOT LESS THAN DAPF 700
C ONE AND/OR WORK(1) IS NOT POSITIVE AND/OR IP IS NOT POSITIVEDAPF 710
C DAPF 720
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED DAPF 730
C NONE DAPF 740
C DAPF 750
C METHOD DAPF 760
C CALCULATION OF THE LEAST SQUARES FITS IS DONE USING DAPF 770
C CHOLESKYS SQUARE ROOT METHOD FOR SYMMETRIC FACTORIZATION. DAPF 780
C THE INCORPORATED TEST ON LOSS OF SIGNIFICANCE MEANS EACH DAPF 790
C RADICAND MUST BE GREATER THAN THE INTERNAL ABSOLUTE DAPF 800
C TOLERANCE TOL. DAPF 810
C IN CASE OF LOSS OF SIGNIFICANCE IN THE ABOVE SENSE ONLY A DAPF 820
C SUBSYSTEM OF THE NORMAL EQUATIONS IS SOLVED. DAPF 830
C IN CASE OF NEGATIVE IOP THE TRIANGULAR FACTORIZATION IS DAPF 840
C TERMINATED PREMATURELY EITHER IF THE SQUARE SUM OF THE DAPF 850
C ERRORS DOES NOT EXCEED ETA*FSQ OR IF THERE IS INDICATION DAPF 860
C FOR LOSS OF SIGNIFICANCE DAPF 870
C DAPF 880
C ..................................................................DAPF 890
C DAPF 900
SUBROUTINE DAPFS(WORK,IP,IRES,IOP,EPS,ETA,IER) DAPF 910
C DAPF 920
C DAPF 930
C DIMENSIONED DUMMY VARIABLES DAPF 940
DIMENSION WORK(1) DAPF 950
DOUBLE PRECISION WORK,SUM,PIV DAPF 960
IRES=0 DAPF 970
C DAPF 980
C TEST OF SPECIFIED DIMENSION DAPF 990
IF(IP)1,1,2 DAPF1000
C DAPF1010
C ERROR RETURN IN CASE OF ILLEGAL DIMENSION DAPF1020
1 IER=-1 DAPF1030
RETURN DAPF1040
C DAPF1050
C INITIALIZE FACTORIZATION PROCESS DAPF1060
2 IPIV=0 DAPF1070
IPP1=IP+1 DAPF1080
IER=1 DAPF1090
ITE=IP*IPP1/2 DAPF1100
IEND=ITE+IPP1 DAPF1110
TOL=ABS(EPS*SNGL(WORK(1))) DAPF1120
TEST=ABS(ETA*SNGL(WORK(IEND))) DAPF1130
C DAPF1140
C START LOOP OVER ALL ROWS OF WORK DAPF1150
DO 11 I=1,IP DAPF1160
IPIV=IPIV+I DAPF1170
JA=IPIV-IRES DAPF1180
JE=IPIV-1 DAPF1190
C DAPF1200
C FORM SCALAR PRODUCT NEEDED TO MODIFY CURRENT ROW ELEMENTS DAPF1210
JK=IPIV DAPF1220
DO 9 K=I,IPP1 DAPF1230
SUM=0.D0 DAPF1240
IF(IRES)5,5,3 DAPF1250
3 JK=JK-IRES DAPF1260
DO 4 J=JA,JE DAPF1270
SUM=SUM+WORK(J)*WORK(JK) DAPF1280
4 JK=JK+1 DAPF1290
5 IF(JK-IPIV)6,6,8 DAPF1300
C DAPF1310
C TEST FOR LOSS OF SIGNIFICANCE DAPF1320
6 SUM=WORK(IPIV)-SUM DAPF1330
IF(SNGL(SUM)-TOL)12,12,7 DAPF1340
7 SUM=DSQRT(SUM) DAPF1350
WORK(IPIV)=SUM DAPF1360
PIV=1.D0/SUM DAPF1370
GOTO 9 DAPF1380
C DAPF1390
C UPDATE OFF-DIAGONAL TERMS DAPF1400
8 SUM=(WORK(JK)-SUM)*PIV DAPF1410
WORK(JK)=SUM DAPF1420
9 JK=JK+K DAPF1430
C DAPF1440
C UPDATE SQUARE SUM OF ERRORS DAPF1450
WORK(IEND)=WORK(IEND)-SUM*SUM DAPF1460
C DAPF1470
C RECORD ADDRESS OF LAST PIVOT ELEMENT DAPF1480
IRES=IRES+1 DAPF1490
IADR=IPIV DAPF1500
C DAPF1510
C TEST FOR TOLERABLE ERROR IF SPECIFIED DAPF1520
IF(IOP)10,11,11 DAPF1530
10 IF(SNGL(WORK(IEND))-TEST)13,13,11 DAPF1540
11 CONTINUE DAPF1550
IF(IOP)12,22,12 DAPF1560
C DAPF1570
C PERFORM BACK SUBSTITUTION IF SPECIFIED DAPF1580
12 IF(IOP)14,23,14 DAPF1590
13 IER=0 DAPF1600
14 IPIV=IRES DAPF1610
15 IF(IPIV)23,23,16 DAPF1620
16 SUM=0.D0 DAPF1630
JA=ITE+IPIV DAPF1640
JJ=IADR DAPF1650
JK=IADR DAPF1660
K=IPIV DAPF1670
DO 19 I=1,IPIV DAPF1680
WORK(JK)=(WORK(JA)-SUM)/WORK(JJ) DAPF1690
IF(K-1)20,20,17 DAPF1700
17 JE=JJ-1 DAPF1710
SUM=0.D0 DAPF1720
DO 18 J=K,IPIV DAPF1730
SUM=SUM+WORK(JK)*WORK(JE) DAPF1740
JK=JK+1 DAPF1750
18 JE=JE+J DAPF1760
JK=JE-IPIV DAPF1770
JA=JA-1 DAPF1780
JJ=JJ-K DAPF1790
19 K=K-1 DAPF1800
20 IF(IOP/2)21,23,21 DAPF1810
21 IADR=IADR-IPIV DAPF1820
IPIV=IPIV-1 DAPF1830
GOTO 15 DAPF1840
C DAPF1850
C NORMAL RETURN DAPF1860
22 IER=0 DAPF1870
23 RETURN DAPF1880
END DAPF1890