share/odepack/fortran/droots.f

   1 *DECK DROOTS
   2       SUBROUTINE DROOTS (NG, HMIN, JFLAG, X0, X1, G0, G1, GX, X, JROOT)
   3       INTEGER NG, JFLAG, JROOT
   4       DOUBLE PRECISION HMIN, X0, X1, G0, G1, GX, X
   5       DIMENSION G0(NG), G1(NG), GX(NG), JROOT(NG)
   6       INTEGER IOWND3, IMAX, LAST, IDUM3
   7       DOUBLE PRECISION ALPHA, X2, RDUM3
   8       COMMON /DLSR01/ ALPHA, X2, RDUM3(3),
   9      1   IOWND3(3), IMAX, LAST, IDUM3(4)
  10 C-----------------------------------------------------------------------
  11 C This subroutine finds the leftmost root of a set of arbitrary
  12 C functions gi(x) (i = 1,...,NG) in an interval (X0,X1).  Only roots
  13 C of odd multiplicity (i.e. changes of sign of the gi) are found.
  14 C Here the sign of X1 - X0 is arbitrary, but is constant for a given
  15 C problem, and -leftmost- means nearest to X0.
  16 C The values of the vector-valued function g(x) = (gi, i=1...NG)
  17 C are communicated through the call sequence of DROOTS.
  18 C The method used is the Illinois algorithm.
  19 C
  20 C Reference:
  21 C Kathie L. Hiebert and Lawrence F. Shampine, Implicitly Defined
  22 C Output Points for Solutions of ODEs, Sandia Report SAND80-0180,
  23 C February 1980.
  24 C
  25 C Description of parameters.
  26 C
  27 C NG     = number of functions gi, or the number of components of
  28 C          the vector valued function g(x).  Input only.
  29 C
  30 C HMIN   = resolution parameter in X.  Input only.  When a root is
  31 C          found, it is located only to within an error of HMIN in X.
  32 C          Typically, HMIN should be set to something on the order of
  33 C               100 * UROUND * MAX(ABS(X0),ABS(X1)),
  34 C          where UROUND is the unit roundoff of the machine.
  35 C
  36 C JFLAG  = integer flag for input and output communication.
  37 C
  38 C          On input, set JFLAG = 0 on the first call for the problem,
  39 C          and leave it unchanged until the problem is completed.
  40 C          (The problem is completed when JFLAG .ge. 2 on return.)
  41 C
  42 C          On output, JFLAG has the following values and meanings:
  43 C          JFLAG = 1 means DROOTS needs a value of g(x).  Set GX = g(X)
  44 C                    and call DROOTS again.
  45 C          JFLAG = 2 means a root has been found.  The root is
  46 C                    at X, and GX contains g(X).  (Actually, X is the
  47 C                    rightmost approximation to the root on an interval
  48 C                    (X0,X1) of size HMIN or less.)
  49 C          JFLAG = 3 means X = X1 is a root, with one or more of the gi
  50 C                    being zero at X1 and no sign changes in (X0,X1).
  51 C                    GX contains g(X) on output.
  52 C          JFLAG = 4 means no roots (of odd multiplicity) were
  53 C                    found in (X0,X1) (no sign changes).
  54 C
  55 C X0,X1  = endpoints of the interval where roots are sought.
  56 C          X1 and X0 are input when JFLAG = 0 (first call), and
  57 C          must be left unchanged between calls until the problem is
  58 C          completed.  X0 and X1 must be distinct, but X1 - X0 may be
  59 C          of either sign.  However, the notion of -left- and -right-
  60 C          will be used to mean nearer to X0 or X1, respectively.
  61 C          When JFLAG .ge. 2 on return, X0 and X1 are output, and
  62 C          are the endpoints of the relevant interval.
  63 C
  64 C G0,G1  = arrays of length NG containing the vectors g(X0) and g(X1),
  65 C          respectively.  When JFLAG = 0, G0 and G1 are input and
  66 C          none of the G0(i) should be zero.
  67 C          When JFLAG .ge. 2 on return, G0 and G1 are output.
  68 C
  69 C GX     = array of length NG containing g(X).  GX is input
  70 C          when JFLAG = 1, and output when JFLAG .ge. 2.
  71 C
  72 C X      = independent variable value.  Output only.
  73 C          When JFLAG = 1 on output, X is the point at which g(x)
  74 C          is to be evaluated and loaded into GX.
  75 C          When JFLAG = 2 or 3, X is the root.
  76 C          When JFLAG = 4, X is the right endpoint of the interval, X1.
  77 C
  78 C JROOT  = integer array of length NG.  Output only.
  79 C          When JFLAG = 2 or 3, JROOT indicates which components
  80 C          of g(x) have a root at X.  JROOT(i) is 1 if the i-th
  81 C          component has a root, and JROOT(i) = 0 otherwise.
  82 C-----------------------------------------------------------------------
  83       INTEGER I, IMXOLD, NXLAST
  84       DOUBLE PRECISION T2, TMAX, ZERO
  85       LOGICAL ZROOT, SGNCHG, XROOT
  86       SAVE ZERO
  87       DATA ZERO/0.0D0/
  88 C
  89       IF (JFLAG .EQ. 1) GO TO 200
  90 C JFLAG .ne. 1.  Check for change in sign of g or zero at X1. ----------
  91       IMAX = 0
  92       TMAX = ZERO
  93       ZROOT = .FALSE.
  94       DO 120 I = 1,NG
  95         IF (ABS(G1(I)) .GT. ZERO) GO TO 110
  96         ZROOT = .TRUE.
  97         GO TO 120
  98 C At this point, G0(i) has been checked and cannot be zero. ------------
  99  110    IF (SIGN(1.0D0,G0(I)) .EQ. SIGN(1.0D0,G1(I))) GO TO 120
 100           T2 = ABS(G1(I)/(G1(I)-G0(I)))
 101           IF (T2 .LE. TMAX) GO TO 120
 102             TMAX = T2
 103             IMAX = I
 104  120    CONTINUE
 105       IF (IMAX .GT. 0) GO TO 130
 106       SGNCHG = .FALSE.
 107       GO TO 140
 108  130  SGNCHG = .TRUE.
 109  140  IF (.NOT. SGNCHG) GO TO 400
 110 C There is a sign change.  Find the first root in the interval. --------
 111       XROOT = .FALSE.
 112       NXLAST = 0
 113       LAST = 1
 114 C
 115 C Repeat until the first root in the interval is found.  Loop point. ---
 116  150  CONTINUE
 117       IF (XROOT) GO TO 300
 118       IF (NXLAST .EQ. LAST) GO TO 160
 119       ALPHA = 1.0D0
 120       GO TO 180
 121  160  IF (LAST .EQ. 0) GO TO 170
 122       ALPHA = 0.5D0*ALPHA
 123       GO TO 180
 124  170  ALPHA = 2.0D0*ALPHA
 125  180  X2 = X1 - (X1-X0)*G1(IMAX)/(G1(IMAX) - ALPHA*G0(IMAX))
 126       IF ((ABS(X2-X0) .LT. HMIN) .AND.
 127      1   (ABS(X1-X0) .GT. 10.0D0*HMIN)) X2 = X0 + 0.1D0*(X1-X0)
 128       JFLAG = 1
 129       X = X2
 130 C Return to the calling routine to get a value of GX = g(X). -----------
 131       RETURN
 132 C Check to see in which interval g changes sign. -----------------------
 133  200  IMXOLD = IMAX
 134       IMAX = 0
 135       TMAX = ZERO
 136       ZROOT = .FALSE.
 137       DO 220 I = 1,NG
 138         IF (ABS(GX(I)) .GT. ZERO) GO TO 210
 139         ZROOT = .TRUE.
 140         GO TO 220
 141 C Neither G0(i) nor GX(i) can be zero at this point. -------------------
 142  210    IF (SIGN(1.0D0,G0(I)) .EQ. SIGN(1.0D0,GX(I))) GO TO 220
 143           T2 = ABS(GX(I)/(GX(I) - G0(I)))
 144           IF (T2 .LE. TMAX) GO TO 220
 145             TMAX = T2
 146             IMAX = I
 147  220    CONTINUE
 148       IF (IMAX .GT. 0) GO TO 230
 149       SGNCHG = .FALSE.
 150       IMAX = IMXOLD
 151       GO TO 240
 152  230  SGNCHG = .TRUE.
 153  240  NXLAST = LAST
 154       IF (.NOT. SGNCHG) GO TO 250
 155 C Sign change between X0 and X2, so replace X1 with X2. ----------------
 156       X1 = X2
 157       CALL DCOPY (NG, GX, 1, G1, 1)
 158       LAST = 1
 159       XROOT = .FALSE.
 160       GO TO 270
 161  250  IF (.NOT. ZROOT) GO TO 260
 162 C Zero value at X2 and no sign change in (X0,X2), so X2 is a root. -----
 163       X1 = X2
 164       CALL DCOPY (NG, GX, 1, G1, 1)
 165       XROOT = .TRUE.
 166       GO TO 270
 167 C No sign change between X0 and X2.  Replace X0 with X2. ---------------
 168  260  CONTINUE
 169       CALL DCOPY (NG, GX, 1, G0, 1)
 170       X0 = X2
 171       LAST = 0
 172       XROOT = .FALSE.
 173  270  IF (ABS(X1-X0) .LE. HMIN) XROOT = .TRUE.
 174       GO TO 150
 175 C
 176 C Return with X1 as the root.  Set JROOT.  Set X = X1 and GX = G1. -----
 177  300  JFLAG = 2
 178       X = X1
 179       CALL DCOPY (NG, G1, 1, GX, 1)
 180       DO 320 I = 1,NG
 181         JROOT(I) = 0
 182         IF (ABS(G1(I)) .GT. ZERO) GO TO 310
 183           JROOT(I) = 1
 184           GO TO 320
 185  310    IF (SIGN(1.0D0,G0(I)) .NE. SIGN(1.0D0,G1(I))) JROOT(I) = 1
 186  320    CONTINUE
 187       RETURN
 188 C
 189 C No sign change in the interval.  Check for zero at right endpoint. ---
 190  400  IF (.NOT. ZROOT) GO TO 420
 191 C
 192 C Zero value at X1 and no sign change in (X0,X1).  Return JFLAG = 3. ---
 193       X = X1
 194       CALL DCOPY (NG, G1, 1, GX, 1)
 195       DO 410 I = 1,NG
 196         JROOT(I) = 0
 197         IF (ABS(G1(I)) .LE. ZERO) JROOT (I) = 1
 198  410  CONTINUE
 199       JFLAG = 3
 200       RETURN
 201 C
 202 C No sign changes in this interval.  Set X = X1, return JFLAG = 4. -----
 203  420  CALL DCOPY (NG, G1, 1, GX, 1)
 204       X = X1
 205       JFLAG = 4
 206       RETURN
 207 C----------------------- End of Subroutine DROOTS ----------------------
 208       END