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
1.7.1
