Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/arat.ssp
There are 2 other files named arat.ssp in the archive. Click here to see a list.
C                                                                       ARAT  10
C     ..................................................................ARAT  20
C                                                                       ARAT  30
C        SUBROUTINE ARAT                                                ARAT  40
C                                                                       ARAT  50
C        PURPOSE                                                        ARAT  60
C           CALCULATE BEST RATIONAL APPROXIMATION OF A DISCRETE         ARAT  70
C           FUNCTION IN THE LEAST SQUARES SENSE                         ARAT  80
C                                                                       ARAT  90
C        USAGE                                                          ARAT 100
C           CALL ARAT(DATI,N,WORK,P,IP,IQ,IER)                          ARAT 110
C                                                                       ARAT 120
C        DESCRIPTION OF PARAMETERS                                      ARAT 130
C           DATI  - TWODIMENSIONAL ARRAY WITH 3 COLUMNS AND N ROWS      ARAT 140
C                   THE FIRST COLUMN MUST CONTAIN THE GIVEN ARGUMENTS,  ARAT 150
C                   THE SECOND COLUMN THE GIVEN FUNCTION VALUES AND     ARAT 160
C                   THE THIRD COLUMN THE GIVEN WEIGHTS IF ANY.          ARAT 170
C                   IF NO WEIGHTS ARE TO BE USED THEN THE THIRD         ARAT 180
C                   COLUMN MAY BE DROPPED , EXCEPT THE FIRST ELEMENT    ARAT 190
C                   WHICH MUST CONTAIN A NONPOSITIVE VALUE              ARAT 200
C           N     - NUMBER OF NODES OF THE GIVEN DISCRETE FUNCTION      ARAT 210
C           WORK  - WORKING STORAGE WHICH IS OF DIMENSION               ARAT 220
C                   (IP+IQ)*(IP+IQ+1)+4*N+1 AT LEAST.                   ARAT 230
C                   ON RETURN THE VALUES OF THE NUMERATOR ARE CONTAINED ARAT 240
C                   IN WORK(N+1) UP TO WORK(2*N), WHILE THE VALUES OF   ARAT 250
C                   THE DENOMINATOR ARE STORED IN WORK(2*N+1) UP TO     ARAT 260
C                   WORK(3*N)                                           ARAT 270
C           P     - RESULTANT COEFFICIENT VECTOR OF DENOMINATOR AND     ARAT 280
C                   NUMERATOR. THE DENOMINATOR IS STORED IN FIRST IQ    ARAT 290
C                   LOCATIONS, THE NUMERATOR IN THE FOLLOWING IP        ARAT 300
C                   LOCATIONS.                                          ARAT 310
C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH.          ARAT 320
C           IP    - DIMENSION OF THE NUMERATOR   (INPUT VALUE)          ARAT 330
C           IQ    - DIMENSION OF THE DENOMINATOR (INPUT VALUE)          ARAT 340
C           IER   - RESULTANT ERROR PARAMETER                           ARAT 350
C                   IER =-1 MEANS FORMAL ERRORS                         ARAT 360
C                   IER = 0 MEANS NO ERRORS                             ARAT 370
C                   IER = 1,2 MEANS POOR CONVERGENCE OF ITERATION       ARAT 380
C                   IER IS ALSO USED AS INPUT VALUE                     ARAT 390
C                   A NONZERO INPUT VALUE INDICATES AVAILABILITY OF AN  ARAT 400
C                   INITIAL APPROXIMATION STORED IN P                   ARAT 410
C                                                                       ARAT 420
C        REMARKS                                                        ARAT 430
C           THE COEFFICIENT VECTORS OF THE DENOMINATOR AND NUMERATOR    ARAT 440
C           OF THE RATIONAL APPROXIMATION ARE BOTH STORED IN P          ARAT 450
C           STARTING WITH LOW POWERS (DENOMINATOR FIRST).               ARAT 460
C           IP+IQ MUST NOT EXCEED N, ALL THREE VALUES MUST BE POSITIVE. ARAT 470
C           SINCE CHEBYSHEV POLYNOMIALS ARE USED AS FUNDAMENTAL         ARAT 480
C           FUNCTIONS, THE ARGUMENTS SHOULD BE REDUCED TO THE INTERVAL  ARAT 490
C           (-1,1). THIS CAN ALWAYS BE ACCOMPLISHED BY MEANS OF A LINEARARAT 500
C           TRANSFORMATION OF THE ORIGINALLY GIVEN ARGUMENTS.           ARAT 510
C           IF A FIT IN OTHER FUNCTIONS IS REQUIRED, CNP AND CNPS MUST  ARAT 520
C           BE REPLACED BY SUBROUTINES WHICH ARE OF ANALOGOUS DESIGN.   ARAT 530
C                                                                       ARAT 540
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  ARAT 550
C           APLL, APFS, FRAT, CNPS, CNP                                 ARAT 560
C           CNP IS REQUIRED WITHIN FRAT                                 ARAT 570
C                                                                       ARAT 580
C        METHOD                                                         ARAT 590
C           THE ITERATIVE SCHEME USED FOR CALCULATION OF THE            ARAT 600
C           APPROXIMATION IS REPEATED SOLUTION OF THE NORMAL EQUATIONS  ARAT 610
C           WHICH ARE OBTAINED BY LINEARIZATION.                        ARAT 620
C           A REFINED TECHNIQUE OF THIS LINEAR LEAST SQUARES APPROACH   ARAT 630
C           IS USED WHICH GUARANTEES THAT THE DENOMINATOR IS FREE OF    ARAT 640
C           ZEROES WITHIN THE APPROXIMATION INTERVAL.                   ARAT 650
C           FOR REFERENCE SEE                                           ARAT 660
C           D.BRAESS, UEBER DAEMPFUNG BEI MINIMALISIERUNGSVERFAHREN,    ARAT 670
C           COMPUTING(1966), VOL.1, ED.3, PP.264-272.                   ARAT 680
C           D.W.MARQUARDT, AN ALGORITHM FOR LEAST-SQUARES ESTIMATION    ARAT 690
C           OF NONLINEAR PARAMETERS,                                    ARAT 700
C           JSIAM(1963), VOL.11, ED.2, PP.431-441.                      ARAT 710
C                                                                       ARAT 720
C     ..................................................................ARAT 730
C                                                                       ARAT 740
      SUBROUTINE ARAT(DATI,N,WORK,P,IP,IQ,IER)                          ARAT 750
C                                                                       ARAT 760
C                                                                       ARAT 770
      EXTERNAL FRAT                                                     ARAT 780
C                                                                       ARAT 790
C        DIMENSIONED LOCAL VARIABLE                                     ARAT 800
      DIMENSION IERV(3)                                                 ARAT 810
C                                                                       ARAT 820
C        DIMENSIONED DUMMY VARIABLES                                    ARAT 830
      DIMENSION DATI(1),WORK(1),P(1)                                    ARAT 840
C                                                                       ARAT 850
C        INITIALIZE TESTVALUES                                          ARAT 860
      LIMIT=20                                                          ARAT 870
      ETA =1.E-11                                                       ARAT 880
      EPS=1.E-5                                                         ARAT 890
C                                                                       ARAT 900
C        CHECK FOR FORMAL ERRORS                                        ARAT 910
      IF(N)4,4,1                                                        ARAT 920
    1 IF(IP)4,4,2                                                       ARAT 930
    2 IF(IQ)4,4,3                                                       ARAT 940
    3 IPQ=IP+IQ                                                         ARAT 950
      IF(N-IPQ)4,5,5                                                    ARAT 960
C                                                                       ARAT 970
C        ERROR RETURN IN CASE OF FORMAL ERRORS                          ARAT 980
    4 IER=-1                                                            ARAT 990
      RETURN                                                            ARAT1000
C                                                                       ARAT1010
C        INITIALIZE ITERATION PROCESS                                   ARAT1020
    5 KOUNT=0                                                           ARAT1030
      IERV(2)=IP                                                        ARAT1040
      IERV(3)=IQ                                                        ARAT1050
      NDP=N+N+1                                                         ARAT1060
      NNE=NDP+NDP                                                       ARAT1070
      IX=IPQ-1                                                          ARAT1080
      IQP1=IQ+1                                                         ARAT1090
      IRHS=NNE+IPQ*IX/2                                                 ARAT1100
      IEND=IRHS+IX                                                      ARAT1110
C                                                                       ARAT1120
C        TEST FOR AVAILABILITY OF AN INITIAL APPROXIMATION              ARAT1130
      IF(IER)8,6,8                                                      ARAT1140
C                                                                       ARAT1150
C        INITIALIZE NUMERATOR AND DENOMINATOR                           ARAT1160
    6 DO 7 I=2,IPQ                                                      ARAT1170
    7 P(I)=0.                                                           ARAT1180
      P(1)=1.                                                           ARAT1190
C                                                                       ARAT1200
C        CALCULATE VALUES OF NUMERATOR AND DENOMINATOR FOR INITIAL      ARAT1210
C        APPROXIMATION                                                  ARAT1220
    8 DO 9 J=1,N                                                        ARAT1230
      T=DATI(J)                                                         ARAT1240
      I=J+N                                                             ARAT1250
      CALL CNPS(WORK(I),T,P(IQP1),IP)                                   ARAT1260
      K=I+N                                                             ARAT1270
    9 CALL CNPS(WORK(K),T,P,IQ)                                         ARAT1280
C                                                                       ARAT1290
C        SET UP NORMAL EQUATIONS (MAIN LOOP OF ITERATION)               ARAT1300
   10 CALL APLL(FRAT,N,IX,WORK,WORK(IEND+1),DATI,IERV)                  ARAT1310
C                                                                       ARAT1320
C        CHECK FOR ZERO DENOMINATOR                                     ARAT1330
      IF(IERV(1))4,11,4                                                 ARAT1340
   11 INCR=0                                                            ARAT1350
      RELAX=2.                                                          ARAT1360
C                                                                       ARAT1370
C        RESTORE MATRIX IN WORKING STORAGE                              ARAT1380
   12 J=IEND                                                            ARAT1390
      DO 13 I=NNE,IEND                                                  ARAT1400
      J=J+1                                                             ARAT1410
   13 WORK(I)=WORK(J)                                                   ARAT1420
      IF(KOUNT)14,14,15                                                 ARAT1430
C                                                                       ARAT1440
C        SAVE SQUARE SUM OF ERRORS                                      ARAT1450
   14 OSUM=WORK(IEND)                                                   ARAT1460
      DIAG=OSUM*EPS                                                     ARAT1470
      K=IQ                                                              ARAT1480
C                                                                       ARAT1490
C        ADD CONSTANT TO DIAGONAL                                       ARAT1500
      IF(WORK(NNE))17,17,19                                             ARAT1510
   15 IF(INCR)19,19,16                                                  ARAT1520
   16 K=IPQ                                                             ARAT1530
   17 J=NNE-1                                                           ARAT1540
      DO 18 I=1,K                                                       ARAT1550
      WORK(J)=WORK(J)+DIAG                                              ARAT1560
   18 J=J+I                                                             ARAT1570
C                                                                       ARAT1580
C        SOLVE NORMAL EQUATIONS                                         ARAT1590
   19 CALL APFS(WORK(NNE),IX,IRES,1,EPS,ETA,IER)                        ARAT1600
C                                                                       ARAT1610
C        CHECK FOR FAILURE OF EQUATION SOLVER                           ARAT1620
      IF(IRES)4,4,20                                                    ARAT1630
C                                                                       ARAT1640
C        TEST FOR DEFECTIVE NORMALEQUATIONS                             ARAT1650
   20 IF(IRES-IX)21,24,24                                               ARAT1660
   21 IF(INCR)22,22,23                                                  ARAT1670
   22 DIAG=DIAG*0.125                                                   ARAT1680
   23 DIAG=DIAG+DIAG                                                    ARAT1690
      INCR=INCR+1                                                       ARAT1700
C                                                                       ARAT1710
C        START WITH OVER RELAXATION                                     ARAT1720
      RELAX=8.                                                          ARAT1730
      IF(INCR-LIMIT)12,45,45                                            ARAT1740
C                                                                       ARAT1750
C        CALCULATE VALUES OF CHANGE OF NUMERATOR AND DENOMINATOR        ARAT1760
   24 L=NDP                                                             ARAT1770
      J=NNE+IRES*(IRES-1)/2-1                                           ARAT1780
      K=J+IQ                                                            ARAT1790
      WORK(J)=0.                                                        ARAT1800
      IRQ=IQ                                                            ARAT1810
      IRP=IRES-IQ+1                                                     ARAT1820
      IF(IRP)25,26,26                                                   ARAT1830
   25 IRQ=IRES+1                                                        ARAT1840
   26 DO 29 I=1,N                                                       ARAT1850
      T=DATI(I)                                                         ARAT1860
      WORK(I)=0.                                                        ARAT1870
      CALL CNPS(WORK(I),T,WORK(K),IRP)                                  ARAT1880
      M=L+N                                                             ARAT1890
      CALL CNPS(WORK(M),T,WORK(J),IRQ)                                  ARAT1900
      IF(WORK(M)*WORK(L))27,29,29                                       ARAT1910
   27 SUM=WORK(L)/WORK(M)                                               ARAT1920
      IF(RELAX+SUM)29,29,28                                             ARAT1930
   28 RELAX=-SUM                                                        ARAT1940
   29 L=L+1                                                             ARAT1950
C                                                                       ARAT1960
C        MODIFY RELAXATION FACTOR IF NECESSARY                          ARAT1970
      SSOE=OSUM                                                         ARAT1980
      ITER=LIMIT                                                        ARAT1990
   30 SUM=0.                                                            ARAT2000
      RELAX=RELAX*0.5                                                   ARAT2010
      DO 32 I=1,N                                                       ARAT2020
      M=I+N                                                             ARAT2030
      K=M+N                                                             ARAT2040
      L=K+N                                                             ARAT2050
      SAVE=DATI(M)-(WORK(M)+RELAX*WORK(I))/(WORK(K)+RELAX*WORK(L))      ARAT2060
      SAVE=SAVE*SAVE                                                    ARAT2070
      IF(DATI(NDP))32,32,31                                             ARAT2080
   31 SAVE=SAVE*DATI(K)                                                 ARAT2090
   32 SUM=SUM+SAVE                                                      ARAT2100
      IF(ITER)45,33,33                                                  ARAT2110
   33 ITER=ITER-1                                                       ARAT2120
      IF(SUM-OSUM)34,37,35                                              ARAT2130
   34 OSUM=SUM                                                          ARAT2140
      GOTO 30                                                           ARAT2150
C                                                                       ARAT2160
C        TEST FOR IMPROVEMENT                                           ARAT2170
   35 IF(OSUM-SSOE)36,30,30                                             ARAT2180
   36 RELAX=RELAX+RELAX                                                 ARAT2190
   37 T=0.                                                              ARAT2200
      SAVE=0.                                                           ARAT2210
      K=IRES+1                                                          ARAT2220
      DO 38 I=2,K                                                       ARAT2230
      J=J+1                                                             ARAT2240
      T=T+ABS(P(I))                                                     ARAT2250
      P(I)=P(I)+RELAX*WORK(J)                                           ARAT2260
   38 SAVE=SAVE+ABS(P(I))                                               ARAT2270
C                                                                       ARAT2280
C        UPDATE CURRENT VALUES OF NUMERATOR AND DENOMINATOR             ARAT2290
      DO 39 I=1,N                                                       ARAT2300
      J=I+N                                                             ARAT2310
      K=J+N                                                             ARAT2320
      L=K+N                                                             ARAT2330
      WORK(J)=WORK(J)+RELAX*WORK(I)                                     ARAT2340
   39 WORK(K)=WORK(K)+RELAX*WORK(L)                                     ARAT2350
C                                                                       ARAT2360
C        TEST FOR CONVERGENCE                                           ARAT2370
      IF(INCR)40,40,42                                                  ARAT2380
   40 IF(SSOE-OSUM-RELAX*EPS*OSUM)46,46,41                              ARAT2390
   41 IF(ABS(T-SAVE)-RELAX*EPS*SAVE)46,46,42                            ARAT2400
   42 IF(OSUM-ETA*SAVE)46,46,43                                         ARAT2410
   43 KOUNT=KOUNT+1                                                     ARAT2420
      IF(KOUNT-LIMIT)10,44,44                                           ARAT2430
C                                                                       ARAT2440
C        ERROR RETURN IN CASE OF POOR CONVERGENCE                       ARAT2450
   44 IER=2                                                             ARAT2460
      RETURN                                                            ARAT2470
   45 IER=1                                                             ARAT2480
      RETURN                                                            ARAT2490
C                                                                       ARAT2500
C        NORMAL RETURN                                                  ARAT2510
   46 IER=0                                                             ARAT2520
      RETURN                                                            ARAT2530
      END                                                               ARAT2540