Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap2_198111 - decus/20-0026/dlbvp.ssp
There are 2 other files named dlbvp.ssp in the archive. Click here to see a list.
C                                                                       DLBV  10
C     ..................................................................DLBV  20
C                                                                       DLBV  30
C        SUBROUTINE DLBVP                                               DLBV  40
C                                                                       DLBV  50
C        PURPOSE                                                        DLBV  60
C           TO SOLVE A LINEAR BOUNDARY VALUE PROBLEM, WHICH CONSISTS OF DLBV  70
C           A SYSTEM OF NDIM LINEAR FIRST ORDER DIFFERENTIAL EQUATIONS  DLBV  80
C                  DY/DX=A(X)*Y(X)+F(X)                                 DLBV  90
C           AND NDIM LINEAR BOUNDARY CONDITIONS                         DLBV 100
C                  B*Y(XL)+C*Y(XU)=R.                                   DLBV 110
C                                                                       DLBV 120
C        USAGE                                                          DLBV 130
C           CALL DLBVP (PRMT,B,C,R,Y,DERY,NDIM,IHLF,AFCT,FCT,DFCT,OUTP, DLBV 140
C                      AUX,A)                                           DLBV 150
C           PARAMETERS AFCT,FCT,DFCT,OUTP REQUIRE AN EXTERNAL STATEMENT.DLBV 160
C                                                                       DLBV 170
C        DESCRIPTION OF PARAMETERS                                      DLBV 180
C           PRMT   - DOUBLE PRECISION INPUT AND OUTPUT VECTOR WITH      DLBV 190
C                    DIMENSION GREATER THAN OR EQUAL TO 5, WHICH        DLBV 200
C                    SPECIFIES THE PARAMETERS OF THE INTERVAL AND OF    DLBV 210
C                    ACCURACY AND WHICH SERVES FOR COMMUNICATION BETWEENDLBV 220
C                    OUTPUT SUBROUTINE (FURNISHED BY THE USER) AND      DLBV 230
C                    SUBROUTINE DLBVP. EXCEPT PRMT(5) THE COMPONENTS    DLBV 240
C                    ARE NOT DESTROYED BY SUBROUTINE DLBVP AND THEY ARE DLBV 250
C           PRMT(1)- LOWER BOUND XL OF THE INTERVAL (INPUT),            DLBV 260
C           PRMT(1)- UPPER BOUND XU OF THE INTERVAL (INPUT),            DLBV 270
C           PRMT(3)- INITIAL INCREMENT OF THE INDEPENDENT VARIABLE      DLBV 280
C                    (INPUT),                                           DLBV 290
C           PRMT(4)- UPPER ERROR BOUND (INPUT). IF RELATIVE ERROR IS    DLBV 300
C                    GREATER THAN PRMT(4), INCREMENT GETS HALVED.       DLBV 310
C                    IF INCREMENT IS LESS THAN PRMT(3) AND RELATIVE     DLBV 320
C                    ERROR LESS THAN PRMT(4)/50, INCREMENT GETS DOUBLED.DLBV 330
C                    THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS        DLBV 340
C                    OUTPUT SUBROUTINE.                                 DLBV 350
C           PRMT(5)- NO INPUT PARAMETER. SUBROUTINE DLBVP INITIALIZES   DLBV 360
C                    PRMT(5)=0. IF THE USER WANTS TO TERMINATE          DLBV 370
C                    SUBROUTINE DLBVP AT ANY OUTPUT POINT, HE HAS TO    DLBV 380
C                    CHANGE PRMT(5) TO NON-ZERO BY MEANS OF SUBROUTINE  DLBV 390
C                    OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE        DLBV 400
C                    FEASIBLE IF ITS DIMENSION IS DEFINED GREATER       DLBV 410
C                    THAN 5. HOWEVER SUBROUTINE DLBVP DOES NOT REQUIRE  DLBV 420
C                    AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL   DLBV 430
C                    FOR HANDING RESULT VALUES TO THE MAIN PROGRAM      DLBV 440
C                    (CALLING DLBVP) WHICH ARE OBTAINED BY SPECIAL      DLBV 450
C                    MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP. DLBV 460
C           B      - DOUBLE PRECISION NDIM BY NDIM INPUT MATRIX         DLBV 470
C                    (DESTROYED). IT IS THE COEFFICIENT MATRIX OF Y(XL) DLBV 480
C                    IN THE BOUNDARY CONDITIONS.                        DLBV 490
C           C      - DOUBLE PRECISION NDIM BY NDIM INPUT MATRIX         DLBV 500
C                    (POSSIBLY DESTROYED). IT IS THE COEFFICIENT MATRIX DLBV 510
C                    OF Y(XU) IN THE BOUNDARY CONDITIONS.               DLBV 520
C           R      - DOUBLE PRECISION INPUT VECTOR WITH DIMENSION NDIM  DLBV 530
C                    (DESTROYED). IT SPECIFIES THE RIGHT HAND SIDE OF   DLBV 540
C                    THE BOUNDARY CONDITIONS.                           DLBV 550
C           Y      - DOUBLE PRECISION AUXILIARY VECTOR WITH             DLBV 560
C                    DIMENSION NDIM. IT IS USED AS STORAGE LOCATION     DLBV 570
C                    FOR THE RESULTING VALUES OF DEPENDENT VARIABLES    DLBV 580
C                    COMPUTED AT INTERMEDIATE POINTS X.                 DLBV 590
C           DERY   - DOUBLE PRECISION INPUT VECTOR OF ERROR WEIGHTS     DLBV 600
C                    (DESTROYED). ITS MAXIMAL COMPONENT SHOULD BE       DLBV 610
C                    EQUAL TO 1. LATERON DERY IS THE VECTOR OF          DLBV 620
C                    DERIVATIVES, WHICH BELONG TO FUNCTION VALUES Y AT  DLBV 630
C                    INTERMEDIATE POINTS X.                             DLBV 640
C           NDIM   - AN INPUT VALUE, WHICH SPECIFIES THE NUMBER OF      DLBV 650
C                    DIFFERENTIAL EQUATIONS IN THE SYSTEM.              DLBV 660
C           IHLF   - AN OUTPUT VALUE, WHICH SPECIFIES THE NUMBER OF     DLBV 670
C                    BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS  DLBV 680
C                    GREATER THAN 10, SUBROUTINE DLBVP RETURNS WITH     DLBV 690
C                    ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM.           DLBV 700
C                    ERROR MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE   DLBV 710
C                    PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)-DLBV 720
C                    PRMT(1)) RESPECTIVELY. FINALLY ERROR MESSAGE       DLBV 730
C                    IHLF=14 INDICATES, THAT THERE IS NO SOLUTION OR    DLBV 740
C                    THAT THERE ARE MORE THAN ONE SOLUTION OF THE       DLBV 750
C                    PROBLEM.                                           DLBV 760
C                    A NEGATIVE VALUE OF IHLF HANDED TO SUBROUTINE OUTP DLBV 770
C                    TOGETHER WITH INITIAL VALUES OF FINALLY GENERATED  DLBV 780
C                    INITIAL VALUE PROBLEM INDICATES, THAT THERE WAS    DLBV 790
C                    POSSIBLE LOSS OF SIGNIFICANCE IN THE SOLUTION OF   DLBV 800
C                    THE SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS FOR    DLBV 810
C                    THESE INITIAL VALUES. THE ABSOLUTE VALUE OF IHLF   DLBV 820
C                    SHOWS, AFTER WHICH ELIMINATION STEP OF GAUSS       DLBV 830
C                    ALGORITHM POSSIBLE LOSS OF SIGNIFICANCE WAS        DLBV 840
C                    DETECTED.                                          DLBV 850
C           AFCT   - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT        DLBV 860
C                    COMPUTES THE COEFFICIENT MATRIX A OF VECTOR Y ON   DLBV 870
C                    THE RIGHT HAND SIDE OF THE SYSTEM OF DIFFERENTIAL  DLBV 880
C                    EQUATIONS FOR A GIVEN X-VALUE. ITS PARAMETER LIST  DLBV 890
C                    MUST BE X,A. SUBROUTINE AFCT SHOULD NOT DESTROY X. DLBV 900
C           FCT    - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT        DLBV 910
C                    COMPUTES VECTOR F (INHOMOGENEOUS PART OF THE       DLBV 920
C                    RIGHT HAND SIDE OF THE SYSTEM OF DIFFERENTIAL      DLBV 930
C                    EQUATIONS) FOR A GIVEN X-VALUE. ITS PARAMETER LIST DLBV 940
C                    MUST BE X,F. SUBROUTINE FCT SHOULD NOT DESTROY X.  DLBV 950
C           DFCT   - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT        DLBV 960
C                    COMPUTES VECTOR DF (DERIVATIVE OF THE INHOMOGENEOUSDLBV 970
C                    PART ON THE RIGHT HAND SIDE OF THE SYSTEM OF       DLBV 980
C                    DIFFERENTIAL EQUATIONS) FOR A GIVEN X-VALUE. ITS   DLBV 990
C                    PARAMETER LIST MUST BE X,DF. SUBROUTINE DFCT       DLBV1000
C                    SHOULD NOT DESTROY X.                              DLBV1010
C           OUTP   - THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED.    DLBV1020
C                    ITS PARAMETER LIST MUST BE X,Y,DERY,IHLF,NDIM,PRMT.DLBV1030
C                    NONE OF THESE PARAMETERS (EXCEPT, IF NECESSARY,    DLBV1040
C                    PRMT(4),PRMT(5),...) SHOULD BE CHANGED BY          DLBV1050
C                    SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO,DLBV1060
C                    SUBROUTINE DLBVP IS TERMINATED.                    DLBV1070
C           AUX    - DOUBLE PRECISION AUXILIARY STORAGE ARRAY WITH 20   DLBV1080
C                    ROWS AND NDIM COLUMNS.                             DLBV1090
C           A      - DOUBLE PRECISION NDIM BY NDIM MATRIX, WHICH IS USEDDLBV1100
C                    AS AUXILIARY STORAGE ARRAY.                        DLBV1110
C                                                                       DLBV1120
C        REMARKS                                                        DLBV1130
C           THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM, IF DLBV1140
C           (1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE    DLBV1150
C               NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE   DLBV1160
C               IHLF=11),                                               DLBV1170
C           (2) INITIAL INCREMENT IS EQUAL TO 0 OR IF IT HAS WRONG SIGN DLBV1180
C               (ERROR MESSAGES IHLF=12 OR IHLF=13),                    DLBV1190
C           (3) THERE IS NO OR MORE THAN ONE SOLUTION OF THE PROBLEM    DLBV1200
C               (ERROR MESSAGE IHLF=14),                                DLBV1210
C           (4) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH,       DLBV1220
C           (5) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO.        DLBV1230
C                                                                       DLBV1240
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DLBV1250
C           SUBROUTINE DGELG     SYSTEM OF LINEAR EQUATIONS.            DLBV1260
C           THE EXTERNAL SUBROUTINES AFCT(X,A), FCT(X,F), DFCT(X,DF),   DLBV1270
C           AND OUTP(X,Y,DERY,IHLF,NDIM,PRMT) MUST BE FURNISHED         DLBV1280
C           BY THE USER.                                                DLBV1290
C                                                                       DLBV1300
C        METHOD                                                         DLBV1310
C           EVALUATION IS DONE USING THE METHOD OF ADJOINT EQUATIONS.   DLBV1320
C           HAMMINGS FOURTH ORDER MODIFIED PREDICTOR-CORRECTOR METHOD   DLBV1330
C           IS USED TO SOLVE THE ADJOINT INITIAL VALUE PROBLEMS AND FI- DLBV1340
C           NALLY TO SOLVE THE GENERATED INITIAL VALUE PROBLEM FOR Y(X).DLBV1350
C           THE INITIAL INCREMENT PRMT(3) IS AUTOMATICALLY ADJUSTED.    DLBV1360
C           FOR COMPUTATION OF INTEGRAL SUM, A FOURTH ORDER HERMITEAN   DLBV1370
C           INTEGRATION FORMULA IS USED.                                DLBV1380
C           FOR REFERENCE, SEE                                          DLBV1390
C           (1) LANCE, NUMERICAL METHODS FOR HIGH SPEED COMPUTERS,      DLBV1400
C               ILIFFE, LONDON, 1960, PP.64-67.                         DLBV1410
C           (2) RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL          DLBV1420
C               COMPUTERS, WILEY, NEW YORK/LONDON, 1960, PP.95-109.     DLBV1430
C           (3) RALSTON, RUNGE-KUTTA METHODS WITH MINIMUM ERROR BOUNDS, DLBV1440
C               MTAC, VOL.16, ISS.80 (1962), PP.431-437.                DLBV1450
C           (4) ZURMUEHL, PRAKTISCHE MATHEMATIK FUER INGENIEURE UND     DLBV1460
C               PHYSIKER, SPRINGER, BERLIN/GOETTINGEN/HEIDELBERG, 1963, DLBV1470
C               PP.227-232.                                             DLBV1480
C                                                                       DLBV1490
C     ..................................................................DLBV1500
C                                                                       DLBV1510
     0SUBROUTINE DLBVP(PRMT,B,C,R,Y,DERY,NDIM,IHLF,AFCT,FCT,DFCT,OUTP,  DLBV1520
     1AUX,A)                                                            DLBV1530
C                                                                       DLBV1540
C                                                                       DLBV1550
      DIMENSION PRMT(1),B(1),C(1),R(1),Y(1),DERY(1),AUX(20,1),A(1)      DLBV1560
      DOUBLE PRECISION PRMT,B,C,R,Y,DERY,AUX,A,H,X,Z,GL,HS,GU,SUM,      DLBV1570
     1DGL,DGU,XST,XEND,DELT                                             DLBV1580
C                                                                       DLBV1590
C     ERROR TEST                                                        DLBV1600
      IF(PRMT(3)*(PRMT(2)-PRMT(1)))2,1,3                                DLBV1610
    1 IHLF=12                                                           DLBV1620
      RETURN                                                            DLBV1630
    2 IHLF=13                                                           DLBV1640
      RETURN                                                            DLBV1650
C                                                                       DLBV1660
C     SEARCH FOR ZERO-COLUMNS IN MATRICES B AND C                       DLBV1670
    3 KK=-NDIM                                                          DLBV1680
      IB=0                                                              DLBV1690
      IC=0                                                              DLBV1700
      DO 7 K=1,NDIM                                                     DLBV1710
      AUX(15,K)=DERY(K)                                                 DLBV1720
      AUX(1,K)=1.D0                                                     DLBV1730
      AUX(17,K)=1.D0                                                    DLBV1740
      KK=KK+NDIM                                                        DLBV1750
      DO 4 I=1,NDIM                                                     DLBV1760
      II=KK+I                                                           DLBV1770
      IF(B(II))5,4,5                                                    DLBV1780
    4 CONTINUE                                                          DLBV1790
      IB=IB+1                                                           DLBV1800
      AUX(1,K)=0.D0                                                     DLBV1810
    5 DO 6 I=1,NDIM                                                     DLBV1820
      II=KK+I                                                           DLBV1830
      IF(C(II))7,6,7                                                    DLBV1840
    6 CONTINUE                                                          DLBV1850
      IC=IC+1                                                           DLBV1860
      AUX(17,K)=0.D0                                                    DLBV1870
    7 CONTINUE                                                          DLBV1880
C                                                                       DLBV1890
C     DETERMINATION OF LOWER AND UPPER BOUND                            DLBV1900
      IF(IC-IB)8,11,11                                                  DLBV1910
    8 H=PRMT(2)                                                         DLBV1920
      PRMT(2)=PRMT(1)                                                   DLBV1930
      PRMT(1)=H                                                         DLBV1940
      PRMT(3)=-PRMT(3)                                                  DLBV1950
      DO 9 I=1,NDIM                                                     DLBV1960
    9 AUX(17,I)=AUX(1,I)                                                DLBV1970
      II=NDIM*NDIM                                                      DLBV1980
      DO 10 I=1,II                                                      DLBV1990
      H=B(I)                                                            DLBV2000
      B(I)=C(I)                                                         DLBV2010
   10 C(I)=H                                                            DLBV2020
C                                                                       DLBV2030
C     PREPARATIONS FOR CONSTRUCTION OF ADJOINT INITIAL VALUE PROBLEMS   DLBV2040
   11 X=PRMT(2)                                                         DLBV2050
      CALL FCT(X,Y)                                                     DLBV2060
      CALL DFCT(X,DERY)                                                 DLBV2070
      DO 12 I=1,NDIM                                                    DLBV2080
      AUX(18,I)=Y(I)                                                    DLBV2090
   12 AUX(19,I)=DERY(I)                                                 DLBV2100
C                                                                       DLBV2110
C     POSSIBLE BREAK-POINT FOR LINKAGE                                  DLBV2120
C                                                                       DLBV2130
C     THE FOLLOWING PART OF SUBROUTINE DLBVP UNTIL NEXT BREAK-POINT FOR DLBV2140
C     LINKAGE HAS TO REMAIN IN CORE DURING THE WHOLE REST OF THE        DLBV2150
C     COMPUTATIONS                                                      DLBV2160
C                                                                       DLBV2170
C     START LOOP FOR GENERATING ADJOINT INITIAL VALUE PROBLEMS          DLBV2180
      K=0                                                               DLBV2190
      KK=0                                                              DLBV2200
  100 K=K+1                                                             DLBV2210
      IF(AUX(17,K))108,108,101                                          DLBV2220
C                                                                       DLBV2230
C     INITIALIZATION OF ADJOINT INITIAL VALUE PROBLEM                   DLBV2240
  101 X=PRMT(2)                                                         DLBV2250
      CALL AFCT(X,A)                                                    DLBV2260
      SUM=0.D0                                                          DLBV2270
      GL=AUX(18,K)                                                      DLBV2280
      DGL=AUX(19,K)                                                     DLBV2290
      II=K                                                              DLBV2300
      DO 104 I=1,NDIM                                                   DLBV2310
      H=-A(II)                                                          DLBV2320
      DERY(I)=H                                                         DLBV2330
      AUX(20,I)=R(I)                                                    DLBV2340
      Y(I)=0.D0                                                         DLBV2350
      IF(I-K)103,102,103                                                DLBV2360
  102 Y(I)=1.D0                                                         DLBV2370
  103 DGL=DGL+H*AUX(18,I)                                               DLBV2380
  104 II=II+NDIM                                                        DLBV2390
      XEND=PRMT(1)                                                      DLBV2400
      H=.0625D0*(XEND-X)                                                DLBV2410
      ISW=0                                                             DLBV2420
      GOTO 400                                                          DLBV2430
C     THIS IS BRANCH TO ADJOINT LINEAR INITIAL VALUE PROBLEM            DLBV2440
C                                                                       DLBV2450
C     THIS IS RETURN FROM ADJOINT LINEAR INITIAL VALUE PROBLEM          DLBV2460
  105 IF(IHLF-10)106,106,117                                            DLBV2470
C                                                                       DLBV2480
C     UPDATING OF COEFFICIENT MATRIX B AND VECTOR R                     DLBV2490
  106 DO 107 I=1,NDIM                                                   DLBV2500
      KK=KK+1                                                           DLBV2510
      H=C(KK)                                                           DLBV2520
      R(I)=AUX(20,I)+H*SUM                                              DLBV2530
      II=I                                                              DLBV2540
      DO 107 J=1,NDIM                                                   DLBV2550
      B(II)=B(II)+H*Y(J)                                                DLBV2560
  107 II=II+NDIM                                                        DLBV2570
      GOTO 109                                                          DLBV2580
  108 KK=KK+NDIM                                                        DLBV2590
  109 IF(K-NDIM)100,110,110                                             DLBV2600
C                                                                       DLBV2610
C                                                                       DLBV2620
C     GENERATION OF LAST INITIAL VALUE PROBLEM                          DLBV2630
  110 EPS=PRMT(4)                                                       DLBV2640
      CALL DGELG(R,B,NDIM,1,EPS,I)                                      DLBV2650
      IF(I)111,112,112                                                  DLBV2660
  111 IHLF=14                                                           DLBV2670
      RETURN                                                            DLBV2680
C                                                                       DLBV2690
  112 PRMT(5)=0.D0                                                      DLBV2700
      IHLF=-I                                                           DLBV2710
      X=PRMT(1)                                                         DLBV2720
      XEND=PRMT(2)                                                      DLBV2730
      H=PRMT(3)                                                         DLBV2740
      DO 113 I=1,NDIM                                                   DLBV2750
  113 Y(I)=R(I)                                                         DLBV2760
      ISW=1                                                             DLBV2770
  114 ISW2=12                                                           DLBV2780
      GOTO 200                                                          DLBV2790
  115 ISW3=-1                                                           DLBV2800
      GOTO 300                                                          DLBV2810
  116 IF(IHLF)400,400,117                                               DLBV2820
C     THIS WAS BRANCH INTO INITIAL VALUE PROBLEM                        DLBV2830
C                                                                       DLBV2840
C     THIS IS RETURN FROM INITIAL VALUE PROBLEM                         DLBV2850
  117 RETURN                                                            DLBV2860
C                                                                       DLBV2870
C     THIS PART OF LINEAR BOUNDARY VALUE PROBLEM COMPUTES THE RIGHT     DLBV2880
C     HAND SIDE DERY OF THE SYSTEM OF ADJOINT LINEAR DIFFERENTIAL       DLBV2890
C     EQUATIONS (IN CASE ISW=0) OR OF THE GIVEN SYSTEM (IN CASE ISW=1). DLBV2900
  200 CALL AFCT(X,A)                                                    DLBV2910
      IF(ISW)201,201,205                                                DLBV2920
C                                                                       DLBV2930
C     ADJOINT SYSTEM                                                    DLBV2940
  201 LL=0                                                              DLBV2950
      DO 203 M=1,NDIM                                                   DLBV2960
      HS=0.D0                                                           DLBV2970
      DO 202 L=1,NDIM                                                   DLBV2980
      LL=LL+1                                                           DLBV2990
  202 HS=HS-A(LL)*Y(L)                                                  DLBV3000
  203 DERY(M)=HS                                                        DLBV3010
  204 GOTO(502,504,506,407,415,418,608,617,632,634,421,115),ISW2        DLBV3020
C                                                                       DLBV3030
C     GIVEN SYSTEM                                                      DLBV3040
  205 CALL FCT(X,DERY)                                                  DLBV3050
      DO 207 M=1,NDIM                                                   DLBV3060
      LL=M-NDIM                                                         DLBV3070
      HS=0.D0                                                           DLBV3080
      DO 206 L=1,NDIM                                                   DLBV3090
      LL=LL+NDIM                                                        DLBV3100
  206 HS=HS+A(LL)*Y(L)                                                  DLBV3110
  207 DERY(M)=HS+DERY(M)                                                DLBV3120
      GOTO 204                                                          DLBV3130
C                                                                       DLBV3140
C     THIS PART OF LINEAR BOUNDARY VALUE PROBLEM COMPUTES THE VALUE OF  DLBV3150
C     INTEGRAL SUM, WHICH IS A PART OF THE OUTPUT OF ADJOINT INITIAL    DLBV3160
C     VALUE PROBLEM (IN CASE ISW=0) OR RECORDS RESULT VALUES OF THE     DLBV3170
C     FINAL INITIAL VALUE PROBLEM (IN CASE ISW=1).                      DLBV3180
  300 IF(ISW)301,301,305                                                DLBV3190
C                                                                       DLBV3200
C     ADJOINT PROBLEM                                                   DLBV3210
  301 CALL FCT(X,R)                                                     DLBV3220
      GU=0.D0                                                           DLBV3230
      DGU=0.D0                                                          DLBV3240
      DO 302 L=1,NDIM                                                   DLBV3250
      GU=GU+Y(L)*R(L)                                                   DLBV3260
  302 DGU=DGU+DERY(L)*R(L)                                              DLBV3270
      CALL DFCT(X,R)                                                    DLBV3280
      DO 303 L=1,NDIM                                                   DLBV3290
  303 DGU=DGU+Y(L)*R(L)                                                 DLBV3300
      SUM=SUM+.5D0*H*((GL+GU)+.16666666666666667D0*H*(DGL-DGU))         DLBV3310
      GL=GU                                                             DLBV3320
      DGL=DGU                                                           DLBV3330
  304 IF(ISW3)116,422,618                                               DLBV3340
C                                                                       DLBV3350
C     GIVEN PROBLEM                                                     DLBV3360
  305 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                DLBV3370
      IF(PRMT(5))117,304,117                                            DLBV3380
C                                                                       DLBV3390
C     POSSIBLE BREAK-POINT FOR LINKAGE                                  DLBV3400
C                                                                       DLBV3410
C     THE FOLLOWING PART OF SUBROUTINE DLBVP SOLVES IN CASE ISW=0 THE   DLBV3420
C     ADJOINT INITIAL VALUE PROBLEM. IT COMPUTES INTEGRAL SUM AND       DLBV3430
C     THE VECTOR Y OF DEPENDENT VARIABLES AT THE LOWER BOUND PRMT(1).   DLBV3440
C     IN CASE ISW=1 IT SOLVES FINALLY GENERATED INITIAL VALUE PROBLEM.  DLBV3450
  400 N=1                                                               DLBV3460
      XST=X                                                             DLBV3470
      IHLF=0                                                            DLBV3480
      DO 401 I=1,NDIM                                                   DLBV3490
      AUX(16,I)=0.D0                                                    DLBV3500
      AUX(1,I)=Y(I)                                                     DLBV3510
  401 AUX(8,I)=DERY(I)                                                  DLBV3520
      ISW1=1                                                            DLBV3530
      GOTO 500                                                          DLBV3540
C                                                                       DLBV3550
  402 X=X+H                                                             DLBV3560
      DO 403 I=1,NDIM                                                   DLBV3570
  403 AUX(2,I)=Y(I)                                                     DLBV3580
C                                                                       DLBV3590
C     INCREMENT H IS TESTED BY MEANS OF BISECTION                       DLBV3600
  404 IHLF=IHLF+1                                                       DLBV3610
      X=X-H                                                             DLBV3620
      DO 405 I=1,NDIM                                                   DLBV3630
  405 AUX(4,I)=AUX(2,I)                                                 DLBV3640
      H=.5D0*H                                                          DLBV3650
      N=1                                                               DLBV3660
      ISW1=2                                                            DLBV3670
      GOTO 500                                                          DLBV3680
C                                                                       DLBV3690
  406 X=X+H                                                             DLBV3700
      ISW2=4                                                            DLBV3710
      GOTO 200                                                          DLBV3720
  407 N=2                                                               DLBV3730
      DO 408 I=1,NDIM                                                   DLBV3740
      AUX(2,I)=Y(I)                                                     DLBV3750
  408 AUX(9,I)=DERY(I)                                                  DLBV3760
      ISW1=3                                                            DLBV3770
      GOTO 500                                                          DLBV3780
C                                                                       DLBV3790
C     TEST ON SATISFACTORY ACCURACY                                     DLBV3800
  409 DO 414 I=1,NDIM                                                   DLBV3810
      Z=DABS(Y(I))                                                      DLBV3820
      IF(Z-1.D0)410,411,411                                             DLBV3830
  410 Z=1.D0                                                            DLBV3840
  411 DELT=.066666666666666667D0*DABS(Y(I)-AUX(4,I))                    DLBV3850
      IF(ISW)413,413,412                                                DLBV3860
  412 DELT=AUX(15,I)*DELT                                               DLBV3870
  413 IF(DELT-Z*PRMT(4))414,414,429                                     DLBV3880
  414 CONTINUE                                                          DLBV3890
C                                                                       DLBV3900
C     SATISFACTORY ACCURACY AFTER LESS THAN 11 BISECTIONS               DLBV3910
      X=X+H                                                             DLBV3920
      ISW2=5                                                            DLBV3930
      GOTO 200                                                          DLBV3940
  415 DO 416 I=1,NDIM                                                   DLBV3950
      AUX(3,I)=Y(I)                                                     DLBV3960
  416 AUX(10,I)=DERY(I)                                                 DLBV3970
      N=3                                                               DLBV3980
      ISW1=4                                                            DLBV3990
      GOTO 500                                                          DLBV4000
C                                                                       DLBV4010
  417 N=1                                                               DLBV4020
      X=X+H                                                             DLBV4030
      ISW2=6                                                            DLBV4040
      GOTO 200                                                          DLBV4050
  418 X=XST                                                             DLBV4060
      DO 419 I=1,NDIM                                                   DLBV4070
      AUX(11,I)=DERY(I)                                                 DLBV4080
  4190Y(I)=AUX(1,I)+H*(.375D0*AUX(8,I)+.7916666666666667D0*AUX(9,I)     DLBV4090
     1-.20833333333333333D0*AUX(10,I)+.041666666666666667D0*DERY(I))    DLBV4100
  420 X=X+H                                                             DLBV4110
      N=N+1                                                             DLBV4120
      ISW2=11                                                           DLBV4130
      GOTO 200                                                          DLBV4140
  421 ISW3=0                                                            DLBV4150
      GOTO 300                                                          DLBV4160
  422 IF(N-4)423,600,600                                                DLBV4170
  423 DO 424 I=1,NDIM                                                   DLBV4180
      AUX(N,I)=Y(I)                                                     DLBV4190
  424 AUX(N+7,I)=DERY(I)                                                DLBV4200
      IF(N-3)425,427,600                                                DLBV4210
C                                                                       DLBV4220
  425 DO 426 I=1,NDIM                                                   DLBV4230
      DELT=AUX(9,I)+AUX(9,I)                                            DLBV4240
      DELT=DELT+DELT                                                    DLBV4250
  426 Y(I)=AUX(1,I)+.33333333333333333D0*H*(AUX(8,I)+DELT+AUX(10,I))    DLBV4260
      GOTO 420                                                          DLBV4270
C                                                                       DLBV4280
  427 DO 428 I=1,NDIM                                                   DLBV4290
      DELT=AUX(9,I)+AUX(10,I)                                           DLBV4300
      DELT=DELT+DELT+DELT                                               DLBV4310
  428 Y(I)=AUX(1,I)+.375D0*H*(AUX(8,I)+DELT+AUX(11,I))                  DLBV4320
      GOTO 420                                                          DLBV4330
C                                                                       DLBV4340
C     NO SATISFACTORY ACCURACY. H MUST BE HALVED.                       DLBV4350
  429 IF(IHLF-10)404,430,430                                            DLBV4360
C                                                                       DLBV4370
C     NO SATISFACTORY ACCURACY AFTER 10 BISECTIONS. ERROR MESSAGE.      DLBV4380
  430 IHLF=11                                                           DLBV4390
      X=X+H                                                             DLBV4400
      IF(ISW)105,105,114                                                DLBV4410
C                                                                       DLBV4420
C     THIS PART OF LINEAR INITIAL VALUE PROBLEM COMPUTES                DLBV4430
C     STARTING VALUES BY MEANS OF RUNGE-KUTTA METHOD.                   DLBV4440
  500 Z=X                                                               DLBV4450
      DO 501 I=1,NDIM                                                   DLBV4460
      X=H*AUX(N+7,I)                                                    DLBV4470
      AUX(5,I)=X                                                        DLBV4480
  501 Y(I)=AUX(N,I)+.4D0*X                                              DLBV4490
C                                                                       DLBV4500
      X=Z+.4D0*H                                                        DLBV4510
      ISW2=1                                                            DLBV4520
      GOTO 200                                                          DLBV4530
  502 DO 503 I=1,NDIM                                                   DLBV4540
      X=H*DERY(I)                                                       DLBV4550
      AUX(6,I)=X                                                        DLBV4560
  503 Y(I)=AUX(N,I)+.29697760924775360D0*AUX(5,I)+.15875964497103583D0*XDLBV4570
C                                                                       DLBV4580
      X=Z+.45573725421878943D0*H                                        DLBV4590
      ISW2=2                                                            DLBV4600
      GOTO 200                                                          DLBV4610
  504 DO 505 I=1,NDIM                                                   DLBV4620
      X=H*DERY(I)                                                       DLBV4630
      AUX(7,I)=X                                                        DLBV4640
  505 Y(I)=AUX(N,I)+.21810038822592047D0*AUX(5,I)-3.0509651486929308D0* DLBV4650
     1AUX(6,I)+3.8328647604670103D0*X                                   DLBV4660
C                                                                       DLBV4670
      X=Z+H                                                             DLBV4680
      ISW2=3                                                            DLBV4690
      GOTO 200                                                          DLBV4700
  506 DO 507 I=1,NDIM                                                   DLBV4710
  5070Y(I)=AUX(N,I)+.17476028226269037D0*AUX(5,I)-.55148066287873294D0* DLBV4720
     1AUX(6,I)+1.2055355993965235D0*AUX(7,I)+.17118478121951903D0*      DLBV4730
     2H*DERY(I)                                                         DLBV4740
      X=Z                                                               DLBV4750
      GOTO(402,406,409,417),ISW1                                        DLBV4760
C                                                                       DLBV4770
C     POSSIBLE BREAK-POINT FOR LINKAGE                                  DLBV4780
C                                                                       DLBV4790
C     STARTING VALUES ARE COMPUTED.                                     DLBV4800
C     NOW START HAMMINGS MODIFIED PREDICTOR-CORRECTOR METHOD.           DLBV4810
  600 ISTEP=3                                                           DLBV4820
  601 IF(N-8)604,602,604                                                DLBV4830
C                                                                       DLBV4840
C     N=8 CAUSES THE ROWS OF AUX TO CHANGE THEIR STORAGE LOCATIONS      DLBV4850
  602 DO 603 N=2,7                                                      DLBV4860
      DO 603 I=1,NDIM                                                   DLBV4870
      AUX(N-1,I)=AUX(N,I)                                               DLBV4880
  603 AUX(N+6,I)=AUX(N+7,I)                                             DLBV4890
      N=7                                                               DLBV4900
C                                                                       DLBV4910
C     N LESS THAN 8 CAUSES N+1 TO GET N                                 DLBV4920
  604 N=N+1                                                             DLBV4930
C                                                                       DLBV4940
C     COMPUTATION OF NEXT VECTOR Y                                      DLBV4950
      DO 605 I=1,NDIM                                                   DLBV4960
      AUX(N-1,I)=Y(I)                                                   DLBV4970
  605 AUX(N+6,I)=DERY(I)                                                DLBV4980
      X=X+H                                                             DLBV4990
  606 ISTEP=ISTEP+1                                                     DLBV5000
      DO 607 I=1,NDIM                                                   DLBV5010
     0DELT=AUX(N-4,I)+1.3333333333333333D0*H*(AUX(N+6,I)+AUX(N+6,I)-    DLBV5020
     1AUX(N+5,I)+AUX(N+4,I)+AUX(N+4,I))                                 DLBV5030
      Y(I)=DELT-.9256198347107438D0*AUX(16,I)                           DLBV5040
  607 AUX(16,I)=DELT                                                    DLBV5050
C     PREDICTOR IS NOW GENERATED IN ROW 16 OF AUX, MODIFIED PREDICTOR   DLBV5060
C     IS GENERATED IN Y. DELT MEANS AN AUXILIARY STORAGE.               DLBV5070
C                                                                       DLBV5080
      ISW2=7                                                            DLBV5090
      GOTO 200                                                          DLBV5100
C     DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY.            DLBV5110
C                                                                       DLBV5120
  608 DO 609 I=1,NDIM                                                   DLBV5130
     0DELT=.125D0*(9.D0*AUX(N-1,I)-AUX(N-3,I)+3.D0*H*(DERY(I)+AUX(N+6,I)DLBV5140
     1+AUX(N+6,I)-AUX(N+5,I)))                                          DLBV5150
      AUX(16,I)=AUX(16,I)-DELT                                          DLBV5160
  609 Y(I)=DELT+.07438016528925620D0*AUX(16,I)                          DLBV5170
C                                                                       DLBV5180
C     TEST WHETHER H MUST BE HALVED OR DOUBLED                          DLBV5190
      DELT=0.D0                                                         DLBV5200
      DO 616 I=1,NDIM                                                   DLBV5210
      Z=DABS(Y(I))                                                      DLBV5220
      IF(Z-1.D0)610,611,611                                             DLBV5230
  610 Z=1.D0                                                            DLBV5240
  611 Z=DABS(AUX(16,I))/Z                                               DLBV5250
      IF(ISW)613,613,612                                                DLBV5260
  612 Z=AUX(15,I)*Z                                                     DLBV5270
  613 IF(Z-PRMT(4))614,614,628                                          DLBV5280
  614 IF(DELT-Z)615,616,616                                             DLBV5290
  615 DELT=Z                                                            DLBV5300
  616 CONTINUE                                                          DLBV5310
C                                                                       DLBV5320
C     H MUST NOT BE HALVED. THAT MEANS Y(I) ARE GOOD.                   DLBV5330
      ISW2=8                                                            DLBV5340
      GOTO 200                                                          DLBV5350
  617 ISW3=1                                                            DLBV5360
      GOTO 300                                                          DLBV5370
  618 IF(H*(X-XEND))619,621,621                                         DLBV5380
  619 IF(DABS(X-XEND)-.1D0*DABS(H))621,620,620                          DLBV5390
  620 IF(DELT-.02D0*PRMT(4))622,622,601                                 DLBV5400
  621 IF(ISW)105,105,117                                                DLBV5410
C                                                                       DLBV5420
C     H COULD BE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE         DLBV5430
C     AVAILABLE.                                                        DLBV5440
  622 IF(IHLF)601,601,623                                               DLBV5450
  623 IF(N-7)601,624,624                                                DLBV5460
  624 IF(ISTEP-4)601,625,625                                            DLBV5470
  625 IMOD=ISTEP/2                                                      DLBV5480
      IF(ISTEP-IMOD-IMOD)601,626,601                                    DLBV5490
  626 H=H+H                                                             DLBV5500
      IHLF=IHLF-1                                                       DLBV5510
      ISTEP=0                                                           DLBV5520
      DO 627 I=1,NDIM                                                   DLBV5530
      AUX(N-1,I)=AUX(N-2,I)                                             DLBV5540
      AUX(N-2,I)=AUX(N-4,I)                                             DLBV5550
      AUX(N-3,I)=AUX(N-6,I)                                             DLBV5560
      AUX(N+6,I)=AUX(N+5,I)                                             DLBV5570
      AUX(N+5,I)=AUX(N+3,I)                                             DLBV5580
      AUX(N+4,I)=AUX(N+1,I)                                             DLBV5590
      DELT=AUX(N+6,I)+AUX(N+5,I)                                        DLBV5600
      DELT=DELT+DELT+DELT                                               DLBV5610
  6270AUX(16,I)=8.962962962962963D0*(Y(I)-AUX(N-3,I))                   DLBV5620
     1-3.3611111111111111D0*H*(DERY(I)+DELT+AUX(N+4,I))                 DLBV5630
      GOTO 601                                                          DLBV5640
C                                                                       DLBV5650
C     H MUST BE HALVED                                                  DLBV5660
  628 IHLF=IHLF+1                                                       DLBV5670
      IF(IHLF-10)630,630,629                                            DLBV5680
  629 IF(ISW)105,105,114                                                DLBV5690
  630 H=.5D0*H                                                          DLBV5700
      ISTEP=0                                                           DLBV5710
      DO 631 I=1,NDIM                                                   DLBV5720
     0Y(I)=.390625D-2*(8.D1*AUX(N-1,I)+135.D0*AUX(N-2,I)+4.D1*AUX(N-3,I)DLBV5730
     1+AUX(N-4,I))-.1171875D0*(AUX(N+6,I)-6.D0*AUX(N+5,I)-AUX(N+4,I))*H DLBV5740
     0AUX(N-4,I)=.390625D-2*(12.D0*AUX(N-1,I)+135.D0*AUX(N-2,I)+        DLBV5750
     1108.D0*AUX(N-3,I)+AUX(N-4,I))-.0234375D0*(AUX(N+6,I)+             DLBV5760
     218.D0*AUX(N+5,I)-9.D0*AUX(N+4,I))*H                               DLBV5770
      AUX(N-3,I)=AUX(N-2,I)                                             DLBV5780
  631 AUX(N+4,I)=AUX(N+5,I)                                             DLBV5790
      DELT=X-H                                                          DLBV5800
      X=DELT-(H+H)                                                      DLBV5810
      ISW2=9                                                            DLBV5820
      GOTO 200                                                          DLBV5830
  632 DO 633 I=1,NDIM                                                   DLBV5840
      AUX(N-2,I)=Y(I)                                                   DLBV5850
      AUX(N+5,I)=DERY(I)                                                DLBV5860
  633 Y(I)=AUX(N-4,I)                                                   DLBV5870
      X=X-(H+H)                                                         DLBV5880
      ISW2=10                                                           DLBV5890
      GOTO 200                                                          DLBV5900
  634 X=DELT                                                            DLBV5910
      DO 635 I=1,NDIM                                                   DLBV5920
      DELT=AUX(N+5,I)+AUX(N+4,I)                                        DLBV5930
      DELT=DELT+DELT+DELT                                               DLBV5940
     0AUX(16,I)=8.962962962962963D0*(AUX(N-1,I)-Y(I))                   DLBV5950
     1-3.3611111111111111D0*H*(AUX(N+6,I)+DELT+DERY(I))                 DLBV5960
  635 AUX(N+3,I)=DERY(I)                                                DLBV5970
      GOTO 606                                                          DLBV5980
C     END OF INITIAL VALUE PROBLEM                                      DLBV5990
      END                                                               DLBV6000