subroutine aptcubs (c0, c1, c2, tol, nroots, xroot, yroot,
& atype, itrun)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCUBS
c
c call aptcubs (c0, c1, c2, tol, nroots, xroot, yroot, atype, itrun)
c
c Version: aptcubs Updated 1997 April 17 15:40.
c aptcubs Originated 1994 December 6 10:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the real and any complex roots z = x + i*y of the cubic
c equation:
c f(z) = c0 + c1 * z + c2 * z**2 + z**3 = 0
c and to return their number (nroots), their values (root),
c and their types (atype), and a trucation indicator (itrun).
c The method uses a tolerance limit (tol) to control and measure
c the accuracy of the results. Accuracy may be reduced if the
c coefficients differ by many orders of magnitude.
c Note that for the three roots x1, x2 and x3, the following
c relationships may be used to check the accuracy of the roots:
c c0 = -z1*z2*z3.
c c1 = z1*z2 + z2*z3 + z3*z1,
c c2 = -(z1 + z2 + z3),
c Also, if one real root x1 is known, the other two roots are the
c roots of the quadratic equation:
c q(z) = c1 + x1 * (c2 + x1) + (c2 + x1) * z + z**2 = 0
c
c Method: The inversion point is translated to the origin by the
c substitution z = w - c2 / 3, which tranforms the equation in z
c to the equation in w = u + i*v:
c g(w) = w**3 + a1 * w + a0 = 0, where
c a1 = c1 - c2**2 / 3,
c a0 = (2*c2**3 - 9*c2*c1 + 27*c0) / 27.
c Note that g(0) = a0, and g'(0) = a1.
c
c The number of real and complex roots is determined by the factor
c qq = a1**3 / 27 + a0**2 / 4.
c If qq < 0, there are three unequal real roots.
c If qq = 0, there is a double real root at umax or umin, and
c another unequal real root, or three equal real roots at w = 0.
c If qq > 0, there is one real root, and two conjugate complex
c roots.
c
c If a1 = 0 and a0 = 0, there is a triple real root at w1 = 0.
c
c If only a1 = 0, there is one real root with the opposite sign
c of a0, and two complex conjugate roots:
c w1 = +/-(abs(a0))**(1/3),
c w2 = -0.5*w1*(1 - i*sqrt(3)),
c w3 = -0.5*w1*(1 + i*sqrt(3)).
c
c If only a0 = 0, there are three roots w1 = 0, w2 = -sqrt (-a1)
c and w3 = sqrt (-a1). The latter two roots are pure imaginary
c if a1 > 0.
c
c If neither a1 or a0 is zero, but a1 < 0, the equation in w has
c a maximum at umax = -sqrt (-a1 / 3) and a minimum at
c umin = sqrt (-a1 / 3). If qq = 0, either umax and umin is a
c a double real root, and a single real root is at
c w = -2 * umax or -2 * umin, respectively. The double real root
c will be tested later for truncation error, and replaced, if
c necessary, by two single real roots or two conjugate complex
c roots.
c
c If no roots have been found yet, Newton's method is used to
c iteratively find a real root, u1, using tol (or 1.e-11, if tol
c is zero) for the relative convergence criterion.
c The initial guess, w0, is as follows:
c a1 < 0, a0 < 0: w0 = wd + 2.0 * sqrt (-a1 / 3),
c a1 < 0, a0 > 0: w0 = -wd - 2.0 * sqrt (-a1 / 3),
c a1 > 0, a0 < 0: w0 = wd,
c a1 > 0, a0 > 0: w0 = -wd.
c where wd = (abs (a0))**(1/3) is the root for a1 = 0, a0 < 0.
c
c For the root w1, the subsidiary quadratic equation is:
c w**2 + w1 * w + w1**2 + a1 = 0, or
c w**2 + w1 * w - a0 / w1 = 0.
c
c If qq < 0, the two real roots of the subsidiary equation are:
c w2 = -0.5 * w1 - 0.5 * sqrt (3 * a0 / w1 - a1),
c w3 = -0.5 * w1 + 0.5 * sqrt (3 * a0 / w1 - a1).
c If w1 = 0, these become w2 = -sqrt (-a1), w3 = sqrt (-a1).
c
c If qq > 0, the two complex roots of the subsidiary equation are:
c w2 = -0.5 * w1 - 0.5 * i * sqrt (a1 - 3 * a0 / w1),
c w3 = -0.5 * w1 + 0.5 * i * sqrt (a1 - 3 * a0 / w1).
c If w1 = 0, these become w2 = -i * sqrt (a1), w3 = i * sqrt (a1).
c
c Finally, all w roots are changed to z roots by subracting
c c2 / 3, and then reordered, starting with real roots in
c ascending numerical order, followed by any pair of conjugate
c complex roots.
c Any single real roots are then iterated using Newton's method
c on the original x equation, to remove any accumulated
c truncation error.
c
c Any double roots found as described above for the case qq = 0,
c for which truncation error is excessive, are then recalculated
c by finding the roots of the subsidiary quadratic equation for
c the known single root, as above. This may replace a double root
c with two real roots or two complex roots, which are then
c reordered as above.
c
c The accuracy criterion tol is used throughout to replace
c calculated quantities with zero when they have values less than
c the estimated error in their calculation, based on tol, and
c if not zero, to determine the convergence of the Newtonian
c iteration method.
c
c Input: c0, c1, c2, tol.
c
c Output: nroots, xroot, yroot, itrun.
c
c Glossary:
c
c c2 Input Coefficient of x**2 in the cubic equation.
c
c atype Output Types of root:
c "single " Single unique real root.
c "dbl max" Double unique real root at maximum.
c "dbl min" Double unique real root at minimum.
c "triple " Triple unique real root.
c "complex" One of two conjugate complex roots.
c "unknown" Unevaluated, preceded by a double or
c triple root.
c The possible combinations of types are:
c "single ", "complex", "complex" (nroots = 1)
c "triple ", "unknown", "unknown" (nroots = 1)
c "single ", "dbl min", "unknown" (nroots = 2)
c "dbl max", "single ", "unknown" (nroots = 2)
c "single ", "single ", "single " (nroots = 3)
c The array size must be at least 3.
c
c c0 Input Coefficient of 1 in the cubic equation.
c
c c1 Input Coefficient of z in the cubic equation.
c
c itrun Output Indicates the degree of truncation error that occurred
c in finding roots of the cubic function.
c 0 if none (within machine accuracy).
c 1 if the truncation error is less than tol.
c 2 if the truncation error exceeds tol. The residual
c value of the function at the roots may be
c excessive (try increasing tol). The coefficients
c recalculated from the roots may differ excessively
c from the original coefficients. This may mean the
c magnitudes of the coefficients differ too much to
c get an accurate result.
c
c nroots Output Number of unique real roots:
c 0: The method failed. Try changing tol.
c 1: One unique real root. If atype(1) = "single",
c there are also two conjugate complex roots.
c If atype(1) = "triple".
c triple root at the inversion point of the
c function, -c2 / 3 (a1 = a0 = 0).
c 2: Two unique real roots, of which one is a double
c root at a maximum or minimum of the function
c (qq = 0, a1*a0 not zero).
c 3: Three unique real roots (qq < 0).
c
c tol Input Numerical tolerance limit and convergence control.
c If zero, a value of 1.e-11 is used for the Newtonian
c iteration method.
c If zero, itrun will almost certainly be 2.
c On 64-bit computers, or 32-bit computers with
c double precision, recommend 1.e-5 to 1.e-11.
c
c xroot Output The real parts of the unique roots of the cubic
c function. The array size must be at least 3.
c Real roots are first, in increasing numerical order,
c followed by any complex roots, in increasing
c numerical order of the imaginary part.
c If atype(1) is "triple", xroot(2) and xroot(3) will
c be -1.e99.
c If atype(1) or atype(2) is "dbl min" or "dbl max",
c xroot(3) will be -1.e99.
c
c yroot Output The imaginary parts of the unique roots of the cubic
c function. The array size must be at least 3.
c Always zero for the first root.
c Nonzero for the second and third roots, only if
c nroots = 1, and atype(1) = "single".
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c____ Roots of the cubic equation c0 + c1 * x + c2 * x**2 + x**3
dimension xroot(3) ! The real part of a root.
dimension yroot(3) ! The imaginary part of a root.
c____ Type of each root: "single", "dbl max", "dbl min" "triple", "complex".
dimension atype(3)
character*8 atype
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcubs finding roots of a cubic. tol=',1pe10.2 /
cbug & ' Equation is c0 + c1 * x + c2 * x**2 + x**3 = 0' /
cbug & ' c0,c1,c2= ',1p3e22.14)
cbug write ( 3, 9901) tol, c0, c1, c2
cbugc***DEBUG ends.
c.... Initialize.
nroots = 0
itrun = 0
xroot(1) = -1.e+99
xroot(2) = -1.e+99
xroot(3) = -1.e+99
yroot(1) = 0.0
yroot(2) = 0.0
yroot(3) = 0.0
atype(1) = 'unknown'
atype(2) = 'unknown'
atype(3) = 'unknown'
c.... See if one or more roots are zero.
if (c0 .eq. 0.0) then ! At least one root is zero.
c.... Find the other two roots.
call aptqrts (0, 1.0, c2, c1, qq, tol,
& nrootx, xnew2, xnew3, itrun)
c.... Classify and order the cubic roots.
if (nrootx .eq. 0) then ! Quadratic roots are complex.
nroots = 1
xroot(1) = 0.0
atype(1) = 'single'
xroot(2) = xnew2
yroot(2) = -xnew3
atype(2) = 'complex'
xroot(3) = xnew2
yroot(3) = xnew3
atype(3) = 'complex'
elseif (nrootx .eq. 1) then ! Quadratic root is real and double.
if (xnew2 .eq. 0.0) then
nroots = 1
xroot(1) = 0.0
atype(1) = 'triple'
elseif (xnew2 .lt. 0.0) then
nroots = 2
xroot(1) = xnew2
if ((c2 + 3.0 * xnew2) .lt. 0.0) then
atype(1) = 'dbl max'
else
atype(1) = 'dbl min'
endif
xroot(2) = 0.0
atype(2) = 'single'
else
nroots = 2
xroot(1) = 0.0
atype(1) = 'single'
xroot(2) = xnew2
if ((c2 + 3.0 * xnew2) .lt. 0.0) then
atype(2) = 'dbl max'
else
atype(2) = 'dbl min'
endif
endif
else ! Quadratic roots are real and single.
if (xnew2 .lt. 0.0) then
xroot(1) = xnew2
atype(1) = 'single'
if (xnew3 .lt. 0.0) then
nroots = 3
xroot(2) = xnew2
atype(2) = 'single'
xroot(3) = 0.0
atype(3) = 'single'
elseif (xnew3 .eq. 0.0) then
nroots = 2
xroot(2) = 0.0
if (c2 .lt. 0.0) then
atype(2) = 'dbl max'
else
atype(2) = 'dbl min'
endif
else
nroots = 3
xroot(2) = 0.0
atype(2) = 'single'
xroot(3) = xnew3
atype(3) = 'single'
endif
elseif (xnew2 .eq. 0.0) then
nroots = 2
xroot(1) = 0.0
if (c2 .lt. 0.0) then
atype(1) = 'dbl max'
else
atype(1) = 'dbl min'
endif
xroot(2) = xnew3
atype(2) = 'single'
else
nroots = 3
xroot(1) = 0.0
atype(1) = 'single'
xroot(2) = xnew2
atype(2) = 'single'
xroot(3) = xnew3
atype(3) = 'single'
endif
endif ! Tested number of real quadratic roots.
go to 700
endif ! At least one root is zero.
c.... Form the reduced equation w**3 + a1 * w + a0 = 0, by replacing
c.... z with w - c2 / 3.
a1 = c1 - c2**2 / 3.0
a1err = tol * (abs (c1) + c2**2 / 3.0)
a0 = 2.0 * c2**3 / 27.0 - c2 * c1 / 3.0 + c0
a0err = tol * (2.0 * abs (c2**3) / 27.0 + abs (c2 * c1) / 3.0 +
& abs ( c0))
cbugcbugc***DEBUG begins.
cbugcbug 9902 format (' a1, a0= ',1p2e22.14 /
cbugcbug & ' a1err,a0err=',1p2e22.14 )
cbugcbug write ( 3, 9902) a1, a0, a1err, a0err
cbugcbugc***DEBUG ends.
c.... Test for zero a1. One real root and two complex roots, or a triple root.
if (abs (a1) .le. a1err) then ! Solve w**3 + a0 = 0.
if (abs (a1) .gt. 0.0) itrun = 1
a1 = 0.0
cbugcbugc***DEBUG begins.
cbugcbug 9903 format (' Standard coefficient a1 = 0. itrun =',i2)
cbugcbug write ( 3, 9903) itrun
cbugcbugc***DEBUG ends.
nroots = 1
if (abs (a0) .le. a0err) then ! One triple root.
if (abs (a0) .gt. 0.0) itrun = 1
a0 = 0.0
xroot(1) = -c2 / 3.0
atype(1) = 'triple'
go to 700
endif
ureal = (abs (a0))**(1.0 / 3.0)
xroot(1) = -sign (ureal, a0) - c2 / 3.0
atype(1) = 'single'
xroot(2) = 0.5 * sign (ureal, a0) - c2 / 3.0
yroot(2) = 0.5 * ureal * sqrt (3.0)
atype(2) = 'complex'
xroot(3) = xroot(2)
yroot(3) = -yroot(2)
atype(3) = 'complex'
rterr = tol * (ureal + abs (c2) / 3.0)
if (abs (xroot(1)) .le. rterr) then
itrun = 1
xroot(1) = 0.0
endif
rterr = tol * (0.5 * ureal + abs (c2) / 3.0)
if (abs (xroot(2)) .le. rterr) then
itrun = 1
xroot(2) = 0.0
xroot(3) = 0.0
endif
go to 700
endif ! Solved y**3 + a0 = 0.
c.... Test for zero a0.
if (abs (a0) .le. a0err) then ! Solve y * (y**2 + a1) = 0.
if (abs (a0) .gt. 0.0) itrun = 1
a0 = 0.0
cbugcbugc***DEBUG begins.
cbugcbug 9904 format (' Standard coefficient a0 = 0. itrun =',i2)
cbugcbug write ( 3, 9904) itrun
cbugcbugc***DEBUG ends.
if (a1 .lt. 0.0) then ! Three unequal real roots.
nroots = 3
ureal = sqrt (-a1)
xroot(1) = -c2 / 3.0 - ureal
xroot(2) = -c2 / 3.0
xroot(3) = -c2 / 3.0 + ureal
atype(1) = 'single'
atype(2) = 'single'
atype(3) = 'single'
rterr = tol * (abs (ureal) + abs (c2) / 3.0)
if (abs (xroot(1)) .le. rterr) then
xroot(1) = 0.0
endif
if (abs (xroot(3)) .le. rterr) then
xroot(3) = 0.0
endif
go to 300
else ! One real root and two complex roots.
nroots = 1
xroot(1) = -c2 / 3.0
atype(1) = 'single'
vroot = sqrt (a1)
xroot(2) = -c2 / 3.0
yroot(2) = -vroot
atype(2) = 'complex'
xroot(3) = -c2 / 3.0
yroot(3) = vroot
atype(3) = 'complex'
go to 700
endif
endif ! Solved y * (y**2 + a1) = 0.
c.... Find the factor determining the number of roots.
qq = a1**3 / 27.0 + 0.25 * a0**2
qqerr = tol * (abs (a1**3) / 27.0 + 0.25 * a0**2) +
& a1**2 * a1err / 9.0 + 0.50 * abs (a0) * a0err
cbugcbugc***DEBUG begins.
cbugcbug 9908 format (' Standard coefficient qq = ',1p2e22.14)
cbugcbug write ( 3, 9908) qq, qqerr
cbugcbugc***DEBUG ends.
if (abs (qq) .le. qqerr) then
if (abs (qq) .gt. 0.0) itrun = 1
cbugcbugc***DEBUG begins.
cbugcbug 9909 format (' Standard coefficient qq = 0. itrun =',i2)
cbugcbug write ( 3, 9909) itrun
cbugcbugc***DEBUG ends.
qq = 0.0
endif
c.... Find y coordinates of any minimum and maximum.
nextr = 0
if (a1 .lt. 0.0) then ! A minimum and a maximum exist.
nextr = 2
umin = sqrt (-a1 / 3.0)
fmin = umin**3 + a1 * umin + a0
fminerr = tol * (abs (umin**3) + abs (a1 * umin) + abs(a0)) +
& a1err * abs (umin) + a0err
umax = -sqrt (-a1 / 3.0)
fmax = umax**3 + a1 * umax + a0
fmaxerr = tol * (abs (umax**3) + abs (a1 * umax) + abs(a0)) +
& a1err * abs (umax) + a0err
cbugcbugc***DEBUG begins.
cbugcbug 9905 format (' umin, umax =',1p2e22.14 /
cbugcbug & ' fmin, fmax =',1p2e22.14 /
cbugcbug & ' errmn,errmx=',1p2e22.14 )
cbugcbug write ( 3, 9905) umin, umax, fmin, fmax, fminerr, fmaxerr
cbugcbugc***DEBUG ends.
if ((qq .eq. 0.0) .and. (abs (fmin) .lt. abs (fmax)) .or.
& (abs(fmin) .le. fminerr)) then ! Double root at umin.
if (abs (fmin) .gt. 0.0) itrun = 1
if (abs (fmin) .gt. fminerr) itrun = 2
cbugcbugc***DEBUG begins.
cbugcbug 9906 format (' Double root at umin. itrun =',i2)
cbugcbug write ( 3, 9906) itrun
cbugcbugc***DEBUG ends.
nroots = 2 ! Double root at umin.
xroot(1) = -2.0 * umin - c2 / 3.0
atype(1) = 'single'
rterr = tol * (2.0 * abs (umin) + abs (c2) / 3.0)
if (abs (xroot(1)) .le. rterr) then
xroot(1) = 0.0
endif
xroot(2) = umin - c2 / 3.0
atype(2) = 'dbl min'
rterr = tol * (abs (umin) + abs (c2) / 3.0)
if (abs (xroot(2)) .le. rterr) then
xroot(2) = 0.0
endif
go to 300
endif
if ((qq .eq. 0.0) .and. (abs (fmax) .le. abs (fmin)) .or.
& (abs(fmax) .le. fmaxerr)) then ! Double root at umax.
if (abs (fmax) .gt. 0.0) itrun = 1
if (abs (fmax) .gt. fmaxerr) itrun = 2
cbugcbugc***DEBUG begins.
cbugcbug 9907 format (' Double root at umax. itrun =',i2)
cbugcbug write ( 3, 9907) itrun
cbugcbugc***DEBUG ends.
nroots = 2 ! Double root at umax.
xroot(1) = umax - c2 / 3.0
atype(1) = 'dbl max'
rterr = tol * (abs (umax) + abs (c2) / 3.0)
if (abs (xroot(1)) .le. rterr) then
xroot(1) = 0.0
endif
xroot(2) = -2.0 * umax - c2 / 3.0
atype(2) = 'single'
rterr = tol * (2.0 * abs (umax) + abs (c2) / 3.0)
if (abs (xroot(2)) .le. rterr) then
xroot(2) = 0.0
endif
go to 300
endif
endif ! Tested a1.
c.... Find number of real roots.
if (qq .lt. 0.0) then ! Three real and unequal roots.
cbugcbugc***DEBUG begins.
cbugcbug 9910 format (' Three real and unequal roots.')
cbugcbug write ( 3, 9910)
cbugcbugc***DEBUG ends.
nroots = 3
elseif (qq .eq. 0.0) then ! IMPOSSIBLE. Already done above.
cbugcbugc***DEBUG begins.
cbugcbug 9911 format (' IMPOSSIBLE LOGIC. qq = 0 again.')
cbugcbug write ( 3, 9911)
cbugcbugc***DEBUG ends.
nroots = 3
go to 700
elseif (qq .gt. 0.0) then ! One real root.
cbugcbugc***DEBUG begins.
cbugcbug 9912 format (' One real root.')
cbugcbug write ( 3, 9912)
cbugcbugc***DEBUG ends.
nroots = 1
endif
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " Iterate to find a single real root.")')
cbugcbugc***DEBUG ends.
c.... Iterate, using Newton's method, to find a single real root.
c.... Choose initial guess for Newton's method.
ureal = (abs (a0))**(1.0 / 3.0)
if (a1 .gt. 0.0) then
utest = 0.0
if (a0 .lt. 0.0) utest = ureal
if (a0 .gt. 0.0) utest = -ureal
elseif (a0 .lt. 0.0) then
utest = 2.0 * umin + ureal
else
utest = 2.0 * umax - ureal
endif
c.... Iterate, using Newton's method.
tolx = tol
if (tol .eq. 0.0) tolx = 1.e-11
nstop = 0
niter = 0
cbugcbugc***DEBUG begins.
cbugcbug 9922 format (' It',i4,' y,f,s',1p3e22.14 )
cbugcbug ftest = utest**3 + a1 * utest + a0
cbugcbug fslope = 3.0 * utest**2 + a1
cbugcbug write ( 3, 9922) niter, utest, ftest, fslope
cbugcbugc***DEBUG ends.
110 niter = niter + 1
ftest = utest**3 + a1 * utest + a0
fterr = tolx * (abs (utest**3) + abs (a1 * utest) + abs (a0))
fslope = 3.0 * utest**2 + a1
fslerr = tolx * (3.0 * utest**2 + abs (a1))
if (abs (fslope) .le. fslerr) then
unew1 = utest
itrun = 2
go to 200
else
unew1 = utest - ftest / fslope
endif
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9922) niter, unew1, ftest, fslope
cbugcbugc***DEBUG ends.
if (nstop .eq. 1) go to 200
if (abs (unew1 - utest) .le. tolx * abs (utest)) nstop = 1
if (abs (ftest) .lt. fterr) nstop = 1 ! TRY
if (niter .gt. 1000) then
itrun = 2
go to 200
endif
utest = unew1
go to 110
200 xroot(1) = unew1 - c2 / 3.0
rterr = tol * (abs (unew1) + abs (c2) / 3.0)
if (abs (xroot(1)) .le. rterr) then
xroot(1) = 0.0
endif
atype(1) = 'single'
cbugcbugc***DEBUG begins.
cbugcbug 9914 format (' Root y, x= ',1p2e22.14)
cbugcbug write ( 3, 9914) unew1, xroot(1)
cbugcbugc***DEBUG ends.
if (nroots .eq. 1) go to 300
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " Solve quadratic for other two roots.")')
cbugcbugc***DEBUG ends.
c.... Find other two roots by factoring out single root, solving quadratic.
g2 = 1.0
g1 = unew1
g0 = -a0 / unew1
call aptqrts (0, g2, g1, g0, qq, tol,
& nrootx, unew2, unew3, itrun)
if (nrootx .ne. 2) then
cbugcbugc***DEBUG begins.
cbugcbug 9918 format (' IMPOSSIBLE LOGIC. aptqrts says not 2 roots.')
cbugcbug write ( 3, 9918)
cbugcbugc***DEBUG ends.
nroots = 0
go to 700
endif
if (unew1 .lt. unew2) then
xroot(1) = unew1 - c2 / 3.0
rterr = tol * (abs (unew1) + abs (c2) / 3.0)
if (abs (xroot(1)) .le. rterr) then
xroot(1) = 0.0
endif
xroot(2) = unew2 - c2 / 3.0
rterr = tol * (abs (unew2) + abs (c2) / 3.0)
if (abs (xroot(2)) .le. rterr) then
xroot(2) = 0.0
endif
xroot(3) = unew3 - c2 / 3.0
rterr = tol * (abs (unew3) + abs (c2) / 3.0)
if (abs (xroot(3)) .le. rterr) then
xroot(3) = 0.0
endif
elseif (unew1 .gt. unew3) then
xroot(1) = unew2 - c2 / 3.0
rterr = tol * (abs (unew2) + abs (c2) / 3.0)
if (abs (xroot(1)) .le. rterr) then
xroot(1) = 0.0
endif
xroot(2) = unew3 - c2 / 3.0
rterr = tol * (abs (unew3) + abs (c2) / 3.0)
if (abs (xroot(2)) .le. rterr) then
xroot(2) = 0.0
endif
xroot(3) = unew1 - c2 / 3.0
rterr = tol * (abs (unew1) + abs (c2) / 3.0)
if (abs (xroot(3)) .le. rterr) then
xroot(3) = 0.0
endif
endif
atype(1) = 'single'
atype(2) = 'single'
atype(3) = 'single'
c.... Use Newtonian iteration to improve real roots not at umax or umin.
300 continue
tolx = tol
if (tol .eq. 0.0) tolx = 1.e-11
do 320 n = 1, nroots
if ((atype(n) .eq. 'dbl max') .or.
& (atype(n) .eq. 'dbl min')) then
ftest = xroot(n)**3 + c2 * xroot(n)**2 + c1 * xroot(n) + c0
fterr = tolx * (abs (xroot(n)**3) + abs (c2) * xroot(n)**2 +
& abs (c1 * xroot(n)) + abs (c0))
if (abs(ftest) .gt. fterr) then
itrun = 2
cbugcbugc***DEBUG begins.
cbugcbug 9931 format (/ ' Set itrun = 2. Bad root ',i2,2x,a8 /
cbugcbug & ' ftest,fterr=',1p2e22.14)
cbugcbug write ( 3, 9931) n, atype(n), ftest, fterr
cbugcbugc***DEBUG ends.
endif
go to 320
endif
xtest = xroot(n)
nstop = 0
niter = 0
cbugcbugc***DEBUG begins.
cbugcbug 9923 format (' It',i4,' x,f,s',1p3e22.14 )
cbugcbug write ( 3, '(/ " Iterate to improve this single root.")')
cbugcbug ftest = xtest**3 + c2 * xtest**2 + c1 * xtest + c0
cbugcbug fslope = 3.0 * xtest**2 + 2.0 * c2 * xtest + c1
cbugcbug write ( 3, 9923) niter, xtest, ftest, fslope
cbugcbugc***DEBUG ends.
310 niter = niter + 1
ftest = xtest**3 + c2 * xtest**2 + c1 * xtest + c0
fterr = tolx * (abs (xtest**3) + abs (c2 * xtest) +
& abs (c1 * xtest) + abs (c0))
fslope = 3.0 * xtest**2 + 2.0 * c2 * xtest + c1
fslerr = tolx * (3.0 * xtest**2 + 2.0 * abs (c2 * xtest) +
& abs (c1))
if (abs (fslope) .le. fslerr) then
itrun = 2
xroot(n) = xtest
go to 320
else
xroot(n) = xtest - ftest / fslope
endif
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9923) niter, xroot(n), ftest, fslope
cbugcbugc***DEBUG ends.
if (nstop .eq. 1) then
ftest = xroot(n)**3 + c2 * xroot(n)**2 + c1 * xroot(n) + c0
fterr = tolx * (abs (xroot(n)**3) + abs (c2) * xroot(n)**2 +
& abs (c1 * xroot(n)) + abs (c0))
if (abs(ftest) .gt. fterr) then
itrun = 2
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9931) n, atype(n), ftest, fterr
cbugcbugc***DEBUG ends.
endif
go to 320
endif
if (abs (xroot(n) - xtest) .le. tolx * abs (xtest)) nstop = 1
if (abs (ftest) .lt. fterr) nstop = 1 ! TRY.
if (niter .gt. 1000) then
itrun = 2
go to 320
endif
xtest = xroot(n)
go to 310
320 continue ! End of loop over roots.
c.... Find any complex roots.
if ((nroots .eq. 1) .and. (atype(1) .eq. 'single')) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " Find two complex roots.")')
cbugcbugc***DEBUG ends.
rreal = -0.5 * unew1 - c2 / 3.0
errx = tol * (0.5 * abs(unew1) + abs (c2) / 3.0)
if (abs (rreal) .le. errx) then
if (itrun .eq. 0) itrun = 1
rreal = 0.0
endif
cbugcbugc***DEBUG begins.
cbugcbug 9933 format (' a1 - 3*a0/u1 =',1pe22.14 /
cbugcbug & ' -u1**2-4*a0/u1=',1pe22.14)
cbugcbug arg1 = a1 - 3.0 * a0 / unew1
cbugcbug arg2 = -(unew1**2 + 4.0 * a0 / unew1)
cbugcbug write ( 3, 9933) arg1, arg2
cbugcbug if ((arg1 .lt. 0.0) .or. (arg2 .lt. 0.0)) then
cbugcbug nroots = 0
cbugcbug go to 700
cbugcbug endif
cbugcbugc***DEBUG ends.
rimag = 0.5 * sqrt (a1 - 3.0 * a0 / unew1)
xroot(2) = rreal
yroot(2) = -rimag
atype(2) = 'complex'
xroot(3) = rreal
yroot(3) = rimag
atype(3) = 'complex'
endif ! Found two complex roots.
c.... Replace any badly truncated double roots.
if ((nroots .eq. 2) .and. (itrun .eq. 2)) then
if (atype(1) .eq. 'single') then ! Replace double root 2.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " Replace double root 2. Solve quadratic.")')
cbugcbugc***DEBUG ends.
g0 = c1 + xroot(1) * (c2 + xroot(1))
g0err = tol * (abs (c1) + abs (c2 * xroot(1)) + xroot(1)**2)
if (xroot(1) .ne. 0.0) then
g0err2 = tol * abs (c0 / xroot(1))
if (g0err2 .le. g0err) then
g0 = -c0 / xroot(1)
g0err = g0err2
endif
endif
if (abs (g0) .le. g0err) g0 = 0.0
g1 = c2 + xroot(1)
g1err = tol * (abs (c2) + abs (xroot(1)))
if (xroot(1) .ne. 0.0) then
g1err2 = tol *
& (abs (c0) + abs (c1 * xroot(1)) / xroot(1)**2)
if (g1err2 .le. g1err) then
g1 = -(c0 + c1 * xroot(1)) / xroot(1)**2
g1err = g1err2
endif
endif
if (abs (g1) .le. g1err) g1 = 0.0
g2 = 1.0
call aptqrts (0, g2, g1, g0, qq, tol,
& nrootx, xnew2, xnew3, itrunx)
if (nrootx .eq. 2) then
if (itrun .eq. 0) itrun = 1
nroots = 3
xroot(2) = xnew2
atype(2) = 'single'
xroot(3) = xnew3
atype(3) = 'single'
elseif (nrootx .eq. 0) then
if (itrun .eq. 0) itrun = 1
nroots = 1
q = sqrt (4.0 * g2 * g0 - g1**2)
rreal = -0.5 * g1 / g2
rimag = 0.5 * q / g2
xroot(2) = rreal
yroot(2) = -rimag
atype(2) = 'complex'
xroot(3) = rreal
yroot(3) = rimag
atype(3) = 'complex'
endif
else ! Replace double root 1.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " Replace double root 1. Solve quadratic.")')
cbugcbugc***DEBUG ends.
g0 = c1 + xroot(2) * (c2 + xroot(2))
g0err = tol * (abs (c1) + abs (c2 * xroot(2)) + xroot(2)**2)
if (xroot(2) .ne. 0.0) then
g0err2 = tol * abs (c0 / xroot(2))
if (g0err2 .le. g0err) then
g0 = -c0 / xroot(2)
g0err = g0err2
endif
endif
if (abs (g0) .le. g0err) g0 = 0.0
g1 = c2 + xroot(2)
g1err = tol * (abs (c2) + abs (xroot(2)))
if (xroot(2) .ne. 0.0) then
g1err2 = tol *
& (abs (c0) + abs (c1 * xroot(2)) / xroot(2)**2)
if (g1err2 .le. g1err) then
g1 = -(c0 + c1 * xroot(2)) / xroot(2)**2
g1err = g1err2
endif
endif
if (abs (g1) .le. g1err) g1 = 0.0
g2 = 1.0
call aptqrts (0, g2, g1, g0, qq, tol,
& nrootx, xnew1, xnew2, itrunx)
if (nrootx .eq. 2) then
if (itrun .eq. 0) itrun = 1
nroots = 3
xroot(3) = xroot(2)
atype(3) = 'single'
xroot(1) = xnew1
atype(1) = 'single'
xroot(2) = xnew2
atype(2) = 'single'
elseif (nrootx .eq. 0) then
if (itrun .eq. 0) itrun = 1
nroots = 1
q = sqrt (4.0 * g2 * g0 - g1**2)
rreal = -0.5 * g1 / g2
rimag = 0.5 * q / g2
xroot(1) = xroot(2)
atype(1) = 'single'
xroot(2) = rreal
yroot(2) = -rimag
atype(2) = 'complex'
xroot(3) = rreal
yroot(3) = rimag
atype(3) = 'complex'
endif
endif
endif
700 continue
c.... See if a single real root at the origin can be made more accurate.
if ((atype(1) .eq. 'single') .and.
& (xroot(1) .eq. 0.0)) then
xroot(1) = -c0 / (xroot(2)*xroot(3) + yroot(2) * yroot(3))
elseif ((atype(2) .eq. 'single') .and.
& (xroot(2) .eq. 0.0)) then
xroot(2) = -c0 / (xroot(3)*xroot(1) + yroot(3) * yroot(1))
elseif ((atype(3) .eq. 'single') .and.
& (xroot(3) .eq. 0.0)) then
xroot(3) = -c0 / (xroot(1)*xroot(2) + yroot(1) * yroot(2))
endif
c.... Check the value of the function at the roots.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " Check the function at the roots.")')
cbugcbugc***DEBUG ends.
do 710 n = 1, 3 ! Loop over roots.
if (atype(n) .eq. 'unknown') go to 710
frres = xroot(n)**3 + c2 * xroot(n)**2 + c1 * xroot(n) + c0 -
& yroot(n)**2 * (3.0 * xroot(n) + c2)
fires = yroot(n) * (3.0 * xroot(n)**2 - yroot(n)**2 +
& 2.0 * c2 * xroot(n) + c1)
frerr = tol * (abs (xroot(n)**3) + abs (c2) * xroot(n)**2 +
& abs (c1 * xroot(n)) + abs (c0) +
& yroot(n)**2 * (3.0 * abs (xroot(n)) + abs (c2)))
fierr = tol * (abs (yroot(n)) * (3.0 * xroot(n)**2) +
& yroot(n)**2 + 2.0 * abs (c2 * xroot(n)) +
& abs (c1))
if (abs (frres) .gt. frerr) itrun = 2
if (abs (fires) .gt. fierr) itrun = 2
710 continue ! End of loop over roots.
c.... Reconstruct the coefficients.
kbad = 0
if (nroots .eq. 1) then
if (atype(1) .eq. 'single') then
c0p = -xroot(1) * (xroot(2) * xroot(3) + yroot(2)**2)
c1p = xroot(1) * (xroot(2) + xroot(3)) +
& xroot(2) * xroot(3) + yroot(2)**2
c2p = -(xroot(1) + xroot(2) + xroot(3))
else
c0p = -xroot(1)**3
c1p = 3.0 * xroot(1)**2
c2p = -3.0 * xroot(1)
endif
elseif (nroots .eq. 2) then
if (atype(1) .eq. 'single') then
c0p = -xroot(1) * xroot(2)**2
c1p = 2.0 * xroot(1) * xroot(2) + xroot(2)**2
c2p = -(xroot(1) + 2.0 * xroot(2))
else
c0p = -xroot(2) * xroot(1)**2
c1p = 2.0 * xroot(1) * xroot(2) + xroot(1)**2
c2p = -(xroot(2) + 2.0 * xroot(1))
endif
elseif (nroots .eq. 3) then
c0p = -xroot(1) * xroot(2) * xroot(3)
c1p = xroot(1) * (xroot(2) + xroot(3)) + xroot(2) * xroot(3)
c2p = -(xroot(1) + xroot(2) + xroot(3))
endif
if (abs (c0 - c0p) .gt. tol * (abs (c0) + abs (c0p))) kbad = 1
if (abs (c1 - c1p) .gt. tol * (abs (c1) + abs (c1p))) kbad = 1
if (abs (c2 - c2p) .gt. tol * (abs (c2) + abs (c2p))) kbad = 1
cbugcbugc***DEBUG begins.
cbugcbug 9930 format (' c0,c1,c2(in)',1p3e22.14 /
cbugcbug & ' c0,c1,c2(re)',1p3e22.14 )
cbugcbug
cbugcbug write ( 3, 9930) c0, c1, c2, c0p, c1p, c2p
cbugcbugc***DEBUG ends.
cbugc***DEBUG begins.
cbug 9920 format (/ 'aptcubs results. nroots=',i2)
cbug write ( 3, 9920) nroots
cbugc***DEBUG ends.
do 990 n = 1, 3
if (atype(n) .eq. 'unknown') go to 990
frres = xroot(n)**3 + c2 * xroot(n)**2 + c1 * xroot(n) + c0 -
& yroot(n)**2 * (3.0 * xroot(n) + c2)
frerr = tol * (abs (xroot(n)**3) + abs (c2) * xroot(n)**2 +
& abs (c1 * xroot(n)) + abs (c0) +
& yroot(n)**2 * (3.0 * abs (xroot(n)) + abs (c2)))
fires = yroot(n) * (3.0 * xroot(n)**2 - yroot(n)**2 +
& 2.0 * c2 * xroot(n) + c1)
fierr = tol * (abs (yroot(n)) * (3.0 * xroot(n)**2) +
& yroot(n)**2 + 2.0 * abs (c2 * xroot(n)) +
& abs (c1))
if (abs (frres) .gt. frerr) kbad = 1
if (abs (fires) .gt. fierr) kbad = 1
cbugc***DEBUG begins.
cbug 9921 format (' Root ',i1,' r,i=',1p2e22.14, 2x,a8 /
cbug & ' cubic r,i=',1p2e22.14 /
cbug & ' allow r,i=',1p2e22.14 )
cbug write ( 3, 9921) n, xroot(n), yroot(n), atype(n), frres, fires,
cbug & frerr, fierr
cbugc***DEBUG ends.
990 continue
if ((kbad .eq. 0) .and. (itrun .eq. 2)) itrun = 1
if ((kbad .eq. 1) .and. (itrun .eq. 0)) itrun = 1
cbugc***DEBUG begins.
cbug 9940 format (/ ' itrun=',i2 )
cbug write ( 3, 9940) itrun
cbugc***DEBUG ends.
return
c.... End of subroutine aptcubs. (+1 line.)
end
UCRL-WEB-209832