This subroutine finds the leftmost root of a set of arbitrary functions gi(x) (i = 1,…,NG) in an interval (X0,X1). Only roots of odd multiplicity (i.e. changes of sign of the gi) are found. Here the sign of X1 - X0 is arbitrary, but is constant for a given problem, and -leftmost- means nearest to X0. The values of the vector-valued function g(x) = (gi, i=1…NG) are communicated through the call sequence of DROOTS. The method used is the Illinois algorithm.
#### Reference: Kathie L. Hiebert and Lawrence F. Shampine, Implicitly Defined Output Points for Solutions of ODEs, Sandia Report SAND80-0180, February 1980.
number of functions gi, or the number of components of the vector valued function g(x). Input only.
resolution parameter in X. Input only. When a root is found, it is located only to within an error of HMIN in X. Typically, HMIN should be set to something on the order of 100 * UROUND * MAX(ABS(X0),ABS(X1)), where UROUND is the unit roundoff of the machine.
integer flag for input and output communication.
On input, set JFLAG = 0 on the first call for the problem, and leave it unchanged until the problem is completed. (The problem is completed when JFLAG .ge. 2 on return.)
JFLAG
:JFLAG has the following values and meanings:
JFLAG = 1 means DROOTS needs a value of g(x). Set GX = g(X)
and call DROOTS again.
JFLAG = 2 means a root has been found. The root is
at X, and GX contains g(X). (Actually, X is the
rightmost approximation to the root on an interval
(X0,X1) of size HMIN or less.)
JFLAG = 3 means X = X1 is a root, with one or more of the gi
being zero at X1 and no sign changes in (X0,X1).
GX contains g(X) on output.
JFLAG = 4 means no roots (of odd multiplicity) were
found in (X0,X1) (no sign changes).
endpoints of the interval where roots are sought. X1 and X0 are input when JFLAG = 0 (first call), and must be left unchanged between calls until the problem is completed. X0 and X1 must be distinct, but X1 - X0 may be of either sign. However, the notion of -left- and -right- will be used to mean nearer to X0 or X1, respectively. When JFLAG .ge. 2 on return, X0 and X1 are output, and are the endpoints of the relevant interval.
arrays of length NG containing the vectors g(X0) and g(X1), respectively. When JFLAG = 0, G0 and G1 are input and none of the G0(i) should be zero. When JFLAG .ge. 2 on return, G0 and G1 are output.
array of length NG containing g(X). GX is input when JFLAG = 1, and output when JFLAG .ge. 2.
independent variable value. Output only. When JFLAG = 1 on output, X is the point at which g(x) is to be evaluated and loaded into GX.
When JFLAG = 2 or 3, X is the root.
When JFLAG = 4, X is the right endpoint of the interval, X1.
integer array of length NG. Output only. When JFLAG = 2 or 3, JROOT indicates which components of g(x) have a root at X. JROOT(i) is 1 if the i-th component has a root, and JROOT(i) = 0 otherwise.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer | :: | Ng | ||||
| real(kind=dp), | intent(in) | :: | Hmin | |||
| integer, | intent(inout) | :: | Jflag | |||
| real(kind=dp), | intent(inout) | :: | X0 | |||
| real(kind=dp), | intent(inout) | :: | X1 | |||
| real(kind=dp) | :: | G0(Ng) | ||||
| real(kind=dp) | :: | G1(Ng) | ||||
| real(kind=dp) | :: | Gx(Ng) | ||||
| real(kind=dp), | intent(out) | :: | X | |||
| integer, | intent(out) | :: | Jroot(Ng) |
| Type | Visibility | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|---|
| real(kind=dp), | public, | parameter | :: | five | = | 5.0d0 | |
| real(kind=dp), | public | :: | fracint | ||||
| real(kind=dp), | public | :: | fracsub | ||||
| real(kind=dp), | public, | parameter | :: | half | = | 0.5d0 | |
| integer, | public | :: | i | ||||
| integer, | public | :: | imxold | ||||
| integer, | public | :: | nxlast | ||||
| logical, | public | :: | sgnchg | ||||
| real(kind=dp), | public | :: | t2 | ||||
| real(kind=dp), | public, | parameter | :: | tenth | = | 0.1d0 | |
| real(kind=dp), | public | :: | tmax | ||||
| logical, | public | :: | xroot | ||||
| real(kind=dp), | public, | parameter | :: | zero | = | 0.0d0 | |
| logical, | public | :: | zroot |
subroutine droots(Ng,Hmin,Jflag,X0,X1,G0,G1,Gx,X,Jroot) ! integer :: Ng real(kind=dp),intent(in) :: Hmin integer,intent(inout) :: Jflag real(kind=dp),intent(inout) :: X0 real(kind=dp),intent(inout) :: X1 real(kind=dp) :: G0(Ng) real(kind=dp) :: G1(Ng) real(kind=dp) :: Gx(Ng) real(kind=dp),intent(out) :: X integer,intent(out) :: Jroot(Ng) ! real(kind=dp) :: fracint , fracsub , t2 , tmax integer :: i , imxold , nxlast logical :: sgnchg , xroot , zroot ! real(kind=dp) , parameter :: five=5.0d0 , half=0.5d0 , tenth=0.1d0 , zero=0.0d0 ! if ( Jflag==1 ) then ! Check to see in which interval g changes sign. ----------------------- imxold = dlsr%imax dlsr%imax = 0 tmax = zero zroot = .false. do i = 1 , Ng if ( abs(Gx(i))<=zero ) then zroot = .true. ! Neither G0(i) nor GX(i) can be zero at this point. ------------------- elseif ( sign(1.0D0,G0(i))/=sign(1.0D0,Gx(i)) ) then t2 = abs(Gx(i)/(Gx(i)-G0(i))) if ( t2>tmax ) then tmax = t2 dlsr%imax = i endif endif enddo if ( dlsr%imax>0 ) then sgnchg = .true. else sgnchg = .false. dlsr%imax = imxold endif nxlast = dlsr%last if ( sgnchg ) then ! Sign change between X0 and X2, so replace X1 with X2. ---------------- X1 = dlsr%x2 !X!call dcopy(Ng,Gx,1,G1,1) G1(1:Ng) = Gx(1:Ng)!X! dlsr%last = 1 xroot = .false. elseif ( .not.zroot ) then ! No sign change between X0 and X2. Replace X0 with X2. --------------- !X!call dcopy(Ng,Gx,1,G0,1) G0(1:Ng) = Gx(1:Ng)!X! X0 = dlsr%x2 dlsr%last = 0 xroot = .false. else ! Zero value at X2 and no sign change in (X0,X2), so X2 is a root. ----- X1 = dlsr%x2 !X!call dcopy(Ng,Gx,1,G1,1) G1(1:Ng) = Gx(1:Ng)!X! xroot = .true. endif if ( abs(X1-X0)<=Hmin ) xroot = .true. else ! JFLAG .ne. 1. Check for change in sign of g or zero at X1. ---------- dlsr%imax = 0 tmax = zero zroot = .false. do i = 1 , Ng if ( abs(G1(i))<=zero ) then zroot = .true. ! At this point, G0(i) has been checked and cannot be zero. ------------ elseif ( sign(1.0D0,G0(i))/=sign(1.0D0,G1(i)) ) then t2 = abs(G1(i)/(G1(i)-G0(i))) if ( t2>tmax ) then tmax = t2 dlsr%imax = i endif endif enddo if ( dlsr%imax>0 ) then sgnchg = .true. else sgnchg = .false. endif if ( .not.sgnchg ) then ! ! No sign change in the interval. Check for zero at right endpoint. --- if ( zroot ) then ! ! Zero value at X1 and no sign change in (X0,X1). Return JFLAG = 3. --- X = X1 !X!call dcopy(Ng,G1,1,Gx,1) Gx(1:Ng) = G1(1:Ng)!X! do i = 1 , Ng Jroot(i) = 0 if ( abs(G1(i))<=zero ) Jroot(i) = 1 enddo Jflag = 3 return endif ! ! No sign changes in this interval. Set X = X1, return JFLAG = 4. ----- !X!call dcopy(Ng,G1,1,Gx,1) Gx(1:Ng) = G1(1:Ng)!X! X = X1 Jflag = 4 return else ! There is a sign change. Find the first root in the interval. -------- xroot = .false. nxlast = 0 dlsr%last = 1 endif endif ! ! Repeat until the first root in the interval is found. Loop point. --- if ( xroot ) then ! ! Return with X1 as the root. Set JROOT. Set X = X1 and GX = G1. ----- Jflag = 2 X = X1 !X!call dcopy(Ng,G1,1,Gx,1) Gx(1:Ng) = G1(1:Ng)!X! do i = 1 , Ng Jroot(i) = 0 if ( abs(G1(i))>zero ) then if ( sign(1.0D0,G0(i))/=sign(1.0D0,G1(i)) ) Jroot(i) = 1 else Jroot(i) = 1 endif enddo else if ( nxlast/=dlsr%last ) then dlsr%alpha = 1.0D0 elseif ( dlsr%last==0 ) then dlsr%alpha = 2.0D0*dlsr%alpha else dlsr%alpha = 0.5D0*dlsr%alpha endif dlsr%x2 = X1 - (X1-X0)*G1(dlsr%imax)/(G1(dlsr%imax)-dlsr%alpha*G0(dlsr%imax)) ! If X2 is too close to X0 or X1, adjust it inward, by a fractional ---- ! distance that is between 0.1 and 0.5. -------------------------------- if ( abs(dlsr%x2-X0)<half*Hmin ) then fracint = abs(X1-X0)/Hmin fracsub = tenth if ( fracint<=five ) fracsub = half/fracint dlsr%x2 = X0 + fracsub*(X1-X0) endif if ( abs(X1-dlsr%x2)<half*Hmin ) then fracint = abs(X1-X0)/Hmin fracsub = tenth if ( fracint<=five ) fracsub = half/fracint dlsr%x2 = X1 - fracsub*(X1-X0) endif Jflag = 1 X = dlsr%x2 ! Return to the calling routine to get a value of GX = g(X). ----------- endif end subroutine droots