Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap2_198111 - decus/20-0026/qatr.ssp
There are 2 other files named qatr.ssp in the archive. Click here to see a list.
C                                                                       QATR  10
C     ..................................................................QATR  20
C                                                                       QATR  30
C        SUBROUTINE QATR                                                QATR  40
C                                                                       QATR  50
C        PURPOSE                                                        QATR  60
C           TO COMPUTE AN APPROXIMATION FOR INTEGRAL(FCT(X), SUMMED     QATR  70
C           OVER X FROM XL TO XU).                                      QATR  80
C                                                                       QATR  90
C        USAGE                                                          QATR 100
C           CALL QATR (XL,XU,EPS,NDIM,FCT,Y,IER,AUX)                    QATR 110
C           PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT.               QATR 120
C                                                                       QATR 130
C        DESCRIPTION OF PARAMETERS                                      QATR 140
C           XL     - THE LOWER BOUND OF THE INTERVAL.                   QATR 150
C           XU     - THE UPPER BOUND OF THE INTERVAL.                   QATR 160
C           EPS    - THE UPPER BOUND OF THE ABSOLUTE ERROR.             QATR 170
C           NDIM   - THE DIMENSION OF THE AUXILIARY STORAGE ARRAY AUX.  QATR 180
C                    NDIM-1 IS THE MAXIMAL NUMBER OF BISECTIONS OF      QATR 190
C                    THE INTERVAL (XL,XU).                              QATR 200
C           FCT    - THE NAME OF THE EXTERNAL FUNCTION SUBPROGRAM USED. QATR 210
C           Y      - THE RESULTING APPROXIMATION FOR THE INTEGRAL VALUE.QATR 220
C           IER    - A RESULTING ERROR PARAMETER.                       QATR 230
C           AUX    - AN AUXILIARY STORAGE ARRAY WITH DIMENSION NDIM.    QATR 240
C                                                                       QATR 250
C        REMARKS                                                        QATR 260
C           ERROR PARAMETER IER IS CODED IN THE FOLLOWING FORM          QATR 270
C           IER=0  - IT WAS POSSIBLE TO REACH THE REQUIRED ACCURACY.    QATR 280
C                    NO ERROR.                                          QATR 290
C           IER=1  - IT IS IMPOSSIBLE TO REACH THE REQUIRED ACCURACY    QATR 300
C                    BECAUSE OF ROUNDING ERRORS.                        QATR 310
C           IER=2  - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE NDIM   QATR 320
C                    IS LESS THAN 5, OR THE REQUIRED ACCURACY COULD NOT QATR 330
C                    BE REACHED WITHIN NDIM-1 STEPS. NDIM SHOULD BE     QATR 340
C                    INCREASED.                                         QATR 350
C                                                                       QATR 360
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  QATR 370
C           THE EXTERNAL FUNCTION SUBPROGRAM FCT(X) MUST BE CODED BY    QATR 380
C           THE USER. ITS ARGUMENT X SHOULD NOT BE DESTROYED.           QATR 390
C                                                                       QATR 400
C        METHOD                                                         QATR 410
C           EVALUATION OF Y IS DONE BY MEANS OF TRAPEZOIDAL RULE IN     QATR 420
C           CONNECTION WITH ROMBERGS PRINCIPLE. ON RETURN Y CONTAINS    QATR 430
C           THE BEST POSSIBLE APPROXIMATION OF THE INTEGRAL VALUE AND   QATR 440
C           VECTOR AUX THE UPWARD DIAGONAL OF ROMBERG SCHEME.           QATR 450
C           COMPONENTS AUX(I) (I=1,2,...,IEND, WITH IEND LESS THAN OR   QATR 460
C           EQUAL TO NDIM) BECOME APPROXIMATIONS TO INTEGRAL VALUE WITH QATR 470
C           DECREASING ACCURACY BY MULTIPLICATION WITH (XU-XL).         QATR 480
C           FOR REFERENCE, SEE                                          QATR 490
C           (1) FILIPPI, DAS VERFAHREN VON ROMBERG-STIEFEL-BAUER ALS    QATR 500
C               SPEZIALFALL DES ALLGEMEINEN PRINZIPS VON RICHARDSON,    QATR 510
C               MATHEMATIK-TECHNIK-WIRTSCHAFT, VOL.11, ISS.2 (1964),    QATR 520
C               PP.49-54.                                               QATR 530
C           (2) BAUER, ALGORITHM 60, CACM, VOL.4, ISS.6 (1961), PP.255. QATR 540
C                                                                       QATR 550
C     ..................................................................QATR 560
C                                                                       QATR 570
      SUBROUTINE QATR(XL,XU,EPS,NDIM,FCT,Y,IER,AUX)                     QATR 580
C                                                                       QATR 590
C                                                                       QATR 600
      DIMENSION AUX(1)                                                  QATR 610
C                                                                       QATR 620
C     PREPARATIONS OF ROMBERG-LOOP                                      QATR 630
      AUX(1)=.5*(FCT(XL)+FCT(XU))                                       QATR 640
      H=XU-XL                                                           QATR 650
      IF(NDIM-1)8,8,1                                                   QATR 660
    1 IF(H)2,10,2                                                       QATR 670
C                                                                       QATR 680
C     NDIM IS GREATER THAN 1 AND H IS NOT EQUAL TO 0.                   QATR 690
    2 HH=H                                                              QATR 700
      E=EPS/ABS(H)                                                      QATR 710
      DELT2=0.                                                          QATR 720
      P=1.                                                              QATR 730
      JJ=1                                                              QATR 740
      DO 7 I=2,NDIM                                                     QATR 750
      Y=AUX(1)                                                          QATR 760
      DELT1=DELT2                                                       QATR 770
      HD=HH                                                             QATR 780
      HH=.5*HH                                                          QATR 790
      P=.5*P                                                            QATR 800
      X=XL+HH                                                           QATR 810
      SM=0.                                                             QATR 820
      DO 3 J=1,JJ                                                       QATR 830
      SM=SM+FCT(X)                                                      QATR 840
    3 X=X+HD                                                            QATR 850
      AUX(I)=.5*AUX(I-1)+P*SM                                           QATR 860
C     A NEW APPROXIMATION OF INTEGRAL VALUE IS COMPUTED BY MEANS OF     QATR 870
C     TRAPEZOIDAL RULE.                                                 QATR 880
C                                                                       QATR 890
C     START OF ROMBERGS EXTRAPOLATION METHOD.                           QATR 900
      Q=1.                                                              QATR 910
      JI=I-1                                                            QATR 920
      DO 4 J=1,JI                                                       QATR 930
      II=I-J                                                            QATR 940
      Q=Q+Q                                                             QATR 950
      Q=Q+Q                                                             QATR 960
    4 AUX(II)=AUX(II+1)+(AUX(II+1)-AUX(II))/(Q-1.)                      QATR 970
C     END OF ROMBERG-STEP                                               QATR 980
C                                                                       QATR 990
      DELT2=ABS(Y-AUX(1))                                               QATR1000
      IF(I-5)7,5,5                                                      QATR1010
    5 IF(DELT2-E)10,10,6                                                QATR1020
    6 IF(DELT2-DELT1)7,11,11                                            QATR1030
    7 JJ=JJ+JJ                                                          QATR1040
    8 IER=2                                                             QATR1050
    9 Y=H*AUX(1)                                                        QATR1060
      RETURN                                                            QATR1070
   10 IER=0                                                             QATR1080
      GO TO 9                                                           QATR1090
   11 IER=1                                                             QATR1100
      Y=H*Y                                                             QATR1110
      RETURN                                                            QATR1120
      END                                                               QATR1130