subroutine aptqrts (noptq, aa, bb, cc, qq, tol,
& nroots, root1, root2, itrun)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQRTS
c
c call aptqrts (noptq, aa, bb, cc, qq, tol,
c & nroots, root1, root2, itrun)
c
c Version: aptqrts Updated 1994 December 15 12:40.
c aptqrts Originated 1989 November 2 14:10.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the real and complex roots of the quadratic equation:
c
c aa * x**2 + bb * x + cc = 0.
c
c The method minimizes truncation error, and indicates when
c truncation error may still be significant, based on tol.
c Option noptq = 1 allows the user to specify the value of
c qq = bb**2 - 4.0 * aa *cc, instead of using the value
c calculated here. This allows any cancellation of terms.
c
c For q = sqrt (qq), the roots are:
c (-bb + q) / (2 * aa) = 2 * cc / (-bb - q)
c (-bb - q) / (2 * aa) = 2 * cc / (-bb + q)
c The most accurate form must be used, to minimize error.
c
c Note: aptqrtv is a vectorized version of aptqrts.
c
c History: 1990 March 21 14:00. Changed to truncate small positive values
c of qq to zero.
c 1994 April 15 16:20. Changed to put two roots in order of
c increasing value.
c 1994 December 15 12:40. Changed to return complex roots.
c
c Input: noptq, aa, bb, cc, qq, tol.
c
c Output: qq, nroots, root1, root2, itrun.
c
c Glossary:
c
c aa Input Coefficient of x**2 in a quadratic equation.
c
c bb Input Coefficient of x in a quadratic equation.
c
c cc Input Coefficient of 1 in a quadratic equation.
c
c itrun Output Truncation error indicator. 0 if insignificant.
c 1 if the magnitude of qq is less than the estimated
c truncation error, based on tol. This indicates
c that the roots are near a minimum or maximum, or
c that the quadratic is almost the square of a linear
c function. The value of qq is truncated to zero.
c 2 if the input value of qq (noptq = 1) differs from
c the calculated value by more than the estimated
c truncation error, based on tol.
c
c noptq Input Option for getting value of qq = bb**2 - 4.0 * aa * cc:
c 0 to calculate from the input values of aa, bb, cc.
c Any value other than 1 has this result.
c 1 to use the (more accurate) input value of qq.
c A nonzero value of tol must be used with this
c option, for comparision of the input with the
c calculated value.
c
c nroots Output Number of real roots:
c -1 if aa = 0 and bb = 0.
c Equation is null (cc = 0) or bad (cc is not 0).
c 0 if qq < 0. The two complex roots are
c root1 + i * root2 and root1 - i * root2.
c 1 if aa = 0 and bb is not 0, or if qq = 0.
c 2 if qq > 0.
c
c qq Input If noptq = 1, qq = bb**2 - 4.0 * aa * cc, but
c calculated more accurately, due to cancellation of
c terms resulting from the composite nature of aa, bb,
c and/or cc.
c
c qq Output If noptq is not 1, qq = bb**2 - 4.0 * aa * cc,
c calculated locally.
c
c root1 Output If nroots = 1, the only real root.
c If nroots = 2, the minimum real root.
c If nroots = 0, the two complex roots are
c root1 + i * root2 and root1 - i * root2.
c
c root2 Output If nroots is 2, the maximum real root,
c If nroots = 0, the two complex roots are
c root1 + i * root2 and root1 - i * root2.
c
c tol Input Numerical tolerance limit.
c Must not be zero if noptq = 1.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c.... (None.)
c.... Local variables.
c---- A very big number.
common /laptqrts/ big
c---- Sqrt (qq), with the sign of bb.
common /laptqrts/ q
c---- Truncation error in qqx.
common /laptqrts/ qqerr
c---- Locally calculated value of qq.
common /laptqrts/ qqx
c---- One of two roots.
common /laptqrts/ rootx
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptqrts finding roots. noptq=',i2,' qq=',1pe22.14 /
cbug & ' aa,bb,cc=',1p3e22.14)
cbug write ( 3, 9901) noptq, qq, aa, bb, cc
cbugc***DEBUG ends.
c.... Initialize.
c---- A very big number.
big = 1.e+99
itrun = 0
root1 = -big
root2 = -big
c.... Find any real roots.
c---- Equation is not a quadratic.
if (aa .eq. 0.0) then
c---- Equation is null or wrong. No roots.
if (bb .eq. 0.0) then
nroots = -1
c---- Equation is linear. 1 root exists.
else
nroots = 1
root1 = -cc / bb
c---- Branch on value of bb.
endif
c---- Equation is quadratic.
else
qqerr = 2.0 * tol * (bb**2 + 4.0 * abs (aa * cc))
qqx = bb**2 - 4.0 * aa * cc
c---- Calculate the value of qq.
if (noptq .ne. 1) then
qq = qqx
c++++ Truncation error significant.
if (abs (qq) .lt. qqerr) then
itrun = 1
qq = 0.0
endif
c---- Use the input value of qq.
else
c++++ Bad input value of qq.
if (abs (qq - qqx) .gt. qqerr) then
itrun = 2
endif
endif
c---- Equation has no real roots. Return real and imaginary parts.
if (qq .lt. 0.0) then
nroots = 0
root1 = -0.5 * bb / aa
root2 = 0.5 * sqrt (-qq) / abs (aa)
c---- Equation has one degenerate root.
elseif (qq .eq. 0.0) then
nroots = 1
root1 = -0.5 * bb / aa
c---- Equation has two real roots.
else
nroots = 2
q = sign (sqrt (qq), bb)
root1 = -2.0 * cc / (q + bb)
root2 = -0.5 * (q + bb) / aa
if (root2 .lt. root1) then
rootx = root1
root1 = root2
root2 = rootx
endif
c---- Branch on value of qq.
endif
c---- Branch on value of aa.
endif
cbugc***DEBUG begins.
cbug 9902 format ('aptqrts results. nroots=',i2,' itrun=',i2 /
cbug & ' root1,root2=',1p2e22.14)
cbug write ( 3, 9902) nroots, itrun, root1, root2
cbugc***DEBUG ends.
return
c.... End of subroutine aptqrts. (+1 line.)
end
UCRL-WEB-209832