Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/dapfs.ssp
There are 2 other files named dapfs.ssp in the archive. Click here to see a list.
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