subroutine aptqrtv (noptq, aa, bb, cc, qq, np, tol,
& nroots, root1, root2, itrun, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQRTV
c
c call aptqrtv (noptq, aa, bb, cc, qq, np, tol,
c & nroots, root1, root2, itrun, nerr)
c
c Version: aptqrtv Updated 1994 December 15 12:40.
c aptqrtv 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 for np sets of coefficients aa, bb, cc.
c The solution is vectorized over the sets of coefficients,
c which are calculated in sections of 64 (or less, for the
c final section).
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: aptqrts is a scalar version of aptqrtv.
c
c History: 1990 January 19 15:00. Fixed bug for a = 0.0 or b = 0.0.
c 1990 March 21 14:00. Changed to allow truncation of small
c positive qq to zero.
c 1994 April 15 16:20. Changed to put two roots in order of
c increasing value.
c
c Input: noptq, aa, bb, cc, qq, np, tol.
c
c Output: qq, nroots, root1, root2, itrun.
c
c Glossary:
c
c aa Input Coefficients of x**2 in a quadratic equation. Size np.
c
c bb Input Coefficients of x in a quadratic equation. Size np.
c
c cc Input Coefficients of 1 in a quadratic equation. Size np.
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 Size np.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c noptq Input Option for getting value of qq = bb**2 - 4.0 * aa * cc:
c 0 to not use input qq, but calculate from 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. Size np.
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 equation is linear:
c aa is 0 and bb is not 0, or if qq = 0.
c 2 if qq > 0.
c
c np Input Number of sets of input data aa, bb, cc, noptq, qq.
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. Size np.
c
c qq Output If noptq is not 1, qq = bb**2 - 4.0 * aa * cc,
c calculated locally. Size np.
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 = 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---- Coefficient of term x**2.
dimension aa (1)
c---- Coefficient of term x.
dimension bb (1)
c---- Coefficient of term 1.
dimension cc (1)
c---- Truncation error indicator.
dimension itrun (1)
c---- Number of real roots.
dimension nroots (1)
c---- Factor bb**2 - 4.0 * aa * cc.
dimension qq (1)
c---- First real root (-1.e+99 if none).
dimension root1 (1)
c---- Second real root (-1.e+99 if none).
dimension root2 (1)
c.... Local variables.
c---- A very big number.
common /laptqrtv/ big
c---- A very small number.
common /laptqrtv/ fuz
c---- Index of coefficient set.
common /laptqrtv/ n
c---- First index of subset of data.
common /laptqrtv/ n1
c---- Last index of subset of data.
common /laptqrtv/ n2
c---- Index in external array.
common /laptqrtv/ nn
c---- Size of current subset of data.
common /laptqrtv/ ns
c---- Sqrt (qq), with the sign of bb.
common /laptqrtv/ q
c---- Truncation error in qqx.
common /laptqrtv/ qqerr
c---- Locally calculated value of qq.
common /laptqrtv/ qqx
c---- One of two roots.
common /laptqrtv/ rootx
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptqrtv finding roots. np=',i3,' noptq=',i2,
cbug & ' tol=',1pe13.5)
cbug 9902 format (
cbug & i3,' aa,bb,cc= ',1p3e22.14 /
cbug & ' qq(input)=',1pe22.14 /
cbug & ' qqx(calc)=',1pe22.14 /
cbug & ' qqerr= ',1pe22.14)
cbug write ( 3, 9901) np, noptq, tol
cbug do 101 n = 1, np
cbug if (noptq .ne. 1) qq(n) = 0.0
cbug qqx = bb(n)**2 - 4.0 * aa(n) * cc(n)
cbug qqerr = 2.0 * tol *
cbug & (bb(n)**2 + 4.0 * abs (aa(n) * cc(n)))
cbug write ( 3, 9902) n, aa(n), bb(n), cc(n), qq(n), qqx, qqerr
cbug 101 continue
cbugc***DEBUG ends.
c.... Initialize.
c---- A very big number.
big = 1.e+99
c---- A very small number.
fuz = 1.e-99
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
c.... Initialize.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Find qq = bb**2 - 4.0 * aa * cc, and q = sqrt (qq).
c---- Loop over coefficient sets.
do 120 nn = n1, n2
qqx = bb(nn)**2 - 4.0 * aa(nn) * cc(nn)
qqerr = 2.0 * tol *
& (bb(nn)**2 + 4.0 * abs (aa(nn) * cc(nn)))
if (abs (qqx) .lt. qqerr) then
qqx = 0.0
endif
if (abs (qqx) .lt. qqerr) then
itrun(nn) = 1
else
itrun(nn) = 0
endif
if (noptq .ne. 1) then
qq(nn) = qqx
endif
if (abs (qq(nn) - qqx) .gt. qqerr) then
itrun(nn) = 2
endif
c---- End of loop over coefficient sets.
120 continue
c.... Find any roots.
c---- Loop over coefficient sets.
do 130 nn = n1, n2
q = sign (sqrt (amax1 (0.0, qq(nn))), bb(nn))
c.... No first root if qq .lt. 0, or if aa = bb = 0.
root1(nn) = -2.0 * cc(nn) / (q + bb(nn) + fuz)
if ((qq(nn) .lt. 0.0) .or.
& ((aa(nn) .eq. 0.0) .and.
& (bb(nn) .eq. 0.0)) ) then
root1(nn) = -big
endif
if (root1(nn) .ne. -big) then
nroots(nn) = 1
else
nroots(nn) = 0
endif
c.... No second real root if qq = 0, qq < 0, or if aa = 0.
root2(nn) = -0.5 * (q + bb(nn)) / (aa(nn) + fuz)
if ((qq(nn) .le. 0.0) .or.
& (aa(nn) .eq. 0.0) ) then
root2(nn) = -big
endif
if (root2(nn) .ne. -big) then
nroots(nn) = nroots(nn) + 1
endif
c.... Flag bad coefficients when aa = bb = 0.
if ((aa(nn) .eq. 0.0) .and.
& (bb(nn) .eq. 0.0) ) then
nroots(nn) = -1
endif
if ((nroots(nn) .eq. 2) .and.
& (root2(nn) .lt. root1(nn))) then
rootx = root1(nn)
root1(nn) = root2(nn)
root2(nn) = rootx
endif
c.... Return real and imaginary parts of two complex roots if qq < 0.
if (qq(nn) .lt. 0.0) then
nroots(nn) = 0
root1(nn) = -0.5 * bb(nn) / (aa(nn) + fuz)
root2(nn) = 0.5 * sqrt (-qq(nn)) / (abs (aa(nn)) + fuz)
endif
c---- End of loop over coefficient sets.
130 continue
c.... See if all sets are done.
c---- Do another data set.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
go to 110
endif
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptqrtv results:' /
cbug & (i3,' nroots=',i2,' itrun=',i2 /
cbug & ' root1,root2=',1p2e22.14))
cbug write ( 3, 9903) (n, nroots(n), itrun(n),
cbug & root1(n), root2(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptqrtv. (+1 line.)
end
UCRL-WEB-209832