share/hompack/fortran/root.f

   1       SUBROUTINE ROOT(T,FT,B,C,RELERR,ABSERR,IFLAG)
   2 C
   3 C  ROOT COMPUTES A ROOT OF THE NONLINEAR EQUATION F(X)=0
   4 C  WHERE F(X) IS A CONTINOUS REAL FUNCTION OF A SINGLE REAL
   5 C  VARIABLE X.  THE METHOD USED IS A COMBINATION OF BISECTION
   6 C  AND THE SECANT RULE.
   7 C
   8 C  NORMAL INPUT CONSISTS OF A CONTINUOS FUNCTION F AND AN
   9 C  INTERVAL (B,C) SUCH THAT F(B)*F(C).LE.0.0.  EACH ITERATION
  10 C  FINDS NEW VALUES OF B AND C SUCH THAT THE INTERVAL(B,C) IS
  11 C  SHRUNK AND F(B)*F(C).LE.0.0.  THE STOPPING CRITERION IS
  12 C
  13 C          DABS(B-C).LE.2.0*(RELERR*DABS(B)+ABSERR)
  14 C
  15 C  WHERE RELERR=RELATIVE ERROR AND ABSERR=ABSOLUTE ERROR ARE
  16 C  INPUT QUANTITIES.  SET THE FLAG, IFLAG, POSITIVE TO INITIALIZE
  17 C  THE COMPUTATION.  AS B,C AND IFLAG ARE USED FOR BOTH INPUT AND
  18 C  OUTPUT, THEY MUST BE VARIABLES IN THE CALLING PROGRAM.
  19 C
  20 C  IF 0 IS A POSSIBLE ROOT, ONE SHOULD NOT CHOOSE ABSERR=0.0.
  21 C
  22 C  THE OUTPUT VALUE OF B IS THE BETTER APPROXIMATION TO A ROOT
  23 C  AS B AND C ARE ALWAYS REDEFINED SO THAT DABS(F(B)).LE.DABS(F(C)).
  24 C
  25 C  TO SOLVE THE EQUATION, ROOT MUST EVALUATE F(X) REPEATEDLY. THIS
  26 C  IS DONE IN THE CALLING PROGRAM.  WHEN AN EVALUATION OF F IS
  27 C  NEEDED AT T, ROOT RETURNS WITH IFLAG NEGATIVE.  EVALUATE FT=F(T)
  28 C  AND CALL ROOT AGAIN.  DO NOT ALTER IFLAG.
  29 C
  30 C  WHEN THE COMPUTATION IS COMPLETE, ROOT RETURNS TO THE CALLING
  31 C  PROGRAM WITH IFLAG POSITIVE=
  32 C
  33 C     IFLAG=1  IF F(B)*F(C).LT.0 AND THE STOPPING CRITERION IS MET.
  34 C
  35 C          =2  IF A VALUE B IS FOUND SUCH THAT THE COMPUTED VALUE
  36 C              F(B) IS EXACTLY ZERO.  THE INTERVAL (B,C) MAY NOT
  37 C              SATISFY THE STOPPING CRITERION.
  38 C
  39 C          =3  IF DABS(F(B)) EXCEEDS THE INPUT VALUES DABS(F(B)),
  40 C              DABS(F(C)).  IN THIS CASE IT IS LIKELY THAT B IS CLOSE
  41 C              TO A POLE OF F.
  42 C
  43 C          =4  IF NO ODD ORDER ROOT WAS FOUND IN THE INTERVAL.  A
  44 C              LOCAL MINIMUM MAY HAVE BEEN OBTAINED.
  45 C
  46 C          =5  IF TOO MANY FUNCTION EVALUATIONS WERE MADE.
  47 C              (AS PROGRAMMED, 500 ARE ALLOWED.)
  48 C
  49 C  THIS CODE IS A MODIFICATION OF THE CODE ZEROIN WHICH IS COMPLETELY
  50 C  EXPLAINED AND DOCUMENTED IN THE TEXT  NUMERICAL COMPUTING:  AN
  51 C  INTRODUCTION,  BY L. F. SHAMPINE AND R. C. ALLEN.
  52 C
  53 C  CALLS  D1MACH .
  54 C
  55       DOUBLE PRECISION A,ABSERR,ACBS,ACMB,AE,B,C,CMB,D1MACH,FA,FB,
  56      1  FC,FT,FX,P,Q,RE,RELERR,T,TOL,U
  57       INTEGER IC,IFLAG,KOUNT
  58       SAVE
  59 C
  60       IF(IFLAG.GE.0) GO TO 100
  61       IFLAG=ABS(IFLAG)
  62       GO TO (200,300,400), IFLAG
  63   100 U=D1MACH(4)
  64       RE=MAX(RELERR,U)
  65       AE=MAX(ABSERR,0.0D0)
  66       IC=0
  67       ACBS=ABS(B-C)
  68       A=C
  69       T=A
  70       IFLAG=-1
  71       RETURN
  72   200 FA=FT
  73       T=B
  74       IFLAG=-2
  75       RETURN
  76   300 FB=FT
  77       FC=FA
  78       KOUNT=2
  79       FX=MAX(ABS(FB),ABS(FC))
  80     1 IF(ABS(FC).GE.ABS(FB))GO TO 2
  81 C
  82 C  INTERCHANGE B AND C SO THAT ABS(F(B)).LE.ABS(F(C)).
  83 C
  84       A=B
  85       FA=FB
  86       B=C
  87       FB=FC
  88       C=A
  89       FC=FA
  90     2 CMB=0.5*(C-B)
  91       ACMB=ABS(CMB)
  92       TOL=RE*ABS(B)+AE
  93 C
  94 C  TEST STOPPING CRITERION AND FUNCTION COUNT.
  95 C
  96       IF(ACMB.LE.TOL)GO TO 8
  97       IF(KOUNT.GE.500)GO TO 12
  98 C
  99 C  CALCULATE NEW ITERATE EXPLICITLY AS B+P/Q
 100 C  WHERE WE ARRANGE P.GE.0.  THE IMPLICIT
 101 C  FORM IS USED TO PREVENT OVERFLOW.
 102 C
 103       P=(B-A)*FB
 104       Q=FA-FB
 105       IF(P.GE.0.0)GO TO 3
 106       P=-P
 107       Q=-Q
 108 C
 109 C  UPDATE A, CHECK IF REDUCTION IN THE SIZE OF BRACKETING
 110 C  INTERVAL IS SATISFACTORY.  IF NOT BISECT UNTIL IT IS.
 111 C
 112     3 A=B
 113       FA=FB
 114       IC=IC+1
 115       IF(IC.LT.4)GO TO 4
 116       IF(8.0*ACMB.GE.ACBS)GO TO 6
 117       IC=0
 118       ACBS=ACMB
 119 C
 120 C  TEST FOR TOO SMALL A CHANGE.
 121 C
 122     4 IF(P.GT.ABS(Q)*TOL)GO TO 5
 123 C
 124 C  INCREMENT BY TOLERANCE
 125 C
 126       B=B+SIGN(TOL,CMB)
 127       GO TO 7
 128 C
 129 C  ROOT OUGHT TO BE BETWEEN B AND (C+B)/2
 130 C
 131     5 IF(P.GE.CMB*Q)GO TO 6
 132 C
 133 C  USE SECANT RULE.
 134 C
 135       B=B+P/Q
 136       GO TO 7
 137 C
 138 C  USE BISECTION.
 139 C
 140     6 B=0.5*(C+B)
 141 C
 142 C  HAVE COMPLETED COMPUTATION FOR NEW ITERATE B.
 143 C
 144     7 T=B
 145       IFLAG=-3
 146       RETURN
 147   400 FB=FT
 148       IF(FB.EQ.0.0)GO TO 9
 149       KOUNT=KOUNT+1
 150       IF(SIGN(1.0D0,FB).NE.SIGN(1.0D0,FC))GO TO 1
 151       C=A
 152       FC=FA
 153       GO TO 1
 154 C
 155 C FINISHED.  SET IFLAG.
 156 C
 157     8 IF(SIGN(1.0D0,FB).EQ.SIGN(1.0D0,FC))GO TO 11
 158       IF(ABS(FB).GT.FX)GO TO 10
 159       IFLAG=1
 160       RETURN
 161     9 IFLAG=2
 162       RETURN
 163    10 IFLAG=3
 164       RETURN
 165    11 IFLAG=4
 166       RETURN
 167    12 IFLAG=5
 168       RETURN
 169       END