A global bracket finder. For a given function it finds values on each side of each of the function's roots within a given range.
| Type | Intent | Optional | Attributes | Name | ||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| real(kind=8), | intent(out), | dimension(:,:) | :: | brackets | A 2 by n real array in which the left and right brackets will be stored and returned. Will find up to n sets of brackets. |
|||||||||||||||
| real(kind=8), | intent(inout) | :: | dx | The initial size of increment to use when examining function. Minimum interval will be 0.01 of this. |
||||||||||||||||
public pure function fnctn(x)The Fortran function for which the brackets will be found. Must take only one real argument, return a real argument. Arguments
Return Value real(kind=8) |
||||||||||||||||||||
| integer, | intent(out) | :: | numroots | An integer value which will return the number of roots for which brackets were found. A negative number indicates that an error occurred. |
||||||||||||||||
| logical, | intent(in) | :: | verbose | A logical variable which specifies whether to print to stdout any bracket values which are found and warning messages when 'dx' set to 'dxmin'. |
||||||||||||||||
| real(kind=8), | intent(in) | :: | xmax | The upper limit of the range on which the subroutine will search for roots and brackets. |
||||||||||||||||
| real(kind=8), | intent(in) | :: | xmin | The lower limit of the range on which the subroutine will search for roots and brackets. |
||||||||||||||||
SUBROUTINE global_brackets ( brackets, dx , fnctn , numroots, &
verbose , xmax , xmin )
!! A global bracket finder. For a given function it finds values
!! on each side of each of the function's roots within a given
!! range.
IMPLICIT NONE
interface
pure real(8) function fnctn ( x )
!!The Fortran function for which the brackets will be
!! found. Must take only one real argument, return a
!! real argument.
real(8), intent(in) :: x
end function fnctn
end interface
! Input and output variables:
logical, intent(in) :: verbose
!! A logical variable which specifies whether to print to
!! stdout any bracket values which are found and warning
!! messages when 'dx' set to 'dxmin'.
real(8), intent(in) :: xmax
!! The upper limit of the range on which the subroutine will
!! search for roots and brackets.
real(8), intent(in) :: xmin
!! The lower limit of the range on which the subroutine will
!! search for roots and brackets.
real(8), intent(inout) :: dx
!! The initial size of increment to use when examining
!! function. Minimum interval will be 0.01 of this.
integer, intent(out) :: numroots
!! An integer value which will return the number of roots
!! for which brackets were found. A negative number
!! indicates that an error occurred.
real(8), intent(out), dimension(:,:) :: brackets
!! A 2 by n real array in which the left and right brackets
!! will be stored and returned. Will find up to n sets of
!! brackets.
! Other variables:
REAL(8) :: dfdx = 0.d0, &
dfdxh = 0.d0, &
dldx = 0.d0, &
dldxh = 0.d0, &
dxmin = 0.d0, &
dxh = 0.d0, &
fleft = 0.d0, &
fmid = 0.d0, &
fright = 0.d0, &
scaleval = 0.d0, &
xleft = 0.d0, &
xright = 0.d0, &
xmid = 0.d0
!------------------------------------------------------------------!
! Check if passed array is big enough to hold returned data
IF ( SIZE(brackets, 1) < 2 ) THEN
numroots = -1
WRITE( 0,2000)
RETURN
END IF
! Initialize variables
dxmin = 1.d-2 * dx
numroots = 0
xleft = xmin
fleft = fnctn(xmin)
! Repeat this loop until have traversed range [xmin:xmax], or
! until brackets array is filled
DO WHILE ( xleft <= xmax )
! Make sure that xright <= xmax
dx = MIN(dx, xmax - xmin)
! Get values and derivatives for next step along the
! function
xright = xleft + dx
fright = fnctn(xright)
scaleval = MAX(1.0d0, ABS(xleft)) / MAX(ABS(fleft), &
ABS(fright))
dfdx = scaleval * (fright - fleft) / dx
dldx = SIGN(1.0d0, dfdx) * SQRT(1.0d0 + dfdx**2.d0)
! Find derivative of secant length (called 'dldx' when found
! above) at midpoint between 'xleft' and 'xright' ('dldxh').
! If difference between 'dldx' and 'dldxh' is too great,
! reduce 'dx' and repeat.
DO
dxh = 5d-1 * dx
xmid = xleft + dxh
fmid = fnctn(xmid)
dfdxh = scaleval * (fmid - fleft) / dxh
dldxh = SIGN(1.0d0, dfdxh) * SQRT(1.0d0 + dfdxh**2)
IF ( ( ABS(dldx - dldxh) > 5.d-1*ABS(dldx) ) &
.AND. ( dxh > dxmin ) ) THEN
dx = dxh
xright = xmid
fright = fmid
dldx = dldxh
ELSE
EXIT
END IF
END DO
! If difference between 'dldx' and 'dldxh' is too small,
! make step-size larger for next iteration
IF ( ABS(dldx - dldxh) < 1.d-1 * ABS(dldx) ) THEN
dx = 1.5d0 * dx
! If diference between 'dldx' and 'dldxh' is too small, but
! can't go smaller without falling below 'dxmin', set 'dx =
! dxmin' and print warning message
ELSE IF ( ( ABS(dldx - dldxh) > 5.d-1 * ABS(dldx) ) &
.AND. ( dxh < dxmin ) ) THEN
dx = dxmin
xright = xleft + dx
fright = fnctn(xright)
IF ( verbose .EQV. .TRUE. ) WRITE( 6,2010) dxmin
END IF
! Otherwise keep 'dx' the same (no action needs to be taken)
! If the signs of 'fleft' and 'fright' are opposite, and if
! signs of 'fright' and the derivative of 'fnctn(xleft)' are
! the same, then 'xleft' and 'xright' bracket a root
IF ( ( fleft * fright < 0 ) .AND. ( fright * (fleft - &
fnctn(xleft - dxmin) ) > 0 ) ) THEN
numroots = numroots + 1
brackets(1:2,numroots) = (/ xleft, xright /)
IF ( verbose .EQV. .TRUE. ) WRITE( 6,2020) xleft, xright
! If have filled 'brackets', then leave loop
IF ( numroots >= SIZE(brackets,2) ) EXIT
END IF
! Update variables for next iteration
xleft = xright
fleft = fright
END DO
!--------------------------------------------------------------!
! Write format statements !
!--------------------------------------------------------------!
2000 FORMAT('GLOBRACK: ERROR: Passed array too small to hold ',&
'returned values. Must have at least 2 columns.')
2010 FORMAT('GLOBRACK: WARNING: Value of dx set to minimum ' &
,'value: ',1PG22.15)
2020 FORMAT('GLOBRACK: Found root bracketed by ',1PG22.15, &
' and',/, &
'GLOBRACK: ',1PG22.15)
!--------------------------------------------------------------!
RETURN
END SUBROUTINE global_brackets