droots.f

Go to the documentation of this file.
00001       SUBROUTINE DROOTS (NG, HMIN, JFLAG, X0, X1, G0, G1, GX, X, JROOT,
00002      *                   IMAX, LAST, ALPHA, X2)
00003 C
00004 C***BEGIN PROLOGUE  DROOTS
00005 C***REFER TO DDASRT
00006 C***ROUTINES CALLED  DCOPY
00007 C***DATE WRITTEN   821001   (YYMMDD)
00008 C***REVISION DATE  900926   (YYMMDD)
00009 C***END PROLOGUE  DROOTS
00010 C
00011       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
00012       INTEGER NG, JFLAG, JROOT, IMAX, LAST
00013       DOUBLE PRECISION HMIN, X0, X1, G0, G1, GX, X, ALPHA, X2
00014       DIMENSION G0(NG), G1(NG), GX(NG), JROOT(NG)
00015 C-----------------------------------------------------------------------
00016 C THIS SUBROUTINE FINDS THE LEFTMOST ROOT OF A SET OF ARBITRARY
00017 C FUNCTIONS GI(X) (I = 1,...,NG) IN AN INTERVAL (X0,X1).  ONLY ROOTS
00018 C OF ODD MULTIPLICITY (I.E. CHANGES OF SIGN OF THE GI) ARE FOUND.
00019 C HERE THE SIGN OF X1 - X0 IS ARBITRARY, BUT IS CONSTANT FOR A GIVEN
00020 C PROBLEM, AND -LEFTMOST- MEANS NEAREST TO X0.
00021 C THE VALUES OF THE VECTOR-VALUED FUNCTION G(X) = (GI, I=1...NG)
00022 C ARE COMMUNICATED THROUGH THE CALL SEQUENCE OF DROOTS.
00023 C THE METHOD USED IS THE ILLINOIS ALGORITHM.
00024 C
00025 C REFERENCE..
00026 C KATHIE L. HIEBERT AND LAWRENCE F. SHAMPINE, IMPLICITLY DEFINED
00027 C OUTPUT POINTS FOR SOLUTIONS OF ODE-S, SANDIA REPORT SAND80-0180,
00028 C FEBRUARY, 1980.
00029 C
00030 C DESCRIPTION OF PARAMETERS.
00031 C
00032 C NG     = NUMBER OF FUNCTIONS GI, OR THE NUMBER OF COMPONENTS OF
00033 C          THE VECTOR VALUED FUNCTION G(X).  INPUT ONLY.
00034 C
00035 C HMIN   = RESOLUTION PARAMETER IN X.  INPUT ONLY.  WHEN A ROOT IS
00036 C          FOUND, IT IS LOCATED ONLY TO WITHIN AN ERROR OF HMIN IN X.
00037 C          TYPICALLY, HMIN SHOULD BE SET TO SOMETHING ON THE ORDER OF
00038 C               100 * UROUND * MAX(ABS(X0),ABS(X1)),
00039 C          WHERE UROUND IS THE UNIT ROUNDOFF OF THE MACHINE.
00040 C
00041 C JFLAG  = INTEGER FLAG FOR INPUT AND OUTPUT COMMUNICATION.
00042 C
00043 C          ON INPUT, SET JFLAG = 0 ON THE FIRST CALL FOR THE PROBLEM,
00044 C          AND LEAVE IT UNCHANGED UNTIL THE PROBLEM IS COMPLETED.
00045 C          (THE PROBLEM IS COMPLETED WHEN JFLAG .GE. 2 ON RETURN.)
00046 C
00047 C          ON OUTPUT, JFLAG HAS THE FOLLOWING VALUES AND MEANINGS..
00048 C          JFLAG = 1 MEANS DROOTS NEEDS A VALUE OF G(X).  SET GX = G(X)
00049 C                    AND CALL DROOTS AGAIN.
00050 C          JFLAG = 2 MEANS A ROOT HAS BEEN FOUND.  THE ROOT IS
00051 C                    AT X, AND GX CONTAINS G(X).  (ACTUALLY, X IS THE
00052 C                    RIGHTMOST APPROXIMATION TO THE ROOT ON AN INTERVAL
00053 C                    (X0,X1) OF SIZE HMIN OR LESS.)
00054 C          JFLAG = 3 MEANS X = X1 IS A ROOT, WITH ONE OR MORE OF THE GI
00055 C                    BEING ZERO AT X1 AND NO SIGN CHANGES IN (X0,X1).
00056 C                    GX CONTAINS G(X) ON OUTPUT.
00057 C          JFLAG = 4 MEANS NO ROOTS (OF ODD MULTIPLICITY) WERE
00058 C                    FOUND IN (X0,X1) (NO SIGN CHANGES).
00059 C
00060 C X0,X1  = ENDPOINTS OF THE INTERVAL WHERE ROOTS ARE SOUGHT.
00061 C          X1 AND X0 ARE INPUT WHEN JFLAG = 0 (FIRST CALL), AND
00062 C          MUST BE LEFT UNCHANGED BETWEEN CALLS UNTIL THE PROBLEM IS
00063 C          COMPLETED.  X0 AND X1 MUST BE DISTINCT, BUT X1 - X0 MAY BE
00064 C          OF EITHER SIGN.  HOWEVER, THE NOTION OF -LEFT- AND -RIGHT-
00065 C          WILL BE USED TO MEAN NEARER TO X0 OR X1, RESPECTIVELY.
00066 C          WHEN JFLAG .GE. 2 ON RETURN, X0 AND X1 ARE OUTPUT, AND
00067 C          ARE THE ENDPOINTS OF THE RELEVANT INTERVAL.
00068 C
00069 C G0,G1  = ARRAYS OF LENGTH NG CONTAINING THE VECTORS G(X0) AND G(X1),
00070 C          RESPECTIVELY.  WHEN JFLAG = 0, G0 AND G1 ARE INPUT AND
00071 C          NONE OF THE G0(I) SHOULD BE BE ZERO.
00072 C          WHEN JFLAG .GE. 2 ON RETURN, G0 AND G1 ARE OUTPUT.
00073 C
00074 C GX     = ARRAY OF LENGTH NG CONTAINING G(X).  GX IS INPUT
00075 C          WHEN JFLAG = 1, AND OUTPUT WHEN JFLAG .GE. 2.
00076 C
00077 C X      = INDEPENDENT VARIABLE VALUE.  OUTPUT ONLY.
00078 C          WHEN JFLAG = 1 ON OUTPUT, X IS THE POINT AT WHICH G(X)
00079 C          IS TO BE EVALUATED AND LOADED INTO GX.
00080 C          WHEN JFLAG = 2 OR 3, X IS THE ROOT.
00081 C          WHEN JFLAG = 4, X IS THE RIGHT ENDPOINT OF THE INTERVAL, X1.
00082 C
00083 C JROOT  = INTEGER ARRAY OF LENGTH NG.  OUTPUT ONLY.
00084 C          WHEN JFLAG = 2 OR 3, JROOT INDICATES WHICH COMPONENTS
00085 C          OF G(X) HAVE A ROOT AT X.  JROOT(I) IS 1 IF THE I-TH
00086 C          COMPONENT HAS A ROOT, AND JROOT(I) = 0 OTHERWISE.
00087 C
00088 C IMAX, LAST, ALPHA, X2 =
00089 C          BOOKKEEPING VARIABLES WHICH MUST BE SAVED FROM CALL
00090 C          TO CALL.  THEY ARE SAVED INSIDE THE CALLING ROUTINE,
00091 C          BUT THEY ARE USED ONLY WITHIN THIS ROUTINE.
00092 C-----------------------------------------------------------------------
00093       INTEGER I, IMXOLD, NXLAST
00094       DOUBLE PRECISION T2, TMAX, ZERO
00095       LOGICAL ZROOT, SGNCHG, XROOT
00096       DATA ZERO/0.0D0/
00097 C
00098       IF (JFLAG .EQ. 1) GO TO 200
00099 C JFLAG .NE. 1.  CHECK FOR CHANGE IN SIGN OF G OR ZERO AT X1. ----------
00100       IMAX = 0
00101       TMAX = ZERO
00102       ZROOT = .FALSE.
00103       DO 120 I = 1,NG
00104         IF (DABS(G1(I)) .GT. ZERO) GO TO 110
00105         ZROOT = .TRUE.
00106         GO TO 120
00107 C AT THIS POINT, G0(I) HAS BEEN CHECKED AND CANNOT BE ZERO. ------------
00108  110    IF (DSIGN(1.0D0,G0(I)) .EQ. DSIGN(1.0D0,G1(I))) GO TO 120
00109           T2 = DABS(G1(I)/(G1(I)-G0(I)))
00110           IF (T2 .LE. TMAX) GO TO 120
00111             TMAX = T2
00112             IMAX = I
00113  120    CONTINUE
00114       IF (IMAX .GT. 0) GO TO 130
00115       SGNCHG = .FALSE.
00116       GO TO 140
00117  130  SGNCHG = .TRUE.
00118  140  IF (.NOT. SGNCHG) GO TO 400
00119 C THERE IS A SIGN CHANGE.  FIND THE FIRST ROOT IN THE INTERVAL. --------
00120       XROOT = .FALSE.
00121       NXLAST = 0
00122       LAST = 1
00123 C
00124 C REPEAT UNTIL THE FIRST ROOT IN THE INTERVAL IS FOUND.  LOOP POINT. ---
00125  150  CONTINUE
00126       IF (XROOT) GO TO 300
00127       IF (NXLAST .EQ. LAST) GO TO 160
00128       ALPHA = 1.0D0
00129       GO TO 180
00130  160  IF (LAST .EQ. 0) GO TO 170
00131       ALPHA = 0.5D0*ALPHA
00132       GO TO 180
00133  170  ALPHA = 2.0D0*ALPHA
00134  180  X2 = X1 - (X1-X0)*G1(IMAX)/(G1(IMAX) - ALPHA*G0(IMAX))
00135       IF ((DABS(X2-X0) .LT. HMIN) .AND.
00136      1   (DABS(X1-X0) .GT. 10.0D0*HMIN)) X2 = X0 + 0.1D0*(X1-X0)
00137       JFLAG = 1
00138       X = X2
00139 C RETURN TO THE CALLING ROUTINE TO GET A VALUE OF GX = G(X). -----------
00140       RETURN
00141 C CHECK TO SEE IN WHICH INTERVAL G CHANGES SIGN. -----------------------
00142  200  IMXOLD = IMAX
00143       IMAX = 0
00144       TMAX = ZERO
00145       ZROOT = .FALSE.
00146       DO 220 I = 1,NG
00147         IF (DABS(GX(I)) .GT. ZERO) GO TO 210
00148         ZROOT = .TRUE.
00149         GO TO 220
00150 C NEITHER G0(I) NOR GX(I) CAN BE ZERO AT THIS POINT. -------------------
00151  210    IF (DSIGN(1.0D0,G0(I)) .EQ. DSIGN(1.0D0,GX(I))) GO TO 220
00152           T2 = DABS(GX(I)/(GX(I) - G0(I)))
00153           IF (T2 .LE. TMAX) GO TO 220
00154             TMAX = T2
00155             IMAX = I
00156  220    CONTINUE
00157       IF (IMAX .GT. 0) GO TO 230
00158       SGNCHG = .FALSE.
00159       IMAX = IMXOLD
00160       GO TO 240
00161  230  SGNCHG = .TRUE.
00162  240  NXLAST = LAST
00163       IF (.NOT. SGNCHG) GO TO 250
00164 C SIGN CHANGE BETWEEN X0 AND X2, SO REPLACE X1 WITH X2. ----------------
00165       X1 = X2
00166       CALL DCOPY (NG, GX, 1, G1, 1)
00167       LAST = 1
00168       XROOT = .FALSE.
00169       GO TO 270
00170  250  IF (.NOT. ZROOT) GO TO 260
00171 C ZERO VALUE AT X2 AND NO SIGN CHANGE IN (X0,X2), SO X2 IS A ROOT. -----
00172       X1 = X2
00173       CALL DCOPY (NG, GX, 1, G1, 1)
00174       XROOT = .TRUE.
00175       GO TO 270
00176 C NO SIGN CHANGE BETWEEN X0 AND X2.  REPLACE X0 WITH X2. ---------------
00177  260  CONTINUE
00178       CALL DCOPY (NG, GX, 1, G0, 1)
00179       X0 = X2
00180       LAST = 0
00181       XROOT = .FALSE.
00182  270  IF (DABS(X1-X0) .LE. HMIN) XROOT = .TRUE.
00183       GO TO 150
00184 C
00185 C RETURN WITH X1 AS THE ROOT.  SET JROOT.  SET X = X1 AND GX = G1. -----
00186  300  JFLAG = 2
00187       X = X1
00188       CALL DCOPY (NG, G1, 1, GX, 1)
00189       DO 320 I = 1,NG
00190         JROOT(I) = 0
00191         IF (DABS(G1(I)) .GT. ZERO) GO TO 310
00192           JROOT(I) = 1
00193           GO TO 320
00194  310    IF (DSIGN(1.0D0,G0(I)) .NE. DSIGN(1.0D0,G1(I))) JROOT(I) = 1
00195  320    CONTINUE
00196       RETURN
00197 C
00198 C NO SIGN CHANGE IN THE INTERVAL.  CHECK FOR ZERO AT RIGHT ENDPOINT. ---
00199  400  IF (.NOT. ZROOT) GO TO 420
00200 C
00201 C ZERO VALUE AT X1 AND NO SIGN CHANGE IN (X0,X1).  RETURN JFLAG = 3. ---
00202       X = X1
00203       CALL DCOPY (NG, G1, 1, GX, 1)
00204       DO 410 I = 1,NG
00205         JROOT(I) = 0
00206         IF (DABS(G1(I)) .LE. ZERO) JROOT (I) = 1
00207  410  CONTINUE
00208       JFLAG = 3
00209       RETURN
00210 C
00211 C NO SIGN CHANGES IN THIS INTERVAL.  SET X = X1, RETURN JFLAG = 4. -----
00212  420  CALL DCOPY (NG, G1, 1, GX, 1)
00213       X = X1
00214       JFLAG = 4
00215       RETURN
00216 C---------------------- END OF SUBROUTINE DROOTS -----------------------
00217       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines