Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap2_198111 - decus/20-0026/hpcl.ssp
There are 2 other files named hpcl.ssp in the archive. Click here to see a list.
C                                                                       HPCL  10
C     ..................................................................HPCL  20
C                                                                       HPCL  30
C        SUBROUTINE HPCL                                                HPCL  40
C                                                                       HPCL  50
C        PURPOSE                                                        HPCL  60
C           TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY LINEAR            HPCL  70
C           DIFFERENTIAL EQUATIONS WITH GIVEN INITIAL VALUES.           HPCL  80
C                                                                       HPCL  90
C        USAGE                                                          HPCL 100
C           CALL HPCL (PRMT,Y,DERY,NDIM,IHLF,AFCT,FCT,OUTP,AUX,A)       HPCL 110
C           PARAMETERS AFCT,FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT. HPCL 120
C                                                                       HPCL 130
C        DESCRIPTION OF PARAMETERS                                      HPCL 140
C           PRMT   - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER  HPCL 150
C                    OR EQUAL TO 5, WHICH SPECIFIES THE PARAMETERS OF   HPCL 160
C                    THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR  HPCL 170
C                    COMMUNICATION BETWEEN OUTPUT SUBROUTINE (FURNISHED HPCL 180
C                    BY THE USER) AND SUBROUTINE HPCL. EXCEPT PRMT(5)   HPCL 190
C                    THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE     HPCL 200
C                    HPCL AND THEY ARE                                  HPCL 210
C           PRMT(1)- LOWER BOUND OF THE INTERVAL (INPUT),               HPCL 220
C           PRMT(2)- UPPER BOUND OF THE INTERVAL (INPUT),               HPCL 230
C           PRMT(3)- INITIAL INCREMENT OF THE INDEPENDENT VARIABLE      HPCL 240
C                    (INPUT),                                           HPCL 250
C           PRMT(4)- UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS    HPCL 260
C                    GREATER THAN PRMT(4), INCREMENT GETS HALVED.       HPCL 270
C                    IF INCREMENT IS LESS THAN PRMT(3) AND ABSOLUTE     HPCL 280
C                    ERROR LESS THAN PRMT(4)/50, INCREMENT GETS DOUBLED.HPCL 290
C                    THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS        HPCL 300
C                    OUTPUT SUBROUTINE.                                 HPCL 310
C           PRMT(5)- NO INPUT PARAMETER. SUBROUTINE HPCL INITIALIZES    HPCL 320
C                    PRMT(5)=0. IF THE USER WANTS TO TERMINATE          HPCL 330
C                    SUBROUTINE HPCL AT ANY OUTPUT POINT, HE HAS TO     HPCL 340
C                    CHANGE PRMT(5) TO NON-ZERO BY MEANS OF SUBROUTINE  HPCL 350
C                    OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE        HPCL 360
C                    FEASIBLE IF ITS DIMENSION IS DEFINED GREATER       HPCL 370
C                    THAN 5. HOWEVER SUBROUTINE HPCL DOES NOT REQUIRE   HPCL 380
C                    AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL   HPCL 390
C                    FOR HANDING RESULT VALUES TO THE MAIN PROGRAM      HPCL 400
C                    (CALLING HPCL) WHICH ARE OBTAINED BY SPECIAL       HPCL 410
C                    MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP. HPCL 420
C           Y      - INPUT VECTOR OF INITIAL VALUES.  (DESTROYED)       HPCL 430
C                    LATERON Y IS THE RESULTING VECTOR OF DEPENDENT     HPCL 440
C                    VARIABLES COMPUTED AT INTERMEDIATE POINTS X.       HPCL 450
C           DERY   - INPUT VECTOR OF ERROR WEIGHTS.  (DESTROYED)        HPCL 460
C                    THE SUM OF ITS COMPONENTS MUST BE EQUAL TO 1.      HPCL 470
C                    LATERON DERY IS THE VECTOR OF DERIVATIVES, WHICH   HPCL 480
C                    BELONG TO FUNCTION VALUES Y AT A POINT X.          HPCL 490
C           NDIM   - AN INPUT VALUE, WHICH SPECIFIES THE NUMBER OF      HPCL 500
C                    EQUATIONS IN THE SYSTEM.                           HPCL 510
C           IHLF   - AN OUTPUT VALUE, WHICH SPECIFIES THE NUMBER OF     HPCL 520
C                    BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS  HPCL 530
C                    GREATER THAN 10, SUBROUTINE HPCL RETURNS WITH      HPCL 540
C                    ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM.           HPCL 550
C                    ERROR MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE   HPCL 560
C                    PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)-HPCL 570
C                    PRMT(1)) RESPECTIVELY.                             HPCL 580
C           AFCT   - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT        HPCL 590
C                    COMPUTES MATRIX A (FACTOR OF VECTOR Y ON THE       HPCL 600
C                    RIGHT HAND SIDE OF THE SYSTEM) FOR A GIVEN X-VALUE.HPCL 610
C                    ITS PARAMETER LIST MUST BE X,A. THE SUBROUTINE     HPCL 620
C                    SHOULD NOT DESTROY X.                              HPCL 630
C           FCT    - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT        HPCL 640
C                    COMPUTES VECTOR F (INHOMOGENEOUS PART OF THE       HPCL 650
C                    RIGHT HAND SIDE OF THE SYSTEM) FOR A GIVEN X-VALUE.HPCL 660
C                    ITS PARAMETER LIST MUST BE X,F. THE SUBROUTINE     HPCL 670
C                    SHOULD NOT DESTROY X.                              HPCL 680
C           OUTP   - THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED.    HPCL 690
C                    ITS PARAMETER LIST MUST BE X,Y,DERY,IHLF,NDIM,PRMT.HPCL 700
C                    NONE OF THESE PARAMETERS (EXCEPT, IF NECESSARY,    HPCL 710
C                    PRMT(4),PRMT(5),...) SHOULD BE CHANGED BY          HPCL 720
C                    SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO,HPCL 730
C                    SUBROUTINE HPCL IS TERMINATED.                     HPCL 740
C           AUX    - AN AUXILIARY STORAGE ARRAY WITH 16 ROWS AND NDIM   HPCL 750
C                    COLUMNS.                                           HPCL 760
C           A      - AN NDIM BY NDIM MATRIX, WHICH IS USED AS AUXILIARY HPCL 770
C                    STORAGE ARRAY.                                     HPCL 780
C                                                                       HPCL 790
C        REMARKS                                                        HPCL 800
C           THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM, IF HPCL 810
C           (1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE    HPCL 820
C               NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE   HPCL 830
C               IHLF=11),                                               HPCL 840
C           (2) INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN       HPCL 850
C               (ERROR MESSAGES IHLF=12 OR IHLF=13),                    HPCL 860
C           (3) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH,       HPCL 870
C           (4) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO.        HPCL 880
C                                                                       HPCL 890
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  HPCL 900
C           THE EXTERNAL SUBROUTINES AFCT(X,A), FCT(X,F) AND            HPCL 910
C           OUTP(X,Y,DERY,IHLF,NDIM,PRMT) MUST BE FURNISHED BY THE USER.HPCL 920
C                                                                       HPCL 930
C        METHOD                                                         HPCL 940
C           EVALUATION IS DONE BY MEANS OF HAMMINGS MODIFIED PREDICTOR- HPCL 950
C           CORRECTOR METHOD. IT IS A FOURTH ORDER METHOD, USING 4      HPCL 960
C           PRECEEDING POINTS FOR COMPUTATION OF A NEW VECTOR Y OF THE  HPCL 970
C           DEPENDENT VARIABLES.                                        HPCL 980
C           FOURTH ORDER RUNGE-KUTTA METHOD SUGGESTED BY RALSTON IS     HPCL 990
C           USED FOR ADJUSTMENT OF THE INITIAL INCREMENT AND FOR        HPCL1000
C           COMPUTATION OF STARTING VALUES.                             HPCL1010
C           SUBROUTINE HPCL AUTOMATICALLY ADJUSTS THE INCREMENT DURING  HPCL1020
C           THE WHOLE COMPUTATION BY HALVING OR DOUBLING.               HPCL1030
C           TO GET FULL FLEXIBILITY IN OUTPUT, AN OUTPUT SUBROUTINE     HPCL1040
C           MUST BE CODED BY THE USER.                                  HPCL1050
C           FOR REFERENCE, SEE                                          HPCL1060
C           (1)  RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL         HPCL1070
C                COMPUTERS, WILEY, NEW YORK/LONDON, 1960, PP.95-109.    HPCL1080
C           (2)  RALSTON, RUNGE-KUTTA METHODS WITH MINIMUM ERROR BOUNDS,HPCL1090
C                MTAC, VOL.16, ISS.80 (1962), PP.431-437.               HPCL1100
C                                                                       HPCL1110
C     ..................................................................HPCL1120
C                                                                       HPCL1130
      SUBROUTINE HPCL(PRMT,Y,DERY,NDIM,IHLF,AFCT,FCT,OUTP,AUX,A)        HPCL1140
C                                                                       HPCL1150
C                                                                       HPCL1160
C     THE FOLLOWING FIRST PART OF SUBROUTINE HPCL (UNTIL FIRST BREAK-   HPCL1170
C     POINT FOR LINKAGE) HAS TO STAY IN CORE DURING THE WHOLE           HPCL1180
C     COMPUTATION                                                       HPCL1190
C                                                                       HPCL1200
      DIMENSION PRMT(1),Y(1),DERY(1),AUX(16,1),A(1)                     HPCL1210
      GOTO 100                                                          HPCL1220
C                                                                       HPCL1230
C     THIS PART OF SUBROUTINE HPCL COMPUTES THE RIGHT HAND SIDE DERY OF HPCL1240
C     THE GIVEN SYSTEM OF LINEAR DIFFERENTIAL EQUATIONS.                HPCL1250
    1 CALL AFCT(X,A)                                                    HPCL1260
      CALL FCT(X,DERY)                                                  HPCL1270
      DO 3 M=1,NDIM                                                     HPCL1280
      LL=M-NDIM                                                         HPCL1290
      HS=0.                                                             HPCL1300
      DO 2 L=1,NDIM                                                     HPCL1310
      LL=LL+NDIM                                                        HPCL1320
    2 HS=HS+A(LL)*Y(L)                                                  HPCL1330
    3 DERY(M)=HS+DERY(M)                                                HPCL1340
      GOTO(105,202,204,206,115,122,125,308,312,327,329,128),ISW2        HPCL1350
C                                                                       HPCL1360
C     POSSIBLE BREAK-POINT FOR LINKAGE                                  HPCL1370
C                                                                       HPCL1380
  100 N=1                                                               HPCL1390
      IHLF=0                                                            HPCL1400
      X=PRMT(1)                                                         HPCL1410
      H=PRMT(3)                                                         HPCL1420
      PRMT(5)=0.                                                        HPCL1430
      DO 101 I=1,NDIM                                                   HPCL1440
      AUX(16,I)=0.                                                      HPCL1450
      AUX(15,I)=DERY(I)                                                 HPCL1460
  101 AUX(1,I)=Y(I)                                                     HPCL1470
      IF(H*(PRMT(2)-X))103,102,104                                      HPCL1480
C                                                                       HPCL1490
C     ERROR RETURNS                                                     HPCL1500
  102 IHLF=12                                                           HPCL1510
      GOTO 104                                                          HPCL1520
  103 IHLF=13                                                           HPCL1530
C                                                                       HPCL1540
C     COMPUTATION OF DERY FOR STARTING VALUES                           HPCL1550
  104 ISW2=1                                                            HPCL1560
      GOTO 1                                                            HPCL1570
C                                                                       HPCL1580
C     RECORDING OF STARTING VALUES                                      HPCL1590
  105 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                HPCL1600
      IF(PRMT(5))107,106,107                                            HPCL1610
  106 IF(IHLF)108,108,107                                               HPCL1620
  107 RETURN                                                            HPCL1630
  108 DO 109 I=1,NDIM                                                   HPCL1640
  109 AUX(8,I)=DERY(I)                                                  HPCL1650
C                                                                       HPCL1660
C     COMPUTATION OF AUX(2,I)                                           HPCL1670
      ISW1=1                                                            HPCL1680
      GOTO 200                                                          HPCL1690
C                                                                       HPCL1700
  110 X=X+H                                                             HPCL1710
      DO 111 I=1,NDIM                                                   HPCL1720
  111 AUX(2,I)=Y(I)                                                     HPCL1730
C                                                                       HPCL1740
C     INCREMENT H IS TESTED BY MEANS OF BISECTION                       HPCL1750
  112 IHLF=IHLF+1                                                       HPCL1760
      X=X-H                                                             HPCL1770
      DO 113 I=1,NDIM                                                   HPCL1780
  113 AUX(4,I)=AUX(2,I)                                                 HPCL1790
      H=.5*H                                                            HPCL1800
      N=1                                                               HPCL1810
      ISW1=2                                                            HPCL1820
      GOTO 200                                                          HPCL1830
C                                                                       HPCL1840
  114 X=X+H                                                             HPCL1850
      ISW2=5                                                            HPCL1860
      GOTO 1                                                            HPCL1870
  115 N=2                                                               HPCL1880
      DO 116 I=1,NDIM                                                   HPCL1890
      AUX(2,I)=Y(I)                                                     HPCL1900
  116 AUX(9,I)=DERY(I)                                                  HPCL1910
      ISW1=3                                                            HPCL1920
      GOTO 200                                                          HPCL1930
C                                                                       HPCL1940
C     COMPUTATION OF TEST VALUE DELT                                    HPCL1950
  117 DELT=0.                                                           HPCL1960
      DO 118 I=1,NDIM                                                   HPCL1970
  118 DELT=DELT+AUX(15,I)*ABS(Y(I)-AUX(4,I))                            HPCL1980
      DELT=.06666667*DELT                                               HPCL1990
      IF(DELT-PRMT(4))121,121,119                                       HPCL2000
  119 IF(IHLF-10)112,120,120                                            HPCL2010
C                                                                       HPCL2020
C     NO SATISFACTORY ACCURACY AFTER 10 BISECTIONS. ERROR MESSAGE.      HPCL2030
  120 IHLF=11                                                           HPCL2040
      X=X+H                                                             HPCL2050
      GOTO 104                                                          HPCL2060
C                                                                       HPCL2070
C     SATISFACTORY ACCURACY AFTER LESS THAN 11 BISECTIONS               HPCL2080
  121 X=X+H                                                             HPCL2090
      ISW2=6                                                            HPCL2100
      GOTO 1                                                            HPCL2110
  122 DO 123 I=1,NDIM                                                   HPCL2120
      AUX(3,I)=Y(I)                                                     HPCL2130
  123 AUX(10,I)=DERY(I)                                                 HPCL2140
      N=3                                                               HPCL2150
      ISW1=4                                                            HPCL2160
      GOTO 200                                                          HPCL2170
C                                                                       HPCL2180
  124 N=1                                                               HPCL2190
      X=X+H                                                             HPCL2200
      ISW2=7                                                            HPCL2210
      GOTO 1                                                            HPCL2220
  125 X=PRMT(1)                                                         HPCL2230
      DO 126 I=1,NDIM                                                   HPCL2240
      AUX(11,I)=DERY(I)                                                 HPCL2250
  1260Y(I)=AUX(1,I)+H*(.375*AUX(8,I)+.7916667*AUX(9,I)                  HPCL2260
     1-.2083333*AUX(10,I)+.04166667*DERY(I))                            HPCL2270
  127 X=X+H                                                             HPCL2280
      N=N+1                                                             HPCL2290
      ISW2=12                                                           HPCL2300
      GOTO 1                                                            HPCL2310
  128 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                HPCL2320
      IF(PRMT(5))107,129,107                                            HPCL2330
  129 IF(N-4)130,300,300                                                HPCL2340
  130 DO 131 I=1,NDIM                                                   HPCL2350
      AUX(N,I)=Y(I)                                                     HPCL2360
  131 AUX(N+7,I)=DERY(I)                                                HPCL2370
      IF(N-3)132,134,300                                                HPCL2380
C                                                                       HPCL2390
  132 DO 133 I=1,NDIM                                                   HPCL2400
      DELT=AUX(9,I)+AUX(9,I)                                            HPCL2410
      DELT=DELT+DELT                                                    HPCL2420
  133 Y(I)=AUX(1,I)+.3333333*H*(AUX(8,I)+DELT+AUX(10,I))                HPCL2430
      GOTO 127                                                          HPCL2440
C                                                                       HPCL2450
  134 DO 135 I=1,NDIM                                                   HPCL2460
      DELT=AUX(9,I)+AUX(10,I)                                           HPCL2470
      DELT=DELT+DELT+DELT                                               HPCL2480
  135 Y(I)=AUX(1,I)+.375*H*(AUX(8,I)+DELT+AUX(11,I))                    HPCL2490
      GOTO 127                                                          HPCL2500
C                                                                       HPCL2510
C     THE FOLLOWING PART OF SUBROUTINE HPCL COMPUTES BY MEANS OF        HPCL2520
C     RUNGE-KUTTA METHOD STARTING VALUES FOR THE NOT SELF-STARTING      HPCL2530
C     PREDICTOR-CORRECTOR METHOD.                                       HPCL2540
  200 Z=X                                                               HPCL2550
      DO 201 I=1,NDIM                                                   HPCL2560
      X=H*AUX(N+7,I)                                                    HPCL2570
      AUX(5,I)=X                                                        HPCL2580
  201 Y(I)=AUX(N,I)+.4*X                                                HPCL2590
C     X IS AN AUXILIARY STORAGE LOCATION                                HPCL2600
C                                                                       HPCL2610
      X=Z+.4*H                                                          HPCL2620
      ISW2=2                                                            HPCL2630
      GOTO 1                                                            HPCL2640
  202 DO 203 I=1,NDIM                                                   HPCL2650
      X=H*DERY(I)                                                       HPCL2660
      AUX(6,I)=X                                                        HPCL2670
  203 Y(I)=AUX(N,I)+.2969776*AUX(5,I)+.1587596*X                        HPCL2680
C                                                                       HPCL2690
      X=Z+.4557372*H                                                    HPCL2700
      ISW2=3                                                            HPCL2710
      GOTO 1                                                            HPCL2720
  204 DO 205 I=1,NDIM                                                   HPCL2730
      X=H*DERY(I)                                                       HPCL2740
      AUX(7,I)=X                                                        HPCL2750
  205 Y(I)=AUX(N,I)+.2181004*AUX(5,I)-3.050965*AUX(6,I)+3.832865*X      HPCL2760
C                                                                       HPCL2770
      X=Z+H                                                             HPCL2780
      ISW2=4                                                            HPCL2790
      GOTO 1                                                            HPCL2800
  206 DO 207 I=1,NDIM                                                   HPCL2810
  2070Y(I)=AUX(N,I)+.1747603*AUX(5,I)-.5514807*AUX(6,I)                 HPCL2820
     1+1.205536*AUX(7,I)+.1711848*H*DERY(I)                             HPCL2830
      X=Z                                                               HPCL2840
      GOTO(110,114,117,124),ISW1                                        HPCL2850
C                                                                       HPCL2860
C     POSSIBLE BREAK-POINT FOR LINKAGE                                  HPCL2870
C                                                                       HPCL2880
C     STARTING VALUES ARE COMPUTED.                                     HPCL2890
C     NOW START HAMMINGS MODIFIED PREDICTOR-CORRECTOR METHOD.           HPCL2900
  300 ISTEP=3                                                           HPCL2910
  301 IF(N-8)304,302,304                                                HPCL2920
C                                                                       HPCL2930
C     N=8 CAUSES THE ROWS OF AUX TO CHANGE THEIR STORAGE LOCATIONS      HPCL2940
  302 DO 303 N=2,7                                                      HPCL2950
      DO 303 I=1,NDIM                                                   HPCL2960
      AUX(N-1,I)=AUX(N,I)                                               HPCL2970
  303 AUX(N+6,I)=AUX(N+7,I)                                             HPCL2980
      N=7                                                               HPCL2990
C                                                                       HPCL3000
C     N LESS THAN 8 CAUSES N+1 TO GET N                                 HPCL3010
  304 N=N+1                                                             HPCL3020
C                                                                       HPCL3030
C     COMPUTATION OF NEXT VECTOR Y                                      HPCL3040
      DO 305 I=1,NDIM                                                   HPCL3050
      AUX(N-1,I)=Y(I)                                                   HPCL3060
  305 AUX(N+6,I)=DERY(I)                                                HPCL3070
      X=X+H                                                             HPCL3080
  306 ISTEP=ISTEP+1                                                     HPCL3090
      DO 307 I=1,NDIM                                                   HPCL3100
     0DELT=AUX(N-4,I)+1.333333*H*(AUX(N+6,I)+AUX(N+6,I)-AUX(N+5,I)+     HPCL3110
     1AUX(N+4,I)+AUX(N+4,I))                                            HPCL3120
      Y(I)=DELT-.9256198*AUX(16,I)                                      HPCL3130
  307 AUX(16,I)=DELT                                                    HPCL3140
C     PREDICTOR IS NOW GENERATED IN ROW 16 OF AUX, MODIFIED PREDICTOR   HPCL3150
C     IS GENERATED IN Y. DELT MEANS AN AUXILIARY STORAGE.               HPCL3160
      ISW2=8                                                            HPCL3170
      GOTO 1                                                            HPCL3180
C     DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY             HPCL3190
C                                                                       HPCL3200
  308 DO 309 I=1,NDIM                                                   HPCL3210
     0DELT=.125*(9.*AUX(N-1,I)-AUX(N-3,I)+3.*H*(DERY(I)+AUX(N+6,I)+     HPCL3220
     1AUX(N+6,I)-AUX(N+5,I)))                                           HPCL3230
      AUX(16,I)=AUX(16,I)-DELT                                          HPCL3240
  309 Y(I)=DELT+.07438017*AUX(16,I)                                     HPCL3250
C                                                                       HPCL3260
C     TEST WHETHER H MUST BE HALVED OR DOUBLED                          HPCL3270
      DELT=0.                                                           HPCL3280
      DO 310 I=1,NDIM                                                   HPCL3290
  310 DELT=DELT+AUX(15,I)*ABS(AUX(16,I))                                HPCL3300
      IF(DELT-PRMT(4))311,324,324                                       HPCL3310
C                                                                       HPCL3320
C     H MUST NOT BE HALVED. THAT MEANS Y(I) ARE GOOD.                   HPCL3330
  311 ISW2=9                                                            HPCL3340
      GOTO 1                                                            HPCL3350
  312 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                HPCL3360
      IF(PRMT(5))314,313,314                                            HPCL3370
  313 IF(IHLF-11)315,314,314                                            HPCL3380
  314 RETURN                                                            HPCL3390
  315 IF(H*(X-PRMT(2)))316,314,314                                      HPCL3400
  316 IF(ABS(X-PRMT(2))-.1*ABS(H))314,317,317                           HPCL3410
  317 IF(DELT-.02*PRMT(4))318,318,301                                   HPCL3420
C                                                                       HPCL3430
C                                                                       HPCL3440
C     H COULD BE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE         HPCL3450
C     AVAILABLE                                                         HPCL3460
  318 IF(IHLF)301,301,319                                               HPCL3470
  319 IF(N-7)301,320,320                                                HPCL3480
  320 IF(ISTEP-4)301,321,321                                            HPCL3490
  321 IMOD=ISTEP/2                                                      HPCL3500
      IF(ISTEP-IMOD-IMOD)301,322,301                                    HPCL3510
  322 H=H+H                                                             HPCL3520
      IHLF=IHLF-1                                                       HPCL3530
      ISTEP=0                                                           HPCL3540
      DO 323 I=1,NDIM                                                   HPCL3550
      AUX(N-1,I)=AUX(N-2,I)                                             HPCL3560
      AUX(N-2,I)=AUX(N-4,I)                                             HPCL3570
      AUX(N-3,I)=AUX(N-6,I)                                             HPCL3580
      AUX(N+6,I)=AUX(N+5,I)                                             HPCL3590
      AUX(N+5,I)=AUX(N+3,I)                                             HPCL3600
      AUX(N+4,I)=AUX(N+1,I)                                             HPCL3610
      DELT=AUX(N+6,I)+AUX(N+5,I)                                        HPCL3620
      DELT=DELT+DELT+DELT                                               HPCL3630
  3230AUX(16,I)=8.962963*(Y(I)-AUX(N-3,I))-3.361111*H*(DERY(I)+DELT     HPCL3640
     1+AUX(N+4,I))                                                      HPCL3650
      GOTO 301                                                          HPCL3660
C                                                                       HPCL3670
C                                                                       HPCL3680
C     H MUST BE HALVED                                                  HPCL3690
  324 IHLF=IHLF+1                                                       HPCL3700
      IF(IHLF-10)325,325,311                                            HPCL3710
  325 H=.5*H                                                            HPCL3720
      ISTEP=0                                                           HPCL3730
      DO 326 I=1,NDIM                                                   HPCL3740
     0Y(I)=.00390625*(80.*AUX(N-1,I)+135.*AUX(N-2,I)+40.*AUX(N-3,I)+    HPCL3750
     1AUX(N-4,I))-.1171875*(AUX(N+6,I)-6.*AUX(N+5,I)-AUX(N+4,I))*H      HPCL3760
     0AUX(N-4,I)=.00390625*(12.*AUX(N-1,I)+135.*AUX(N-2,I)+             HPCL3770
     1108.*AUX(N-3,I)+AUX(N-4,I))-.0234375*(AUX(N+6,I)+18.*AUX(N+5,I)-  HPCL3780
     29.*AUX(N+4,I))*H                                                  HPCL3790
      AUX(N-3,I)=AUX(N-2,I)                                             HPCL3800
  326 AUX(N+4,I)=AUX(N+5,I)                                             HPCL3810
      DELT=X-H                                                          HPCL3820
      X=DELT-(H+H)                                                      HPCL3830
      ISW2=10                                                           HPCL3840
      GOTO 1                                                            HPCL3850
  327 DO 328 I=1,NDIM                                                   HPCL3860
      AUX(N-2,I)=Y(I)                                                   HPCL3870
      AUX(N+5,I)=DERY(I)                                                HPCL3880
  328 Y(I)=AUX(N-4,I)                                                   HPCL3890
      X=X-(H+H)                                                         HPCL3900
      ISW2=11                                                           HPCL3910
      GOTO 1                                                            HPCL3920
  329 X=DELT                                                            HPCL3930
      DO 330 I=1,NDIM                                                   HPCL3940
      DELT=AUX(N+5,I)+AUX(N+4,I)                                        HPCL3950
      DELT=DELT+DELT+DELT                                               HPCL3960
     0AUX(16,I)=8.962963*(AUX(N-1,I)-Y(I))-3.361111*H*(AUX(N+6,I)+DELT  HPCL3970
     1+DERY(I))                                                         HPCL3980
  330 AUX(N+3,I)=DERY(I)                                                HPCL3990
      GOTO 306                                                          HPCL4000
      END                                                               HPCL4010