1 SUBROUTINE DROOTS (NG, HMIN, JFLAG, X0, X1, G0, G1, GX, X, JROOT, 2 * IMAX, LAST, ALPHA, X2) 3C 4C***BEGIN PROLOGUE DROOTS 5C***REFER TO DDASRT 6C***ROUTINES CALLED DCOPY 7C***DATE WRITTEN 821001 (YYMMDD) 8C***REVISION DATE 900926 (YYMMDD) 9C***END PROLOGUE DROOTS 10C 11 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 12 INTEGER NG, JFLAG, JROOT, IMAX, LAST 13 DOUBLE PRECISION HMIN, X0, X1, G0, G1, GX, X, ALPHA, X2 14 DIMENSION G0(NG), G1(NG), GX(NG), JROOT(NG) 15C----------------------------------------------------------------------- 16C THIS SUBROUTINE FINDS THE LEFTMOST ROOT OF A SET OF ARBITRARY 17C FUNCTIONS GI(X) (I = 1,...,NG) IN AN INTERVAL (X0,X1). ONLY ROOTS 18C OF ODD MULTIPLICITY (I.E. CHANGES OF SIGN OF THE GI) ARE FOUND. 19C HERE THE SIGN OF X1 - X0 IS ARBITRARY, BUT IS CONSTANT FOR A GIVEN 20C PROBLEM, AND -LEFTMOST- MEANS NEAREST TO X0. 21C THE VALUES OF THE VECTOR-VALUED FUNCTION G(X) = (GI, I=1...NG) 22C ARE COMMUNICATED THROUGH THE CALL SEQUENCE OF DROOTS. 23C THE METHOD USED IS THE ILLINOIS ALGORITHM. 24C 25C REFERENCE.. 26C KATHIE L. HIEBERT AND LAWRENCE F. SHAMPINE, IMPLICITLY DEFINED 27C OUTPUT POINTS FOR SOLUTIONS OF ODE-S, SANDIA REPORT SAND80-0180, 28C FEBRUARY, 1980. 29C 30C DESCRIPTION OF PARAMETERS. 31C 32C NG = NUMBER OF FUNCTIONS GI, OR THE NUMBER OF COMPONENTS OF 33C THE VECTOR VALUED FUNCTION G(X). INPUT ONLY. 34C 35C HMIN = RESOLUTION PARAMETER IN X. INPUT ONLY. WHEN A ROOT IS 36C FOUND, IT IS LOCATED ONLY TO WITHIN AN ERROR OF HMIN IN X. 37C TYPICALLY, HMIN SHOULD BE SET TO SOMETHING ON THE ORDER OF 38C 100 * UROUND * MAX(ABS(X0),ABS(X1)), 39C WHERE UROUND IS THE UNIT ROUNDOFF OF THE MACHINE. 40C 41C JFLAG = INTEGER FLAG FOR INPUT AND OUTPUT COMMUNICATION. 42C 43C ON INPUT, SET JFLAG = 0 ON THE FIRST CALL FOR THE PROBLEM, 44C AND LEAVE IT UNCHANGED UNTIL THE PROBLEM IS COMPLETED. 45C (THE PROBLEM IS COMPLETED WHEN JFLAG .GE. 2 ON RETURN.) 46C 47C ON OUTPUT, JFLAG HAS THE FOLLOWING VALUES AND MEANINGS.. 48C JFLAG = 1 MEANS DROOTS NEEDS A VALUE OF G(X). SET GX = G(X) 49C AND CALL DROOTS AGAIN. 50C JFLAG = 2 MEANS A ROOT HAS BEEN FOUND. THE ROOT IS 51C AT X, AND GX CONTAINS G(X). (ACTUALLY, X IS THE 52C RIGHTMOST APPROXIMATION TO THE ROOT ON AN INTERVAL 53C (X0,X1) OF SIZE HMIN OR LESS.) 54C JFLAG = 3 MEANS X = X1 IS A ROOT, WITH ONE OR MORE OF THE GI 55C BEING ZERO AT X1 AND NO SIGN CHANGES IN (X0,X1). 56C GX CONTAINS G(X) ON OUTPUT. 57C JFLAG = 4 MEANS NO ROOTS (OF ODD MULTIPLICITY) WERE 58C FOUND IN (X0,X1) (NO SIGN CHANGES). 59C 60C X0,X1 = ENDPOINTS OF THE INTERVAL WHERE ROOTS ARE SOUGHT. 61C X1 AND X0 ARE INPUT WHEN JFLAG = 0 (FIRST CALL), AND 62C MUST BE LEFT UNCHANGED BETWEEN CALLS UNTIL THE PROBLEM IS 63C COMPLETED. X0 AND X1 MUST BE DISTINCT, BUT X1 - X0 MAY BE 64C OF EITHER SIGN. HOWEVER, THE NOTION OF -LEFT- AND -RIGHT- 65C WILL BE USED TO MEAN NEARER TO X0 OR X1, RESPECTIVELY. 66C WHEN JFLAG .GE. 2 ON RETURN, X0 AND X1 ARE OUTPUT, AND 67C ARE THE ENDPOINTS OF THE RELEVANT INTERVAL. 68C 69C G0,G1 = ARRAYS OF LENGTH NG CONTAINING THE VECTORS G(X0) AND G(X1), 70C RESPECTIVELY. WHEN JFLAG = 0, G0 AND G1 ARE INPUT AND 71C NONE OF THE G0(I) SHOULD BE BE ZERO. 72C WHEN JFLAG .GE. 2 ON RETURN, G0 AND G1 ARE OUTPUT. 73C 74C GX = ARRAY OF LENGTH NG CONTAINING G(X). GX IS INPUT 75C WHEN JFLAG = 1, AND OUTPUT WHEN JFLAG .GE. 2. 76C 77C X = INDEPENDENT VARIABLE VALUE. OUTPUT ONLY. 78C WHEN JFLAG = 1 ON OUTPUT, X IS THE POINT AT WHICH G(X) 79C IS TO BE EVALUATED AND LOADED INTO GX. 80C WHEN JFLAG = 2 OR 3, X IS THE ROOT. 81C WHEN JFLAG = 4, X IS THE RIGHT ENDPOINT OF THE INTERVAL, X1. 82C 83C JROOT = INTEGER ARRAY OF LENGTH NG. OUTPUT ONLY. 84C WHEN JFLAG = 2 OR 3, JROOT INDICATES WHICH COMPONENTS 85C OF G(X) HAVE A ROOT AT X. JROOT(I) IS 1 IF THE I-TH 86C COMPONENT HAS A ROOT, AND JROOT(I) = 0 OTHERWISE. 87C 88C IMAX, LAST, ALPHA, X2 = 89C BOOKKEEPING VARIABLES WHICH MUST BE SAVED FROM CALL 90C TO CALL. THEY ARE SAVED INSIDE THE CALLING ROUTINE, 91C BUT THEY ARE USED ONLY WITHIN THIS ROUTINE. 92C----------------------------------------------------------------------- 93 INTEGER I, IMXOLD, NXLAST 94 DOUBLE PRECISION T2, TMAX, ZERO 95 LOGICAL ZROOT, SGNCHG, XROOT 96 DATA ZERO/0.0D0/ 97C 98 IF (JFLAG .EQ. 1) GO TO 200 99C JFLAG .NE. 1. CHECK FOR CHANGE IN SIGN OF G OR ZERO AT X1. ---------- 100 IMAX = 0 101 TMAX = ZERO 102 ZROOT = .FALSE. 103 DO 120 I = 1,NG 104 IF (DABS(G1(I)) .GT. ZERO) GO TO 110 105 ZROOT = .TRUE. 106 GO TO 120 107C AT THIS POINT, G0(I) HAS BEEN CHECKED AND CANNOT BE ZERO. ------------ 108 110 IF (DSIGN(1.0D0,G0(I)) .EQ. DSIGN(1.0D0,G1(I))) GO TO 120 109 T2 = DABS(G1(I)/(G1(I)-G0(I))) 110 IF (T2 .LE. TMAX) GO TO 120 111 TMAX = T2 112 IMAX = I 113 120 CONTINUE 114 IF (IMAX .GT. 0) GO TO 130 115 SGNCHG = .FALSE. 116 GO TO 140 117 130 SGNCHG = .TRUE. 118 140 IF (.NOT. SGNCHG) GO TO 400 119C THERE IS A SIGN CHANGE. FIND THE FIRST ROOT IN THE INTERVAL. -------- 120 XROOT = .FALSE. 121 NXLAST = 0 122 LAST = 1 123C 124C REPEAT UNTIL THE FIRST ROOT IN THE INTERVAL IS FOUND. LOOP POINT. --- 125 150 CONTINUE 126 IF (XROOT) GO TO 300 127 IF (NXLAST .EQ. LAST) GO TO 160 128 ALPHA = 1.0D0 129 GO TO 180 130 160 IF (LAST .EQ. 0) GO TO 170 131 ALPHA = 0.5D0*ALPHA 132 GO TO 180 133 170 ALPHA = 2.0D0*ALPHA 134 180 X2 = X1 - (X1-X0)*G1(IMAX)/(G1(IMAX) - ALPHA*G0(IMAX)) 135 IF ((DABS(X2-X0) .LT. HMIN) .AND. 136 1 (DABS(X1-X0) .GT. 10.0D0*HMIN)) X2 = X0 + 0.1D0*(X1-X0) 137 JFLAG = 1 138 X = X2 139C RETURN TO THE CALLING ROUTINE TO GET A VALUE OF GX = G(X). ----------- 140 RETURN 141C CHECK TO SEE IN WHICH INTERVAL G CHANGES SIGN. ----------------------- 142 200 IMXOLD = IMAX 143 IMAX = 0 144 TMAX = ZERO 145 ZROOT = .FALSE. 146 DO 220 I = 1,NG 147 IF (DABS(GX(I)) .GT. ZERO) GO TO 210 148 ZROOT = .TRUE. 149 GO TO 220 150C NEITHER G0(I) NOR GX(I) CAN BE ZERO AT THIS POINT. ------------------- 151 210 IF (DSIGN(1.0D0,G0(I)) .EQ. DSIGN(1.0D0,GX(I))) GO TO 220 152 T2 = DABS(GX(I)/(GX(I) - G0(I))) 153 IF (T2 .LE. TMAX) GO TO 220 154 TMAX = T2 155 IMAX = I 156 220 CONTINUE 157 IF (IMAX .GT. 0) GO TO 230 158 SGNCHG = .FALSE. 159 IMAX = IMXOLD 160 GO TO 240 161 230 SGNCHG = .TRUE. 162 240 NXLAST = LAST 163 IF (.NOT. SGNCHG) GO TO 250 164C SIGN CHANGE BETWEEN X0 AND X2, SO REPLACE X1 WITH X2. ---------------- 165 X1 = X2 166 CALL DCOPY (NG, GX, 1, G1, 1) 167 LAST = 1 168 XROOT = .FALSE. 169 GO TO 270 170 250 IF (.NOT. ZROOT) GO TO 260 171C ZERO VALUE AT X2 AND NO SIGN CHANGE IN (X0,X2), SO X2 IS A ROOT. ----- 172 X1 = X2 173 CALL DCOPY (NG, GX, 1, G1, 1) 174 XROOT = .TRUE. 175 GO TO 270 176C NO SIGN CHANGE BETWEEN X0 AND X2. REPLACE X0 WITH X2. --------------- 177 260 CONTINUE 178 CALL DCOPY (NG, GX, 1, G0, 1) 179 X0 = X2 180 LAST = 0 181 XROOT = .FALSE. 182 270 IF (DABS(X1-X0) .LE. HMIN) XROOT = .TRUE. 183 GO TO 150 184C 185C RETURN WITH X1 AS THE ROOT. SET JROOT. SET X = X1 AND GX = G1. ----- 186 300 JFLAG = 2 187 X = X1 188 CALL DCOPY (NG, G1, 1, GX, 1) 189 DO 320 I = 1,NG 190 JROOT(I) = 0 191 IF (DABS(G1(I)) .GT. ZERO) GO TO 310 192 JROOT(I) = 1 193 GO TO 320 194 310 IF (DSIGN(1.0D0,G0(I)) .NE. DSIGN(1.0D0,G1(I))) JROOT(I) = 1 195 320 CONTINUE 196 RETURN 197C 198C NO SIGN CHANGE IN THE INTERVAL. CHECK FOR ZERO AT RIGHT ENDPOINT. --- 199 400 IF (.NOT. ZROOT) GO TO 420 200C 201C ZERO VALUE AT X1 AND NO SIGN CHANGE IN (X0,X1). RETURN JFLAG = 3. --- 202 X = X1 203 CALL DCOPY (NG, G1, 1, GX, 1) 204 DO 410 I = 1,NG 205 JROOT(I) = 0 206 IF (DABS(G1(I)) .LE. ZERO) JROOT (I) = 1 207 410 CONTINUE 208 JFLAG = 3 209 RETURN 210C 211C NO SIGN CHANGES IN THIS INTERVAL. SET X = X1, RETURN JFLAG = 4. ----- 212 420 CALL DCOPY (NG, G1, 1, GX, 1) 213 X = X1 214 JFLAG = 4 215 RETURN 216C---------------------- END OF SUBROUTINE DROOTS ----------------------- 217 END 218