Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/drtmi.ssp
There are 2 other files named drtmi.ssp in the archive. Click here to see a list.
C                                                                       DRTM  10
C     ..................................................................DRTM  20
C                                                                       DRTM  30
C        SUBROUTINE DRTMI                                               DRTM  40
C                                                                       DRTM  50
C        PURPOSE                                                        DRTM  60
C           TO SOLVE GENERAL NONLINEAR EQUATIONS OF THE FORM FCT(X)=0   DRTM  70
C           BY MEANS OF MUELLER-S ITERATION METHOD.                     DRTM  80
C                                                                       DRTM  90
C        USAGE                                                          DRTM 100
C           CALL DRTMI (X,F,FCT,XLI,XRI,EPS,IEND,IER)                   DRTM 110
C           PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT.               DRTM 120
C                                                                       DRTM 130
C        DESCRIPTION OF PARAMETERS                                      DRTM 140
C           X      - DOUBLE PRECISION RESULTANT ROOT OF EQUATION        DRTM 150
C                    FCT(X)=0.                                          DRTM 160
C           F      - DOUBLE PRECISION RESULTANT FUNCTION VALUE          DRTM 170
C                    AT ROOT X.                                         DRTM 180
C           FCT    - NAME OF THE EXTERNAL DOUBLE PRECISION FUNCTION     DRTM 190
C                    SUBPROGRAM USED.                                   DRTM 200
C           XLI    - DOUBLE PRECISION INPUT VALUE WHICH SPECIFIES THE   DRTM 210
C                    INITIAL LEFT BOUND OF THE ROOT X.                  DRTM 220
C           XRI    - DOUBLE PRECISION INPUT VALUE WHICH SPECIFIES THE   DRTM 230
C                    INITIAL RIGHT BOUND OF THE ROOT X.                 DRTM 240
C           EPS    - SINGLE PRECISION INPUT VALUE WHICH SPECIFIES THE   DRTM 250
C                    UPPER BOUND OF THE ERROR OF RESULT X.              DRTM 260
C           IEND   - MAXIMUM NUMBER OF ITERATION STEPS SPECIFIED.       DRTM 270
C           IER    - RESULTANT ERROR PARAMETER CODED AS FOLLOWS         DRTM 280
C                     IER=0 - NO ERROR,                                 DRTM 290
C                     IER=1 - NO CONVERGENCE AFTER IEND ITERATION STEPS DRTM 300
C                             FOLLOWED BY IEND SUCCESSIVE STEPS OF      DRTM 310
C                             BISECTION,                                DRTM 320
C                     IER=2 - BASIC ASSUMPTION FCT(XLI)*FCT(XRI) LESS   DRTM 330
C                             THAN OR EQUAL TO ZERO IS NOT SATISFIED.   DRTM 340
C                                                                       DRTM 350
C        REMARKS                                                        DRTM 360
C           THE PROCEDURE ASSUMES THAT FUNCTION VALUES AT INITIAL       DRTM 370
C           BOUNDS XLI AND XRI HAVE NOT THE SAME SIGN. IF THIS BASIC    DRTM 380
C           ASSUMPTION IS NOT SATISFIED BY INPUT VALUES XLI AND XRI, THEDRTM 390
C           PROCEDURE IS BYPASSED AND GIVES THE ERROR MESSAGE IER=2.    DRTM 400
C                                                                       DRTM 410
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DRTM 420
C           THE EXTERNAL DOUBLE PRECISION FUNCTION SUBPROGRAM FCT(X)    DRTM 430
C           MUST BE FURNISHED BY THE USER.                              DRTM 440
C                                                                       DRTM 450
C        METHOD                                                         DRTM 460
C           SOLUTION OF EQUATION FCT(X)=0 IS DONE BY MEANS OF MUELLER-S DRTM 470
C           ITERATION METHOD OF SUCCESSIVE BISECTIONS AND INVERSE       DRTM 480
C           PARABOLIC INTERPOLATION, WHICH STARTS AT THE INITIAL BOUNDS DRTM 490
C           XLI AND XRI. CONVERGENCE IS QUADRATIC IF THE DERIVATIVE OF  DRTM 500
C           FCT(X) AT ROOT X IS NOT EQUAL TO ZERO. ONE ITERATION STEP   DRTM 510
C           REQUIRES TWO EVALUATIONS OF FCT(X). FOR TEST ON SATISFACTORYDRTM 520
C           ACCURACY SEE FORMULAE (3,4) OF MATHEMATICAL DESCRIPTION.    DRTM 530
C           FOR REFERENCE, SEE G. K. KRISTIANSEN, ZERO OF ARBITRARY     DRTM 540
C           FUNCTION, BIT, VOL. 3 (1963), PP.205-206.                   DRTM 550
C                                                                       DRTM 560
C     ..................................................................DRTM 570
C                                                                       DRTM 580
      SUBROUTINE DRTMI(X,F,FCT,XLI,XRI,EPS,IEND,IER)                    DRTM 590
C                                                                       DRTM 600
C                                                                       DRTM 610
      DOUBLE PRECISION X,F,FCT,XLI,XRI,XL,XR,FL,FR,TOL,TOLF,A,DX,XM,FM  DRTM 620
C                                                                       DRTM 630
C     PREPARE ITERATION                                                 DRTM 640
      IER=0                                                             DRTM 650
      XL=XLI                                                            DRTM 660
      XR=XRI                                                            DRTM 670
      X=XL                                                              DRTM 680
      TOL=X                                                             DRTM 690
      F=FCT(TOL)                                                        DRTM 700
      IF(F)1,16,1                                                       DRTM 710
    1 FL=F                                                              DRTM 720
      X=XR                                                              DRTM 730
      TOL=X                                                             DRTM 740
      F=FCT(TOL)                                                        DRTM 750
      IF(F)2,16,2                                                       DRTM 760
    2 FR=F                                                              DRTM 770
      IF(DSIGN(1.D0,FL)+DSIGN(1.D0,FR))25,3,25                          DRTM 780
C                                                                       DRTM 790
C     BASIC ASSUMPTION FL*FR LESS THAN 0 IS SATISFIED.                  DRTM 800
C     GENERATE TOLERANCE FOR FUNCTION VALUES.                           DRTM 810
    3 I=0                                                               DRTM 820
      TOLF=100.*EPS                                                     DRTM 830
C                                                                       DRTM 840
C                                                                       DRTM 850
C     START ITERATION LOOP                                              DRTM 860
    4 I=I+1                                                             DRTM 870
C                                                                       DRTM 880
C     START BISECTION LOOP                                              DRTM 890
      DO 13 K=1,IEND                                                    DRTM 900
      X=.5D0*(XL+XR)                                                    DRTM 910
      TOL=X                                                             DRTM 920
      F=FCT(TOL)                                                        DRTM 930
      IF(F)5,16,5                                                       DRTM 940
    5 IF(DSIGN(1.D0,F)+DSIGN(1.D0,FR))7,6,7                             DRTM 950
C                                                                       DRTM 960
C     INTERCHANGE XL AND XR IN ORDER TO GET THE SAME SIGN IN F AND FR   DRTM 970
    6 TOL=XL                                                            DRTM 980
      XL=XR                                                             DRTM 990
      XR=TOL                                                            DRTM1000
      TOL=FL                                                            DRTM1010
      FL=FR                                                             DRTM1020
      FR=TOL                                                            DRTM1030
    7 TOL=F-FL                                                          DRTM1040
      A=F*TOL                                                           DRTM1050
      A=A+A                                                             DRTM1060
      IF(A-FR*(FR-FL))8,9,9                                             DRTM1070
    8 IF(I-IEND)17,17,9                                                 DRTM1080
    9 XR=X                                                              DRTM1090
      FR=F                                                              DRTM1100
C                                                                       DRTM1110
C     TEST ON SATISFACTORY ACCURACY IN BISECTION LOOP                   DRTM1120
      TOL=EPS                                                           DRTM1130
      A=DABS(XR)                                                        DRTM1140
      IF(A-1.D0)11,11,10                                                DRTM1150
   10 TOL=TOL*A                                                         DRTM1160
   11 IF(DABS(XR-XL)-TOL)12,12,13                                       DRTM1170
   12 IF(DABS(FR-FL)-TOLF)14,14,13                                      DRTM1180
   13 CONTINUE                                                          DRTM1190
C     END OF BISECTION LOOP                                             DRTM1200
C                                                                       DRTM1210
C     NO CONVERGENCE AFTER IEND ITERATION STEPS FOLLOWED BY IEND        DRTM1220
C     SUCCESSIVE STEPS OF BISECTION OR STEADILY INCREASING FUNCTION     DRTM1230
C     VALUES AT RIGHT BOUNDS. ERROR RETURN.                             DRTM1240
      IER=1                                                             DRTM1250
   14 IF(DABS(FR)-DABS(FL))16,16,15                                     DRTM1260
   15 X=XL                                                              DRTM1270
      F=FL                                                              DRTM1280
   16 RETURN                                                            DRTM1290
C                                                                       DRTM1300
C     COMPUTATION OF ITERATED X-VALUE BY INVERSE PARABOLIC INTERPOLATIONDRTM1310
   17 A=FR-F                                                            DRTM1320
      DX=(X-XL)*FL*(1.D0+F*(A-TOL)/(A*(FR-FL)))/TOL                     DRTM1330
      XM=X                                                              DRTM1340
      FM=F                                                              DRTM1350
      X=XL-DX                                                           DRTM1360
      TOL=X                                                             DRTM1370
      F=FCT(TOL)                                                        DRTM1380
      IF(F)18,16,18                                                     DRTM1390
C                                                                       DRTM1400
C     TEST ON SATISFACTORY ACCURACY IN ITERATION LOOP                   DRTM1410
   18 TOL=EPS                                                           DRTM1420
      A=DABS(X)                                                         DRTM1430
      IF(A-1.D0)20,20,19                                                DRTM1440
   19 TOL=TOL*A                                                         DRTM1450
   20 IF(DABS(DX)-TOL)21,21,22                                          DRTM1460
   21 IF(DABS(F)-TOLF)16,16,22                                          DRTM1470
C                                                                       DRTM1480
C     PREPARATION OF NEXT BISECTION LOOP                                DRTM1490
   22 IF(DSIGN(1.D0,F)+DSIGN(1.D0,FL))24,23,24                          DRTM1500
   23 XR=X                                                              DRTM1510
      FR=F                                                              DRTM1520
      GO TO 4                                                           DRTM1530
   24 XL=X                                                              DRTM1540
      FL=F                                                              DRTM1550
      XR=XM                                                             DRTM1560
      FR=FM                                                             DRTM1570
      GO TO 4                                                           DRTM1580
C     END OF ITERATION LOOP                                             DRTM1590
C                                                                       DRTM1600
C                                                                       DRTM1610
C     ERROR RETURN IN CASE OF WRONG INPUT DATA                          DRTM1620
   25 IER=2                                                             DRTM1630
      RETURN                                                            DRTM1640
      END                                                               DRTM1650