Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/dqsf.ssp
There are 2 other files named dqsf.ssp in the archive. Click here to see a list.
C                                                                       DQSF  10
C     ..................................................................DQSF  20
C                                                                       DQSF  30
C        SUBROUTINE DQSF                                                DQSF  40
C                                                                       DQSF  50
C        PURPOSE                                                        DQSF  60
C           TO COMPUTE THE VECTOR OF INTEGRAL VALUES FOR A GIVEN        DQSF  70
C           EQUIDISTANT TABLE OF FUNCTION VALUES.                       DQSF  80
C                                                                       DQSF  90
C        USAGE                                                          DQSF 100
C           CALL DQSF (H,Y,Z,NDIM)                                      DQSF 110
C                                                                       DQSF 120
C        DESCRIPTION OF PARAMETERS                                      DQSF 130
C           H      - DOUBLE PRECISION INCREMENT OF ARGUMENT VALUES.     DQSF 140
C           Y      - DOUBLE PRECISION INPUT VECTOR OF FUNCTION VALUES.  DQSF 150
C           Z      - RESULTING DOUBLE PRECISION VECTOR OF INTEGRAL      DQSF 160
C                    VALUES. Z MAY BE IDENTICAL WITH Y.                 DQSF 170
C           NDIM   - THE DIMENSION OF VECTORS Y AND Z.                  DQSF 180
C                                                                       DQSF 190
C        REMARKS                                                        DQSF 200
C           NO ACTION IN CASE NDIM LESS THAN 3.                         DQSF 210
C                                                                       DQSF 220
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DQSF 230
C           NONE                                                        DQSF 240
C                                                                       DQSF 250
C        METHOD                                                         DQSF 260
C           BEGINNING WITH Z(1)=0, EVALUATION OF VECTOR Z IS DONE BY    DQSF 270
C           MEANS OF SIMPSONS RULE TOGETHER WITH NEWTONS 3/8 RULE OR A  DQSF 280
C           COMBINATION OF THESE TWO RULES. TRUNCATION ERROR IS OF      DQSF 290
C           ORDER H**5 (I.E. FOURTH ORDER METHOD). ONLY IN CASE NDIM=3  DQSF 300
C           TRUNCATION ERROR OF Z(2) IS OF ORDER H**4.                  DQSF 310
C           FOR REFERENCE, SEE                                          DQSF 320
C           (1) F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,     DQSF 330
C               MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.71-76.   DQSF 340
C           (2) R.ZURMUEHL, PRAKTISCHE MATHEMATIK FUER INGENIEURE UND   DQSF 350
C               PHYSIKER, SPRINGER, BERLIN/GOETTINGEN/HEIDELBERG, 1963, DQSF 360
C               PP.214-221.                                             DQSF 370
C                                                                       DQSF 380
C     ..................................................................DQSF 390
C                                                                       DQSF 400
      SUBROUTINE DQSF(H,Y,Z,NDIM)                                       DQSF 410
C                                                                       DQSF 420
C                                                                       DQSF 430
      DIMENSION Y(1),Z(1)                                               DQSF 440
      DOUBLE PRECISION Y,Z,H,HT,SUM1,SUM2,AUX,AUX1,AUX2                 DQSF 450
C                                                                       DQSF 460
      HT=.33333333333333333D0*H                                         DQSF 470
      IF(NDIM-5)7,8,1                                                   DQSF 480
C                                                                       DQSF 490
C     NDIM IS GREATER THAN 5. PREPARATIONS OF INTEGRATION LOOP          DQSF 500
    1 SUM1=Y(2)+Y(2)                                                    DQSF 510
      SUM1=SUM1+SUM1                                                    DQSF 520
      SUM1=HT*(Y(1)+SUM1+Y(3))                                          DQSF 530
      AUX1=Y(4)+Y(4)                                                    DQSF 540
      AUX1=AUX1+AUX1                                                    DQSF 550
      AUX1=SUM1+HT*(Y(3)+AUX1+Y(5))                                     DQSF 560
      AUX2=HT*(Y(1)+3.875D0*(Y(2)+Y(5))+2.625D0*(Y(3)+Y(4))+Y(6))       DQSF 570
      SUM2=Y(5)+Y(5)                                                    DQSF 580
      SUM2=SUM2+SUM2                                                    DQSF 590
      SUM2=AUX2-HT*(Y(4)+SUM2+Y(6))                                     DQSF 600
      Z(1)=0.D0                                                         DQSF 610
      AUX=Y(3)+Y(3)                                                     DQSF 620
      AUX=AUX+AUX                                                       DQSF 630
      Z(2)=SUM2-HT*(Y(2)+AUX+Y(4))                                      DQSF 640
      Z(3)=SUM1                                                         DQSF 650
      Z(4)=SUM2                                                         DQSF 660
      IF(NDIM-6)5,5,2                                                   DQSF 670
C                                                                       DQSF 680
C     INTEGRATION LOOP                                                  DQSF 690
    2 DO 4 I=7,NDIM,2                                                   DQSF 700
      SUM1=AUX1                                                         DQSF 710
      SUM2=AUX2                                                         DQSF 720
      AUX1=Y(I-1)+Y(I-1)                                                DQSF 730
      AUX1=AUX1+AUX1                                                    DQSF 740
      AUX1=SUM1+HT*(Y(I-2)+AUX1+Y(I))                                   DQSF 750
      Z(I-2)=SUM1                                                       DQSF 760
      IF(I-NDIM)3,6,6                                                   DQSF 770
    3 AUX2=Y(I)+Y(I)                                                    DQSF 780
      AUX2=AUX2+AUX2                                                    DQSF 790
      AUX2=SUM2+HT*(Y(I-1)+AUX2+Y(I+1))                                 DQSF 800
    4 Z(I-1)=SUM2                                                       DQSF 810
    5 Z(NDIM-1)=AUX1                                                    DQSF 820
      Z(NDIM)=AUX2                                                      DQSF 830
      RETURN                                                            DQSF 840
    6 Z(NDIM-1)=SUM2                                                    DQSF 850
      Z(NDIM)=AUX1                                                      DQSF 860
      RETURN                                                            DQSF 870
C     END OF INTEGRATION LOOP                                           DQSF 880
C                                                                       DQSF 890
    7 IF(NDIM-3)12,11,8                                                 DQSF 900
C                                                                       DQSF 910
C     NDIM IS EQUAL TO 4 OR 5                                           DQSF 920
    8 SUM2=1.125D0*HT*(Y(1)+Y(2)+Y(2)+Y(2)+Y(3)+Y(3)+Y(3)+Y(4))         DQSF 930
      SUM1=Y(2)+Y(2)                                                    DQSF 940
      SUM1=SUM1+SUM1                                                    DQSF 950
      SUM1=HT*(Y(1)+SUM1+Y(3))                                          DQSF 960
      Z(1)=0.D0                                                         DQSF 970
      AUX1=Y(3)+Y(3)                                                    DQSF 980
      AUX1=AUX1+AUX1                                                    DQSF 990
      Z(2)=SUM2-HT*(Y(2)+AUX1+Y(4))                                     DQSF1000
      IF(NDIM-5)10,9,9                                                  DQSF1010
    9 AUX1=Y(4)+Y(4)                                                    DQSF1020
      AUX1=AUX1+AUX1                                                    DQSF1030
      Z(5)=SUM1+HT*(Y(3)+AUX1+Y(5))                                     DQSF1040
   10 Z(3)=SUM1                                                         DQSF1050
      Z(4)=SUM2                                                         DQSF1060
      RETURN                                                            DQSF1070
C                                                                       DQSF1080
C     NDIM IS EQUAL TO 3                                                DQSF1090
   11 SUM1=HT*(1.25D0*Y(1)+Y(2)+Y(2)-.25D0*Y(3))                        DQSF1100
      SUM2=Y(2)+Y(2)                                                    DQSF1110
      SUM2=SUM2+SUM2                                                    DQSF1120
      Z(3)=HT*(Y(1)+SUM2+Y(3))                                          DQSF1130
      Z(1)=0.D0                                                         DQSF1140
      Z(2)=SUM1                                                         DQSF1150
   12 RETURN                                                            DQSF1160
      END                                                               DQSF1170