Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-02 - 43,50145/darat.ssp
There are 2 other files named darat.ssp in the archive. Click here to see a list.
C                                                                       DARA  10
C     ..................................................................DARA  20
C                                                                       DARA  30
C        SUBROUTINE DARAT                                               DARA  40
C                                                                       DARA  50
C        PURPOSE                                                        DARA  60
C           CALCULATE BEST RATIONAL APPROXIMATION OF A DISCRETE         DARA  70
C           FUNCTION IN THE LEAST SQUARES SENSE                         DARA  80
C                                                                       DARA  90
C        USAGE                                                          DARA 100
C           CALL DARAT(DATI,N,WORK,P,IP,IQ,IER)                         DARA 110
C                                                                       DARA 120
C        DESCRIPTION OF PARAMETERS                                      DARA 130
C           DATI  - TWODIMENSIONAL ARRAY WITH 3 COLUMNS AND N ROWS      DARA 140
C                   THE FIRST COLUMN MUST CONTAIN THE GIVEN ARGUMENTS,  DARA 150
C                   THE SECOND COLUMN THE GIVEN FUNCTION VALUES AND     DARA 160
C                   THE THIRD COLUMN THE GIVEN WEIGHTS IF ANY.          DARA 170
C                   IF NO WEIGHTS ARE TO BE USED THEN THE THIRD         DARA 180
C                   COLUMN MAY BE DROPPED , EXCEPT THE FIRST ELEMENT    DARA 190
C                   WHICH MUST CONTAIN A NONPOSITIVE VALUE              DARA 200
C                   DATI MUST BE OF DOUBLE PRECISION                    DARA 210
C           N     - NUMBER OF NODES OF THE GIVEN DISCRETE FUNCTION      DARA 220
C           WORK  - WORKING STORAGE WHICH IS OF DIMENSION               DARA 230
C                   (IP+IQ)*(IP+IQ+1)+4*N+1 AT LEAST.                   DARA 240
C                   ON RETURN THE VALUES OF THE NUMERATOR ARE CONTAINED DARA 250
C                   IN WORK(N+1) UP TO WORK(2*N), WHILE THE VALUES OF   DARA 260
C                   THE DENOMINATOR ARE STORED IN WORK(2*N+1) UP TO     DARA 270
C                   WORK(3*N)                                           DARA 280
C                   WORK MUST BE OF DOUBLE PRECISION                    DARA 290
C           P     - RESULTANT COEFFICIENT VECTOR OF DENOMINATOR AND     DARA 300
C                   NUMERATOR. THE DENOMINATOR IS STORED IN FIRST IQ    DARA 310
C                   LOCATIONS, THE NUMERATOR IN THE FOLLOWING IP        DARA 320
C                   LOCATIONS.                                          DARA 330
C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH.          DARA 340
C                   P MUST BE OF DOUBLE PRECISION                       DARA 350
C           IP    - DIMENSION OF THE NUMERATOR   (INPUT VALUE)          DARA 360
C           IQ    - DIMENSION OF THE DENOMINATOR (INPUT VALUE)          DARA 370
C           IER   - RESULTANT ERROR PARAMETER                           DARA 380
C                   IER =-1 MEANS FORMAL ERRORS                         DARA 390
C                   IER = 0 MEANS NO ERRORS                             DARA 400
C                   IER = 1,2 MEANS POOR CONVERGENCE OF ITERATION       DARA 410
C                   IER IS ALSO USED AS INPUT VALUE                     DARA 420
C                   A NONZERO INPUT VALUE INDICATES AVAILABILITY OF AN  DARA 430
C                   INITIAL APPROXIMATION STORED IN P                   DARA 440
C                                                                       DARA 450
C        REMARKS                                                        DARA 460
C           THE COEFFICIENT VECTORS OF THE DENOMINATOR AND NUMERATOR    DARA 470
C           OF THE RATIONAL APPROXIMATION ARE BOTH STORED IN P          DARA 480
C           STARTING WITH LOW POWERS (DENOMINATOR FIRST).               DARA 490
C           IP+IQ MUST NOT EXCEED N, ALL THREE VALUES MUST BE POSITIVE. DARA 500
C           SINCE CHEBYSHEV POLYNOMIALS ARE USED AS FUNDAMENTAL         DARA 510
C           FUNCTIONS, THE ARGUMENTS SHOULD BE REDUCED TO THE INTERVAL  DARA 520
C           (-1,1). THIS CAN ALWAYS BE ACCOMPLISHED BY MEANS OF A LINEARDARA 530
C           TRANSFORMATION OF THE ORIGINALLY GIVEN ARGUMENTS.           DARA 540
C           IF A FIT IN OTHER FUNCTIONS IS REQUIRED, DCNP AND DCNPS MUSTDARA 550
C           BE REPLACED BY SUBROUTINES WHICH ARE OF ANALOGOUS DESIGN.   DARA 560
C                                                                       DARA 570
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DARA 580
C           DAPLL, DAPFS, DFRAT, DCNPS, DCNP                            DARA 590
C           DCNP IS REQUIRED WITHIN DFRAT                               DARA 600
C                                                                       DARA 610
C        METHOD                                                         DARA 620
C           THE ITERATIVE SCHEME USED FOR CALCULATION OF THE            DARA 630
C           APPROXIMATION IS REPEATED SOLUTION OF THE NORMAL EQUATIONS  DARA 640
C           WHICH ARE OBTAINED BY LINEARIZATION.                        DARA 650
C           A REFINED TECHNIQUE OF THIS LINEAR LEAST SQUARES APPROACH   DARA 660
C           IS USED WHICH GUARANTEES THAT THE DENOMINATOR IS FREE OF    DARA 670
C           ZEROES WITHIN THE APPROXIMATION INTERVAL.                   DARA 680
C           FOR REFERENCE SEE                                           DARA 690
C           D.BRAESS, UEBER DAEMPFUNG BEI MINIMALISIERUNGSVERFAHREN,    DARA 700
C           COMPUTING(1966), VOL.1, ED.3, PP.264-272.                   DARA 710
C           D.W.MARQUARDT, AN ALGORITHM FOR LEAST-SQUARES ESTIMATION    DARA 720
C           OF NONLINEAR PARAMETERS,                                    DARA 730
C           JSIAM(1963), VOL.11, ED.2, PP.431-441.                      DARA 740
C                                                                       DARA 750
C     ..................................................................DARA 760
C                                                                       DARA 770
      SUBROUTINE DARAT(DATI,N,WORK,P,IP,IQ,IER)                         DARA 780
C                                                                       DARA 790
C                                                                       DARA 800
      EXTERNAL DFRAT                                                    DARA 810
C                                                                       DARA 820
C        DIMENSIONED LOCAL VARIABLE                                     DARA 830
      DIMENSION IERV(3)                                                 DARA 840
C                                                                       DARA 850
C        DIMENSIONED DUMMY VARIABLES                                    DARA 860
      DIMENSION DATI(1),WORK(1),P(1)                                    DARA 870
      DOUBLE PRECISION DATI,WORK,P,T,OSUM,DIAG,RELAX,SUM,SSOE,SAVE      DARA 880
C                                                                       DARA 890
C        INITIALIZE TESTVALUES                                          DARA 900
      LIMIT=20                                                          DARA 910
      ETA=1.E-29                                                        DARA 920
      EPS=1.E-14                                                        DARA 930
C                                                                       DARA 940
C        CHECK FOR FORMAL ERRORS                                        DARA 950
      IF(N)4,4,1                                                        DARA 960
    1 IF(IP)4,4,2                                                       DARA 970
    2 IF(IQ)4,4,3                                                       DARA 980
    3 IPQ=IP+IQ                                                         DARA 990
      IF(N-IPQ)4,5,5                                                    DARA1000
C                                                                       DARA1010
C        ERROR RETURN IN CASE OF FORMAL ERRORS                          DARA1020
    4 IER=-1                                                            DARA1030
      RETURN                                                            DARA1040
C                                                                       DARA1050
C        INITIALIZE ITERATION PROCESS                                   DARA1060
    5 KOUNT=0                                                           DARA1070
      IERV(2)=IP                                                        DARA1080
      IERV(3)=IQ                                                        DARA1090
      NDP=N+N+1                                                         DARA1100
      NNE=NDP+NDP                                                       DARA1110
      IX=IPQ-1                                                          DARA1120
      IQP1=IQ+1                                                         DARA1130
      IRHS=NNE+IPQ*IX/2                                                 DARA1140
      IEND=IRHS+IX                                                      DARA1150
C                                                                       DARA1160
C        TEST FOR AVAILABILITY OF AN INITIAL APPROXIMATION              DARA1170
      IF(IER)8,6,8                                                      DARA1180
C                                                                       DARA1190
C        INITIALIZE NUMERATOR AND DENOMINATOR                           DARA1200
    6 DO 7 I=2,IPQ                                                      DARA1210
    7 P(I)=0.D0                                                         DARA1220
      P(1)=1.D0                                                         DARA1230
C                                                                       DARA1240
C        CALCULATE VALUES OF NUMERATOR AND DENOMINATOR FOR INITIAL      DARA1250
C        APPROXIMATION                                                  DARA1260
    8 DO 9 J=1,N                                                        DARA1270
      T=DATI(J)                                                         DARA1280
      I=J+N                                                             DARA1290
      CALL DCNPS(WORK(I),T,P(IQP1),IP)                                  DARA1300
      K=I+N                                                             DARA1310
    9 CALL DCNPS(WORK(K),T,P,IQ)                                        DARA1320
C                                                                       DARA1330
C        SET UP NORMAL EQUATIONS (MAIN LOOP OF ITERATION)               DARA1340
   10 CALL DAPLL(DFRAT,N,IX,WORK,WORK(IEND+1),DATI,IERV)                DARA1350
C                                                                       DARA1360
C        CHECK FOR ZERO DENOMINATOR                                     DARA1370
      IF(IERV(1))4,11,4                                                 DARA1380
   11 INCR=0                                                            DARA1390
      RELAX=2.D0                                                        DARA1400
C                                                                       DARA1410
C        RESTORE MATRIX IN WORKING STORAGE                              DARA1420
   12 J=IEND                                                            DARA1430
      DO 13 I=NNE,IEND                                                  DARA1440
      J=J+1                                                             DARA1450
   13 WORK(I)=WORK(J)                                                   DARA1460
      IF(KOUNT)14,14,15                                                 DARA1470
C                                                                       DARA1480
C        SAVE SQUARE SUM OF ERRORS                                      DARA1490
   14 OSUM=WORK(IEND)                                                   DARA1500
      DIAG=OSUM*EPS                                                     DARA1510
      K=IQ                                                              DARA1520
C                                                                       DARA1530
C        ADD CONSTANT TO DIAGONAL                                       DARA1540
      IF(WORK(NNE))17,17,19                                             DARA1550
   15 IF(INCR)19,19,16                                                  DARA1560
   16 K=IPQ                                                             DARA1570
   17 J=NNE-1                                                           DARA1580
      DO 18 I=1,K                                                       DARA1590
      WORK(J)=WORK(J)+DIAG                                              DARA1600
   18 J=J+I                                                             DARA1610
C                                                                       DARA1620
C        SOLVE NORMAL EQUATIONS                                         DARA1630
   19 CALL DAPFS(WORK(NNE),IX,IRES,1,EPS,ETA,IER)                       DARA1640
C                                                                       DARA1650
C        CHECK FOR FAILURE OF EQUATION SOLVER                           DARA1660
      IF(IRES)4,4,20                                                    DARA1670
C                                                                       DARA1680
C        TEST FOR DEFECTIVE NORMALEQUATIONS                             DARA1690
   20 IF(IRES-IX)21,24,24                                               DARA1700
   21 IF(INCR)22,22,23                                                  DARA1710
   22 DIAG=DIAG*0.125D0                                                 DARA1720
   23 DIAG=DIAG+DIAG                                                    DARA1730
      INCR=INCR+1                                                       DARA1740
C                                                                       DARA1750
C        START WITH OVER RELAXATION                                     DARA1760
      RELAX=8.D0                                                        DARA1770
      IF(INCR-LIMIT)12,45,45                                            DARA1780
C                                                                       DARA1790
C        CALCULATE VALUES OF CHANGE OF NUMERATOR AND DENOMINATOR        DARA1800
   24 L=NDP                                                             DARA1810
      J=NNE+IRES*(IRES-1)/2-1                                           DARA1820
      K=J+IQ                                                            DARA1830
      WORK(J)=0.D0                                                      DARA1840
      IRQ=IQ                                                            DARA1850
      IRP=IRES-IQ+1                                                     DARA1860
      IF(IRP)25,26,26                                                   DARA1870
   25 IRQ=IRES+1                                                        DARA1880
   26 DO 29 I=1,N                                                       DARA1890
      T=DATI(I)                                                         DARA1900
      WORK(I)=0.D0                                                      DARA1910
      CALL DCNPS(WORK(I),T,WORK(K),IRP)                                 DARA1920
      M=L+N                                                             DARA1930
      CALL DCNPS(WORK(M),T,WORK(J),IRQ)                                 DARA1940
      IF(WORK(M)*WORK(L))27,29,29                                       DARA1950
   27 SUM=WORK(L)/WORK(M)                                               DARA1960
      IF(RELAX+SUM)29,29,28                                             DARA1970
   28 RELAX=-SUM                                                        DARA1980
   29 L=L+1                                                             DARA1990
C                                                                       DARA2000
C        MODIFY RELAXATION FACTOR IF NECESSARY                          DARA2010
      SSOE=OSUM                                                         DARA2020
      ITER=LIMIT                                                        DARA2030
   30 SUM=0.D0                                                          DARA2040
      RELAX=RELAX*0.5D0                                                 DARA2050
      DO 32 I=1,N                                                       DARA2060
      M=I+N                                                             DARA2070
      K=M+N                                                             DARA2080
      L=K+N                                                             DARA2090
      SAVE=DATI(M)-(WORK(M)+RELAX*WORK(I))/(WORK(K)+RELAX*WORK(L))      DARA2100
      SAVE=SAVE*SAVE                                                    DARA2110
      IF(DATI(NDP))32,32,31                                             DARA2120
   31 SAVE=SAVE*DATI(K)                                                 DARA2130
   32 SUM=SUM+SAVE                                                      DARA2140
      IF(ITER)45,33,33                                                  DARA2150
   33 ITER=ITER-1                                                       DARA2160
      IF(SUM-OSUM)34,37,35                                              DARA2170
   34 OSUM=SUM                                                          DARA2180
      GOTO 30                                                           DARA2190
C                                                                       DARA2200
C        TEST FOR IMPROVEMENT                                           DARA2210
   35 IF(OSUM-SSOE)36,30,30                                             DARA2220
   36 RELAX=RELAX+RELAX                                                 DARA2230
   37 T=0.                                                              DARA2240
      SAVE=0.D0                                                         DARA2250
      K=IRES+1                                                          DARA2260
      DO 38 I=2,K                                                       DARA2270
      J=J+1                                                             DARA2280
      T=T+DABS(P(I))                                                    DARA2290
      P(I)=P(I)+RELAX*WORK(J)                                           DARA2300
   38 SAVE=SAVE+DABS(P(I))                                              DARA2310
C                                                                       DARA2320
C        UPDATE CURRENT VALUES OF NUMERATOR AND DENOMINATOR             DARA2330
      DO 39 I=1,N                                                       DARA2340
      J=I+N                                                             DARA2350
      K=J+N                                                             DARA2360
      L=K+N                                                             DARA2370
      WORK(J)=WORK(J)+RELAX*WORK(I)                                     DARA2380
   39 WORK(K)=WORK(K)+RELAX*WORK(L)                                     DARA2390
C                                                                       DARA2400
C        TEST FOR CONVERGENCE                                           DARA2410
      IF(INCR)40,40,42                                                  DARA2420
   40 IF(SSOE-OSUM-RELAX*OSUM*DBLE(EPS))46,46,41                        DARA2430
   41 IF(DABS(T-SAVE)-RELAX*SAVE*DBLE(EPS))46,46,42                     DARA2440
   42 IF(OSUM-SAVE*DBLE(ETA))46,46,43                                   DARA2450
   43 KOUNT=KOUNT+1                                                     DARA2460
      IF(KOUNT-LIMIT)10,44,44                                           DARA2470
C                                                                       DARA2480
C        ERROR RETURN IN CASE OF POOR CONVERGENCE                       DARA2490
   44 IER=2                                                             DARA2500
      RETURN                                                            DARA2510
   45 IER=1                                                             DARA2520
      RETURN                                                            DARA2530
C                                                                       DARA2540
C        NORMAL RETURN                                                  DARA2550
   46 IER=0                                                             DARA2560
      RETURN                                                            DARA2570
      END                                                               DARA2580