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