droots.f

Go to the documentation of this file.

      SUBROUTINE droots (NG, HMIN, JFLAG, X0, X1, G0, G1, GX, X, JROOT,
     *                   imax, last, alpha, x2)
C
C***BEGIN PROLOGUE  DROOTS
C***REFER TO DDASRT
C***ROUTINES CALLED  DCOPY
C***DATE WRITTEN   821001   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C***END PROLOGUE  DROOTS
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      INTEGER ng, jflag, jroot, imax, last
      DOUBLE PRECISION hmin, x0, x1, g0, g1, gx, x, alpha, x2
      dimension g0(ng), g1(ng), gx(ng), jroot(ng)
C-----------------------------------------------------------------------
C THIS SUBROUTINE FINDS THE LEFTMOST ROOT OF A SET OF ARBITRARY
C FUNCTIONS GI(X) (I = 1,...,NG) IN AN INTERVAL (X0,X1).  ONLY ROOTS
C OF ODD MULTIPLICITY (I.E. CHANGES OF SIGN OF THE GI) ARE FOUND.
C HERE THE SIGN OF X1 - X0 IS ARBITRARY, BUT IS CONSTANT FOR A GIVEN
C PROBLEM, AND -LEFTMOST- MEANS NEAREST TO X0.
C THE VALUES OF THE VECTOR-VALUED FUNCTION G(X) = (GI, I=1...NG)
C ARE COMMUNICATED THROUGH THE CALL SEQUENCE OF DROOTS.
C THE METHOD USED IS THE ILLINOIS ALGORITHM.
C
C REFERENCE..
C KATHIE L. HIEBERT AND LAWRENCE F. SHAMPINE, IMPLICITLY DEFINED
C OUTPUT POINTS FOR SOLUTIONS OF ODE-S, SANDIA REPORT SAND80-0180,
C FEBRUARY, 1980.
C
C DESCRIPTION OF PARAMETERS.
C
C NG     = NUMBER OF FUNCTIONS GI, OR THE NUMBER OF COMPONENTS OF
C          THE VECTOR VALUED FUNCTION G(X).  INPUT ONLY.
C
C HMIN   = RESOLUTION PARAMETER IN X.  INPUT ONLY.  WHEN A ROOT IS
C          FOUND, IT IS LOCATED ONLY TO WITHIN AN ERROR OF HMIN IN X.
C          TYPICALLY, HMIN SHOULD BE SET TO SOMETHING ON THE ORDER OF
C               100 * UROUND * MAX(ABS(X0),ABS(X1)),
C          WHERE UROUND IS THE UNIT ROUNDOFF OF THE MACHINE.
C
C JFLAG  = INTEGER FLAG FOR INPUT AND OUTPUT COMMUNICATION.
C
C          ON INPUT, SET JFLAG = 0 ON THE FIRST CALL FOR THE PROBLEM,
C          AND LEAVE IT UNCHANGED UNTIL THE PROBLEM IS COMPLETED.
C          (THE PROBLEM IS COMPLETED WHEN JFLAG .GE. 2 ON RETURN.)
C
C          ON OUTPUT, JFLAG HAS THE FOLLOWING VALUES AND MEANINGS..
C          JFLAG = 1 MEANS DROOTS NEEDS A VALUE OF G(X).  SET GX = G(X)
C                    AND CALL DROOTS AGAIN.
C          JFLAG = 2 MEANS A ROOT HAS BEEN FOUND.  THE ROOT IS
C                    AT X, AND GX CONTAINS G(X).  (ACTUALLY, X IS THE
C                    RIGHTMOST APPROXIMATION TO THE ROOT ON AN INTERVAL
C                    (X0,X1) OF SIZE HMIN OR LESS.)
C          JFLAG = 3 MEANS X = X1 IS A ROOT, WITH ONE OR MORE OF THE GI
C                    BEING ZERO AT X1 AND NO SIGN CHANGES IN (X0,X1).
C                    GX CONTAINS G(X) ON OUTPUT.
C          JFLAG = 4 MEANS NO ROOTS (OF ODD MULTIPLICITY) WERE
C                    FOUND IN (X0,X1) (NO SIGN CHANGES).
C
C X0,X1  = ENDPOINTS OF THE INTERVAL WHERE ROOTS ARE SOUGHT.
C          X1 AND X0 ARE INPUT WHEN JFLAG = 0 (FIRST CALL), AND
C          MUST BE LEFT UNCHANGED BETWEEN CALLS UNTIL THE PROBLEM IS
C          COMPLETED.  X0 AND X1 MUST BE DISTINCT, BUT X1 - X0 MAY BE
C          OF EITHER SIGN.  HOWEVER, THE NOTION OF -LEFT- AND -RIGHT-
C          WILL BE USED TO MEAN NEARER TO X0 OR X1, RESPECTIVELY.
C          WHEN JFLAG .GE. 2 ON RETURN, X0 AND X1 ARE OUTPUT, AND
C          ARE THE ENDPOINTS OF THE RELEVANT INTERVAL.
C
C G0,G1  = ARRAYS OF LENGTH NG CONTAINING THE VECTORS G(X0) AND G(X1),
C          RESPECTIVELY.  WHEN JFLAG = 0, G0 AND G1 ARE INPUT AND
C          NONE OF THE G0(I) SHOULD BE BE ZERO.
C          WHEN JFLAG .GE. 2 ON RETURN, G0 AND G1 ARE OUTPUT.
C
C GX     = ARRAY OF LENGTH NG CONTAINING G(X).  GX IS INPUT
C          WHEN JFLAG = 1, AND OUTPUT WHEN JFLAG .GE. 2.
C
C X      = INDEPENDENT VARIABLE VALUE.  OUTPUT ONLY.
C          WHEN JFLAG = 1 ON OUTPUT, X IS THE POINT AT WHICH G(X)
C          IS TO BE EVALUATED AND LOADED INTO GX.
C          WHEN JFLAG = 2 OR 3, X IS THE ROOT.
C          WHEN JFLAG = 4, X IS THE RIGHT ENDPOINT OF THE INTERVAL, X1.
C
C JROOT  = INTEGER ARRAY OF LENGTH NG.  OUTPUT ONLY.
C          WHEN JFLAG = 2 OR 3, JROOT INDICATES WHICH COMPONENTS
C          OF G(X) HAVE A ROOT AT X.  JROOT(I) IS 1 IF THE I-TH
C          COMPONENT HAS A ROOT, AND JROOT(I) = 0 OTHERWISE.
C
C IMAX, LAST, ALPHA, X2 =
C          BOOKKEEPING VARIABLES WHICH MUST BE SAVED FROM CALL
C          TO CALL.  THEY ARE SAVED INSIDE THE CALLING ROUTINE,
C          BUT THEY ARE USED ONLY WITHIN THIS ROUTINE.
C-----------------------------------------------------------------------
      INTEGER i, imxold, nxlast
      DOUBLE PRECISION t2, tmax, zero
      LOGICAL zroot, sgnchg, xroot
      DATA zero/0.0d0/
C
      IF (jflag .EQ. 1) go to 200
C JFLAG .NE. 1.  CHECK FOR CHANGE IN SIGN OF G OR ZERO AT X1. ----------
      imax = 0
      tmax = zero
      zroot = .false.
      DO 120 i = 1,ng
        IF (dabs(g1(i)) .GT. zero) go to 110
        zroot = .true.
        go to 120
C AT THIS POINT, G0(I) HAS BEEN CHECKED AND CANNOT BE ZERO. ------------
 110    IF (dsign(1.0d0,g0(i)) .EQ. dsign(1.0d0,g1(i))) go to 120
          t2 = dabs(g1(i)/(g1(i)-g0(i)))
          IF (t2 .LE. tmax) go to 120
            tmax = t2
            imax = i
 120    CONTINUE
      IF (imax .GT. 0) go to 130
      sgnchg = .false.
      go to 140
 130  sgnchg = .true.
 140  IF (.NOT. sgnchg) go to 400
C THERE IS A SIGN CHANGE.  FIND THE FIRST ROOT IN THE INTERVAL. --------
      xroot = .false.
      nxlast = 0
      last = 1
C
C REPEAT UNTIL THE FIRST ROOT IN THE INTERVAL IS FOUND.  LOOP POINT. ---
 150  CONTINUE
      IF (xroot) go to 300
      IF (nxlast .EQ. last) go to 160
      alpha = 1.0d0
      go to 180
 160  IF (last .EQ. 0) go to 170
      alpha = 0.5d0*alpha
      go to 180
 170  alpha = 2.0d0*alpha
 180  x2 = x1 - (x1-x0)*g1(imax)/(g1(imax) - alpha*g0(imax))
      IF ((dabs(x2-x0) .LT. hmin) .AND.
     1   (dabs(x1-x0) .GT. 10.0d0*hmin)) x2 = x0 + 0.1d0*(x1-x0)
      jflag = 1
      x = x2
C RETURN TO THE CALLING ROUTINE TO GET A VALUE OF GX = G(X). -----------
      RETURN
C CHECK TO SEE IN WHICH INTERVAL G CHANGES SIGN. -----------------------
 200  imxold = imax
      imax = 0
      tmax = zero
      zroot = .false.
      DO 220 i = 1,ng
        IF (dabs(gx(i)) .GT. zero) go to 210
        zroot = .true.
        go to 220
C NEITHER G0(I) NOR GX(I) CAN BE ZERO AT THIS POINT. -------------------
 210    IF (dsign(1.0d0,g0(i)) .EQ. dsign(1.0d0,gx(i))) go to 220
          t2 = dabs(gx(i)/(gx(i) - g0(i)))
          IF (t2 .LE. tmax) go to 220
            tmax = t2
            imax = i
 220    CONTINUE
      IF (imax .GT. 0) go to 230
      sgnchg = .false.
      imax = imxold
      go to 240
 230  sgnchg = .true.
 240  nxlast = last
      IF (.NOT. sgnchg) go to 250
C SIGN CHANGE BETWEEN X0 AND X2, SO REPLACE X1 WITH X2. ----------------
      x1 = x2
      CALL dcopy(ng, gx, 1, g1, 1)
      last = 1
      xroot = .false.
      go to 270
 250  IF (.NOT. zroot) go to 260
C ZERO VALUE AT X2 AND NO SIGN CHANGE IN (X0,X2), SO X2 IS A ROOT. -----
      x1 = x2
      CALL dcopy(ng, gx, 1, g1, 1)
      xroot = .true.
      go to 270
C NO SIGN CHANGE BETWEEN X0 AND X2.  REPLACE X0 WITH X2. ---------------
 260  CONTINUE
      CALL dcopy(ng, gx, 1, g0, 1)
      x0 = x2
      last = 0
      xroot = .false.
 270  IF (dabs(x1-x0) .LE. hmin) xroot = .true.
      go to 150
C
C RETURN WITH X1 AS THE ROOT.  SET JROOT.  SET X = X1 AND GX = G1. -----
 300  jflag = 2
      x = x1
      CALL dcopy(ng, g1, 1, gx, 1)
      DO 320 i = 1,ng
        jroot(i) = 0
        IF (dabs(g1(i)) .GT. zero) go to 310
          jroot(i) = 1
          go to 320
 310    IF (dsign(1.0d0,g0(i)) .NE. dsign(1.0d0,g1(i))) jroot(i) = 1
 320    CONTINUE
      RETURN
C
C NO SIGN CHANGE IN THE INTERVAL.  CHECK FOR ZERO AT RIGHT ENDPOINT. ---
 400  IF (.NOT. zroot) go to 420
C
C ZERO VALUE AT X1 AND NO SIGN CHANGE IN (X0,X1).  RETURN JFLAG = 3. ---
      x = x1
      CALL dcopy(ng, g1, 1, gx, 1)
      DO 410 i = 1,ng
        jroot(i) = 0
        IF (dabs(g1(i)) .LE. zero) jroot(i) = 1
 410  CONTINUE
      jflag = 3
      RETURN
C
C NO SIGN CHANGES IN THIS INTERVAL.  SET X = X1, RETURN JFLAG = 4. -----
 420  CALL dcopy(ng, g1, 1, gx, 1)
      x = x1
      jflag = 4
      RETURN
C---------------------- END OF SUBROUTINE DROOTS -----------------------
      END