subroutine aptquar (c0, c1, c2, c3, tol, nroot, nrooti,
& xroot, yroot, atype, itrun)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQUAR
c
c call aptquar (c0, c1, c2, c3, tol, nroot, nrooti,
c xroot, yroot, atype, itrun)
c
c Version: aptquar Updated 1997 April 4 17:00.
c aptquar Originated 1997 March 14 18: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
c quartic equation:
c f(z) = c0 + c1 * z + c2 * z**2 + c3 * z**3 + z**4 = 0
c and to return their number (nroot and nrooti), their values
c (xroot and yroot), and their types (atype), and a trucation
c 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 Any complex roots occur in conjugate pairs, so the quartic
c equation may have only 0, 2 or 4 complex roots.
c
c Note that for the four roots z1, z2, z3 and z4, the following
c relationships may be used to check the accuracy of the roots:
c c0 = z1*z2*z3*z4,
c c1 = -(z1*z2*z3 + z2*z3*z4 + z3*z4*z1 + z4*z1*z2),
c c2 = z1*z2 + z2*z3 + z3*z4 + z4*z1 + z1*z3 + z2*z4,
c c3 = -(z1 + z2 + z3 + z4).
c
c Also, if one real root x1 is known, the other three roots are
c the roots of the cubic equation:
c q(z) = d0 + d1 * z + d2 * z**2 + z**3,
c where d0 = c1 + c2 * x1 + c3 * x1**2 + x1**3
c = -c0 / x1 (x1 nonzero),
c d1 = c2 + c3 * x1 + x1**2
c = -(c0 + c1 * x1) / x1**2 (x1 nonzero),
c d2 = c3 + x1
c = -(c0 + c1 * x1 + c2 * x1**2) / x1**3 (x1 nonzero)
c (for each coefficient, use the expression which minimizes the
c numerical truncation error)
c Either one or three of the roots of the cubic equation are real.
c If three are real, one may be a double or triple root.
c
c If a pair of roots z1 and z2 are known, the other two roots are
c the roots of the quadratic equation:
c q(z) = d0 + d1 * z + z**2,
c where d0 = c0 / (z1 * z2),
c and d1 = c3 + Z1 + z2.
c There are either two real roots or a pair of conjugate complex
c roots. The two real roots may be the same.
c
c If significant truncation error occurs (itrun = 2), try calling
c aptquar with the coefficients c0' = c3 / c0, c1' = c2 / c0,
c c2' = c1 / c0, c3' = 1 / c0, and use the reciprocals of the
c resulting roots. This may make very small roots more accurate.
c It may not be possible to get accurate values of complex roots
c for which the real and imaginary parts differ by many orders of
c magnitude.
c
c Method: The method in the "Handbook of tables for Mathematics", 3rd Ed.,
c 1967, page 133 is used.
c
c Any real root of maximum magnitude is then factored out of the
c quartic equation, and the resulting cubic equation is solved.
c
c Any zero roots are recalculated to eliminate as much of the
c truncation error as possible.
c
c Finally, all roots are reordered, starting with real roots in
c ascending numerical order, followed by any pair of conjugate
c complex roots.
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.
c
c Input: c0, c1, c2, c3, tol.
c
c Output: nroot, nrooti, xroot, yroot, itrun.
c
c Glossary:
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 "quadrpl" Quadruple unique real root.
c "complex" One of two or four unique complex roots.
c "dbl com" One of two double unique complex roots.
c "unknown" Unevaluated.
c The array size must be at least 4.
c
c c0 Input Coefficient of 1 in the quartic equation.
c
c c1 Input Coefficient of z in the quartic equation.
c
c c2 Input Coefficient of z**2 in the quartic equation.
c
c c3 Input Coefficient of z**3 in the quartic equation.
c
c itrun Output Indicates the degree of truncation error that occurred
c in finding roots of the quartic 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 nroot Output Number of unique real roots. Possible values:
c 0 (nrooti = 0):
c The method failed. Try changing tol.
c 0 (nrooti = 2)
c "dbl com", "dbl com"
c 0 (nrooti = 4):
c "complex", "complex", "complex", "complex".
c 1 (nrooti = 0):
c "quadrpl" (at the minimum point).
c 1 (nrooti = 2):
c "dbl min", "complex", "complex".
c 2 (nrooti = 0):
c "single", "triple" (at an inversion point), or
c "dbl min", "dbl min".
c 2 (nrooti = 2):
c "single", "single", "complex", "complex".
c 3 (nrooti = 0):
c "dbl max", "single", "single", or
c "dbl min", "single", "single".
c 4 (nrooti = 0):
c "single", "single", "single", "single".
c
c nrooti Output The number of complex roots, which must be either 0, 2
c or 4, as they only occur in conjugate complex pairs.
c See nroot.
c
c tol Input Numerical tolerance limit and convergence control.
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 quartic
c function. The array size must be at least 4.
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 May truncate to zero if yroot is many orders of
c magnitude greater than xroot.
c
c yroot Output The imaginary parts of the complex roots of the quartic
c function. The array size must be at least 4.
c Always zero for the first nroot roots.
c Nonzero for roots (nroot + 1) through
c (nroot + nrooti).
c May truncate to zero if xroot is many orders of
c magnitude greater than yroot.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c.... Roots of the quartic equation z**4 + c3*z**3 + c2*z**2 + c1*z + c0 = 0.
dimension xroot(4) ! The real part of a root.
dimension yroot(4) ! The imaginary part of a root.
c.... Quadric root types:
c.... "single", "dbl max", "dbl min", "triple", "quadrpl", "complex".
dimension atype(4) ! Root type for quartic.
character*8 atype ! Root type for quartic.
c.... Cubic root types:
c.... "single", "dbl max", "dbl min", "triple", "complex".
dimension atcub(3) ! Root type for cubic.
character*8 atcub ! Root type for cubic.
dimension xcubr(3) ! The real part of a root of a cubic.
dimension xcubi(3) ! The imaginary part of a root of a cubic.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptquar finding roots of a quartic. tol=',1pe10.2 /
cbug & ' Equation is c0 + c1 * z + c2 * z**2 + c3 * z**2 + z**4 = 0' /
cbug & ' c0,c1,c2= ',1p3e22.14 /
cbug & ' c3= ',1pe22.14)
cbug write ( 3, 9901) tol, c0, c1, c2, c3
cbugc***DEBUG ends.
c.... Initialize.
nroot = 0
nrooti = 0
itrun = 0
xroot(1) = -1.e+99
xroot(2) = -1.e+99
xroot(3) = -1.e+99
xroot(4) = -1.e+99
yroot(1) = 0.0
yroot(2) = 0.0
yroot(3) = 0.0
yroot(4) = 0.0
atype(1) = 'unknown'
atype(2) = 'unknown'
atype(3) = 'unknown'
atype(4) = 'unknown'
c.... See if one of the roots is zero.
if (c0 .eq. 0.0) then ! At least one root is zero.
xroot(1) = 0.0
yroot(1) = 0.0
call aptcubs (c1, c2, c3, tol, nroots, xroot(2), yroot(2),
& atype(2), itrun)
if (nroots .eq. 0) then
kbad = 1
itrun = 2
go to 395
elseif (atype(2) .eq. 'triple') then
xroot(3) = xroot(2)
yroot(3) = yroot(2)
xroot(4) = xroot(2)
yroot(4) = yroot(2)
elseif ((atype(2) .eq. 'dbl min') .or.
& (atype(2) .eq. 'dbl max')) then
xroot(4) = xroot(2)
yroot(4) = yroot(2)
elseif ((atype(3) .eq. 'dbl min') .or.
& (atype(3) .eq. 'dbl max')) then
xroot(4) = xroot(3)
yroot(4) = yroot(3)
endif
go to 200
endif ! At least one root is zero.
c.... Find the coefficients of the resolvent cubic equation.
c.... f(u) = a0 + a1 * u + a2 * u**2 + u**3
c.... = (-c3**2 * c0 + 4 * c2 * c0 - c**2) +
c.... (c3 * c1 - 4 * c0) * u - c2 * u**2 + u**3 = 0
a0 = -c3**2 * c0 + 4.0 * c2 * c0 - c1**2
a0err = tol * (c3**2 * abs (c0) + 4.0 * abs (c2 * c0) + c1**2)
if (abs (a0) .le. a0err) then
if ((itrun .eq. 0) .and. (a0 .ne. 0.0)) itrun = 1
a0 = 0.0
endif
a1 = c3 * c1 - 4.0 * c0
a1err = tol * (abs (c3 * c1) + 4.0 * abs (c0))
if (abs (a1) .le. a1err) then
if ((itrun .eq. 0) .and. (a1 .ne. 0.0)) itrun = 1
a1 = 0.0
endif
a2 = -c2
c.... Find the most positive real root of the resolvent cubic equation.
call aptcubs (a0, a1, a2, tol,
& nroots, xcubr, xcubi, atcub, itrun)
u1 = xcubr(nroots)
cbugc***DEBUG begins.
cbug 9701 format ('aptquar DEBUG. Resolvent root u1 = ',1pe22.14)
cbug write ( 3, 9701) u1
cbugc***DEBUG ends.
c.... Find R = sqrt (rsq) = rreal + i * rimag, where i = sqrt (-1),
c.... and rsq = c3**2 / 4 - c2 + u1.
c.... 1 / R = rreal / abs (rsq) + i * rimag / abs (rsq).
rsq = 0.25 * c3**2 - c2 + u1
rsqerr = tol * (0.25 * c3**2 + abs (c2) + abs (u1))
if (abs (rsq) .le. rsqerr) then
if ((itrun .eq. 0) .and. (rsq .ne. 0.0)) itrun = 1
rsq = 0.0
endif
if (rsq .lt. 0.0) then
rreal = 0.0
rimag = sqrt (-rsq)
else
rreal = sqrt (rsq)
rimag = 0.0
endif
cbugc***DEBUG begins.
cbug 9702 format ('aptquar DEBUG. R = ',1pe22.14,' + i *',1pe22.14)
cbug write ( 3, 9702) rreal, rimag
cbugc***DEBUG ends.
rmag = sqrt (abs (rsq))
c.... Find dcap = sqrt (dsq) = sqrt (ff + gg) and
c.... ecap = sqrt (esq) = sqrt (ff - gg),
c.... where ff = 0.75 * c3**2 - 2.0 * c2 - rsq, and
c.... and gg = 2 * sqrt (u1**2 - 4 * c0) (R = 0),
c.... or gg = (c3 * c2 - 2 * c1 - 0.25 * c3**3) / R (R nonzero).
c.... gg = ggreal + i * ggimag
ff = 0.75 * c3**2 - 2.0 * c2 - rsq
fferr = tol * (0.75 * c3**2 + 2.0 * abs (c2) + abs (rsq))
if (abs (ff) .le. fferr) then
if ((itrun .eq. 0) .and. (ff .ne. 0.0)) itrun = 1
ff = 0.0
endif
cbugc***DEBUG begins.
cbug 9703 format ('aptquar DEBUG. ff = ',1pe22.14)
cbug write ( 3, 9703) ff
cbugc***DEBUG ends.
if (rmag .eq. 0.0) then
ff = 0.75 * c3**2 - 2.0 * c2
fferr = tol * (0.75 * c3**2 + 2.0 * abs (c2))
if (abs (ff) .le. fferr) then
if ((itrun .eq. 0) .and. (ff .ne. 0.0)) itrun = 1
ff = 0.0
endif
ggarg = u1**2 - 4.0 * c0
ggargerr = tol * (u1**2 + 4.0 * abs (c0))
if (abs (ggarg) .le. ggargerr) then
if ((itrun .eq. 0) .and. (ggarg .ne. 0.0)) itrun = 1
ggarg = 0.0
endif
if (ggarg .lt. 0.0) then
ggreal = 0.0
ggimag = 2.0 * sqrt (-ggarg)
else
ggreal = 2.0 * sqrt (ggarg)
ggimag = 0.0
endif
else
ff = 0.5 * c3**2 - c2 - u1
fferr = tol * (0.5 * c3**2 + abs (c2) + abs (u1))
if (abs (ff) .le. fferr) then
if ((itrun .eq. 0) .and. (ff .ne. 0.0)) itrun = 1
ff = 0.0
endif
fact = c3 * c2 - 2.0 * c1 - 0.25 * c3**3
facterr = tol * (abs (c3 * c2) + 2.0 * abs (c1) +
& 0.25 * abs (c3**3))
if (abs (fact) .le. facterr) then
if ((itrun .eq. 0) .and. (fact .ne. 0.0)) itrun = 1
fact = 0.0
endif
ggreal = fact * rreal / abs (rsq)
ggimag = -fact * rimag / abs (rsq)
endif
cbugc***DEBUG begins.
cbug 9704 format ('aptquar DEBUG. gg = ',1pe22.14,' + i *',1pe22.14)
cbug write ( 3, 9704) ggreal, ggimag
cbugc***DEBUG ends.
dsqreal = ff + ggreal
dsqerr = tol * (abs (ff) + abs (ggreal))
if (abs (dsqreal) .le. dsqerr) then
if ((itrun .eq. 0) .and. (dsqreal .ne. 0.0)) itrun = 1
dsqreal = 0.0
endif
dsqimag = ggimag
esqreal = ff - ggreal
if (abs (esqreal) .le. dsqerr) then
if ((itrun .eq. 0) .and. (esqreal .ne. 0.0)) itrun = 1
esqreal = 0.0
endif
esqimag = -ggimag
cbugc***DEBUG begins.
cbug 9705 format ('aptquar DEBUG. Dsq = ',1pe22.14,' + i *',1pe22.14)
cbug 9706 format ('aptquar DEBUG. Esq = ',1pe22.14,' + i *',1pe22.14)
cbug write ( 3, 9705) dsqreal, dsqimag
cbug write ( 3, 9706) esqreal, esqimag
cbugc***DEBUG ends.
if (dsqimag .eq. 0.0) then ! dsq is real.
if (dsqreal .gt. 0.0) then
dcapreal = sqrt (dsqreal)
dcapimag = 0.0
elseif (dsqreal .eq. 0.0) then
dcapreal = 0.0
dcapimag = 0.0
else
dcapreal = 0.0
dcapimag = sqrt (-dsqreal)
endif
else ! dsq is complex.
dmag = sqrt (dsqreal**2 + dsqimag**2)
if (abs (dsqreal) .gt. abs (dsqimag)) then
if (dsqreal .gt. 0.0) then
cosfun = 1.0 + dsqreal / dmag
sinth = dsqimag / dmag
cosdang = sqrt (0.5 * cosfun)
sindang = sinth / sqrt (2.0 * cosfun)
else
cosfun = 1.0 - dsqreal / dmag
sinth = dsqimag / dmag
cosdang = sinth / sqrt (2.0 * cosfun)
sindang = sqrt (0.5 * cosfun)
endif
else
dang = 0.5 * atan2 (dsqimag, dsqreal)
cosdang = cos (dang)
sindang = sin (dang)
endif
dmagrt = sqrt (dmag)
dcapreal = dmagrt * cosdang
dcapimag = dmagrt * sindang
endif
if (esqimag .eq. 0.0) then ! esq is real.
if (esqreal .gt. 0.0) then
ecapreal = sqrt (esqreal)
ecapimag = 0.0
elseif (esqreal .eq. 0.0) then
ecapreal = 0.0
ecapimag = 0.0
else
ecapreal = 0.0
ecapimag = sqrt (-esqreal)
endif
else ! esq is complex.
emag = sqrt (esqreal**2 + esqimag**2)
if (abs (esqreal) .gt. abs (esqimag)) then
if (esqreal .gt. 0.0) then
cosfun = 1.0 + esqreal / emag
sinth = esqimag / emag
coseang = sqrt (0.5 * cosfun)
sineang = sinth / sqrt (2.0 * cosfun)
else
cosfun = 1.0 - esqreal / emag
sinth = esqimag / emag
coseang = sinth / sqrt (2.0 * cosfun)
sineang = sqrt (0.5 * cosfun)
endif
else
eang = 0.5 * atan2 (esqimag, esqreal)
coseang = cos (eang)
sineang = sin (eang)
endif
emagrt = sqrt (emag)
ecapreal = emagrt * coseang
ecapimag = emagrt * sineang
endif
cbugc***DEBUG begins.
cbug 9707 format ('aptquar DEBUG. dcap= ',1pe22.14,' + i *',1pe22.14)
cbug 9708 format ('aptquar DEBUG. ecap= ',1pe22.14,' + i *',1pe22.14)
cbug write ( 3, 9707) dcapreal, dcapimag
cbug write ( 3, 9708) ecapreal, ecapimag
cbugc***DEBUG ends.
c.... Find the four roots of the quartic equation.
xroot(1) = -0.25 * c3 + 0.5 * (-rreal - ecapreal)
xrerr = tol * (0.25 * abs (c3) +
& 0.5 * (abs (rreal) + abs (ecapreal)))
if (abs (xroot(1)) .le. xrerr) then
if ((itrun .eq. 0) .and. (xroot(1) .ne. 0.0)) itrun = 1
xroot(1) = 0.0
endif
yroot(1) = 0.5 * (-rimag - ecapimag)
yrerr = tol * (0.5 * (abs (rimag) + abs (ecapimag)))
if (abs (yroot(1)) .le. yrerr) then
if ((itrun .eq. 0) .and. (yroot(1) .ne. 0.0)) itrun = 1
yroot(1) = 0.0
endif
xroot(2) = -0.25 * c3 + 0.5 * (-rreal + ecapreal)
if (abs (xroot(2)) .le. xrerr) then
if ((itrun .eq. 0) .and. (xroot(2) .ne. 0.0)) itrun = 1
xroot(2) = 0.0
endif
yroot(2) = 0.5 * (-rimag + ecapimag)
if (abs (yroot(2)) .le. yrerr) then
if ((itrun .eq. 0) .and. (yroot(2) .ne. 0.0)) itrun = 1
yroot(2) = 0.0
endif
xroot(3) = -0.25 * c3 + 0.5 * ( rreal - dcapreal)
xrerr = tol * (0.25 * abs (c3) +
& 0.5 * (abs (rreal) + abs (dcapreal)))
if (abs (xroot(3)) .le. xrerr) then
if ((itrun .eq. 0) .and. (xroot(3) .ne. 0.0)) itrun = 1
xroot(3) = 0.0
endif
yroot(3) = 0.5 * ( rimag - dcapimag)
yrerr = tol * (0.5 * (abs (rimag) + abs (dcapimag)))
if (abs (yroot(3)) .le. yrerr) then
if ((itrun .eq. 0) .and. (yroot(3) .ne. 0.0)) itrun = 1
yroot(3) = 0.0
endif
xroot(4) = -0.25 * c3 + 0.5 * ( rreal + dcapreal)
if (abs (xroot(4)) .le. xrerr) then
if ((itrun .eq. 0) .and. (xroot(4) .ne. 0.0)) itrun = 1
xroot(4) = 0.0
endif
yroot(4) = 0.5 * ( rimag + dcapimag)
if (abs (yroot(4)) .le. yrerr) then
if ((itrun .eq. 0) .and. (yroot(4) .ne. 0.0)) itrun = 1
yroot(4) = 0.0
endif
cbugc***DEBUG begins.
cbug 9810 format (' Root ',i1,' r,i=',1p2e22.14)
cbug write ( 3, 9810) (n, xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.
c.... Find the root of maximum magnitude.
nmax = 0
rsqmax = 0.0
do 180 n = 1, 4
rsqn = xroot(n)**2 + yroot(n)**2
if (rsqn .gt. rsqmax) then
nmax = n
rsqmax = rsqn
endif
180 continue
cbugc***DEBUG begins.
cbug 9812 format (' Root ',i1,' (max) =',1p2e22.14,' (real, imag)' )
cbug
cbug write ( 3, 9812) nmax, xroot(nmax), yroot(nmax)
cbugc***DEBUG ends.
c.... Factor the maximum root(s) out of the quartic equation.
if (yroot(nmax) .eq. 0.0) then ! Maximum root is real.
c.... Factor out one real root.
x1 = xroot(nmax)
x1a = abs (x1)
c0a = abs (c0)
c1a = abs (c1)
c2a = abs (c2)
c3a = abs (c3)
d0err = tol * (c1a + x1a * (c2a + x1a * (c3a + x1a)))
d0err2 = tol * (c0a / x1a)
if (d0err2 .le. d0err) then
d0 = -c0 / x1
else
d0 = c1 + x1 * (c2 + x1 * (c3 + x1))
endif
d1err = tol * (c2a + x1a * (c3a + x1a))
d1err2 = tol * (c0a + c1a * x1a) / x1**2
if (d1err2 .le. d1err) then
d1 = -(c0 + c1 * x1) / x1**2
else
d1 = c2 + x1 * (c3 + x1)
endif
d2err = tol * abs (c3a + x1a)
d2err2 = tol * (c0a + x1a * (c1 * x1a + c2 * x1) / x1a**3)
if (d2err2 .le. d2err) then
d2 = -(c0 + x1 * (c1 + c2 * x1)) / x1**3
else
d2 = c3 + x1
endif
cbugc***DEBUG begins.
cbug 9808 format (' d0,d1,d2= ',1p3e22.14)
cbug write ( 3, 9808) d0, d1, d2
cbugc***DEBUG ends.
xroot(1) = x1
yroot(1) = 0.0
call aptcubs (d0, d1, d2, tol, nroots, xroot(2), yroot(2),
& atype(2), itrun)
if (nroots .eq. 0) then
kbad = 1
itrun = 2
go to 395
elseif (atype(2) .eq. 'triple') then
xroot(3) = xroot(2)
yroot(3) = yroot(2)
xroot(4) = xroot(2)
yroot(4) = yroot(2)
elseif ((atype(2) .eq. 'dbl min') .or.
& (atype(2) .eq. 'dbl max')) then
xroot(4) = xroot(2)
yroot(4) = yroot(2)
elseif ((atype(3) .eq. 'dbl min') .or.
& (atype(3) .eq. 'dbl max')) then
xroot(4) = xroot(3)
yroot(4) = yroot(3)
endif
cbugc***DEBUG begins.
cbug 9809 format ('aptquar DEBUG. From aptcubs: nroots=',i2)
cbug write ( 3, 9809) nroots
cbug 9811 format (' Root ',i1,' (',a7,')',1p2e22.14,' (real, imag)' )
cbug write ( 3, 9811) (n, atype(n), xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.
cbugc***DEBUG begins.
cbug write ( 3, 9810) (n, xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.
else ! Maximum roots are complex conjugates.
c.... Factor out two conjugate complex roots.
x1 = xroot(nmax)
y1 = yroot(nmax)
d0 = c0 / (x1**2 + y1**2)
d1 = c3 + 2.0 * x1
d1err = tol * (abs (c3) + 2.0 * abs (x1))
if (abs (d1) .le. d1err) d1 = 0.0
d2 = 1.0
cbugc***DEBUG begins.
cbug 9910 format (' d0,d1 = ',1p2e22.14)
cbug write ( 3, 9910) d0, d1
cbugc***DEBUG ends.
call aptqrts (0, d2, d1, d0, qq, tol,
& nrootx, xnew1, xnew2, itrun)
xroot(1) = x1
yroot(1) = -abs (y1)
xroot(2) = x1
yroot(2) = abs (y1)
if (nrootx .eq. 1) then
xroot(3) = xnew1
yroot(3) = 0.0
xroot(4) = xnew1
yroot(4) = 0.0
elseif (nrootx .eq. 2) then
xroot(3) = xnew1
yroot(3) = 0.0
xroot(4) = xnew2
yroot(4) = 0.0
else
xroot(3) = xnew1
yroot(3) = -abs (xnew2)
xroot(4) = xnew1
yroot(4) = abs (xnew2)
endif
cbugc***DEBUG begins.
cbug write ( 3, 9810) (n, xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.
endif ! Factored out largest real root.
c.... Reorder the roots in ascending order of absolute magnitude.
nzero = 0
do 192 n = 1, 4
xsqn = xroot(n)**2 + yroot(n)**2
if (xsqn .eq. 0.0) nzero = nzero + 1
192 continue
do 195 n = 1, 3
do 190 nn = n + 1, 4
xsqn = xroot(n)**2 + yroot(n)**2
xsqnn = xroot(nn)**2 + yroot(nn)**2
if (xsqnn .lt. xsqn) then
xrootx = xroot(n)
yrootx = yroot(n)
xroot(n) = xroot(nn)
yroot(n) = yroot(nn)
xroot(nn) = xrootx
yroot(nn) = yrootx
endif
190 continue
195 continue
cbugc***DEBUG begins.
cbug 9815 format ('aptquar DEBUG. nzero=',i1)
cbug write ( 3, 9815) nzero
cbug write ( 3, 9810) (n, xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.
c.... See if there is a double root at zero. If so, recalculate.
if (nzero .eq. 2) then
d0 = c0 / (xroot(3) * xroot(4) - yroot(3) * yroot(4))
d1 = c3 + xroot(3) + xroot(4)
d1err = tol * (abs (c3) + abs (xroot(3)) + abs (xroot(4)))
if (abs (d1) .le. d1err) then
if ((d1 .ne. 0.0) .and. (itrun .eq. 0)) itrun = 1
d1 = 0.0
endif
d2 = 1.0
call aptqrts (0, d2, d1, d0, qq, tol,
& nrootx, xnew1, xnew2, itrun)
if (nrootx .eq. 2) then
xroot(1) = xnew1
yroot(1) = 0.0
xroot(2) = xnew2
yroot(2) = 0.0
else
xroot(1) = xnew1
yroot(1) = -xnew2
xroot(2) = xnew1
yroot(2) = xnew2
endif
cbugc***DEBUG begins.
cbug 9818 format ('aptquar DEBUG. Replaced double root at zero.')
cbug write ( 3, 9818)
cbug write ( 3, 9810) (n, xroot(n), yroot(n), n = 1, 4)
cbugc***DEBUG ends.
endif
c.... Reorder the roots in ascending order of real, then complex. Find types.
200 do 220 n = 1, 3 ! Loop over roots
do 210 nn = n + 1, 4 ! Loop over roots.
c.... Put real roots in ascending order.
if ((yroot(n) .eq. 0.0) .and.
& (yroot(nn) .eq. 0.0)) then
atype(n) = 'single'
atype(nn) = 'single'
if (xroot(n) .gt. xroot(nn)) then
xrootx = xroot(n)
xroot(n) = xroot(nn)
xroot(nn) = xrootx
endif
c.... Put complex roots in ascending order of real, then imaginary parts.
elseif ((yroot(n) .ne. 0.0) .and.
& (yroot(nn) .ne. 0.0)) then
atype(n) = 'complex'
atype(nn) = 'complex'
if (xroot(n) .gt. xroot(nn)) then
xrootx = xroot(n)
xroot(n) = xroot(nn)
xroot(nn) = xrootx
yrootx = yroot(n)
yroot(n) = yroot(nn)
yroot(nn) = yrootx
elseif (xroot(n) .eq. xroot(nn)) then
if (abs (yroot(n)) .gt. abs (yroot(nn))) then
xrootx = xroot(n)
xroot(n) = xroot(nn)
xroot(nn) = xrootx
yrootx = yroot(n)
yroot(n) = yroot(nn)
yroot(nn) = yrootx
elseif ((abs (yroot(n)) .eq. abs (yroot(nn))) .and.
& (yroot(nn) .lt. 0.0)) then
xrootx = xroot(n)
xroot(n) = xroot(nn)
xroot(nn) = xrootx
yrootx = yroot(n)
yroot(n) = yroot(nn)
yroot(nn) = yrootx
endif
endif
c.... Put real roots before complex roots.
elseif ((yroot(n) .ne. 0.0) .and.
& (yroot(nn) .eq. 0.0)) then
atype(n) = 'single'
atype(nn) = 'complex'
xrootx = xroot(n)
xroot(n) = xroot(nn)
xroot(nn) = xrootx
yrootx = yroot(n)
yroot(n) = yroot(nn)
yroot(nn) = yrootx
elseif ((yroot(n) .eq. 0.0) .and.
& (yroot(nn) .ne. 0.0)) then
atype(n) = 'single'
atype(nn) = 'complex'
endif
210 continue ! Loop over roots.
220 continue ! Loop over roots
c.... Reconstruct the coefficients.
kbad = 0
c0p = xroot(1) * xroot(2) * xroot(3) * xroot(4) -
& (xroot(1) * xroot(2) * yroot(3) * yroot(4) +
& xroot(2) * xroot(3) * yroot(4) * yroot(1) +
& xroot(3) * xroot(4) * yroot(1) * yroot(2) +
& xroot(4) * xroot(1) * yroot(2) * yroot(3) +
& xroot(1) * xroot(3) * yroot(2) * yroot(4) +
& xroot(2) * xroot(4) * yroot(1) * yroot(3)) +
& yroot(1) * yroot(2) * yroot(3) * yroot(4)
c1p = -(xroot(1) * xroot(2) * (xroot(3) + xroot(4)) +
& xroot(3) * xroot(4) * (xroot(1) + xroot(2))) +
& xroot(1) * (yroot(2) * (yroot(3) + yroot(4)) +
& yroot(3) * yroot(4) ) +
& xroot(2) * (yroot(3) * (yroot(4) + yroot(1)) +
& yroot(4) * yroot(1) ) +
& xroot(3) * (yroot(4) * (yroot(1) + yroot(2)) +
& yroot(1) * yroot(2) ) +
& xroot(4) * (yroot(1) * (yroot(2) + yroot(3)) +
& yroot(2) * yroot(3) )
c2p = xroot(1) * xroot(2) + xroot(2) * xroot(3) +
& xroot(3) * xroot(4) + xroot(4) * xroot(1) +
& xroot(1) * xroot(3) + xroot(2) * xroot(4) -
& (yroot(1) * yroot(2) + yroot(2) * yroot(3) +
& yroot(3) * yroot(4) + yroot(4) * yroot(1) +
& yroot(1) * yroot(3) + yroot(2) * yroot(4) )
c3p = -(xroot(1) + xroot(2) + xroot(3) + xroot(4))
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
if (abs (c3 - c3p) .gt. tol * (abs (c3) + abs (c3p))) kbad = 1
cbugc***DEBUG begins.
cbug 9930 format ('aptquar DEBUG. itrun =',i2, /
cbug & ' c0,c1(input) =',1p2e22.14 /
cbug & ' c0,c1(recalc)=',1p2e22.14 //
cbug & ' c2,c3(input) =',1p2e22.14 /
cbug & ' c2,c3(recalc)=',1p2e22.14 )
cbug write ( 3, 9930) itrun, c0, c1, c0p, c1p, c2, c3, c2p, c3p
cbugC***DEBUG ends.
c.... Find and classify unique roots.
nroots = 1
do 230 n = 2, 4
xdiff = abs (xroot(n) - xroot(n-1))
ydiff = abs (yroot(n) - yroot(n-1))
xderr = tol * (abs (xroot(n)) + abs (xroot(n-1)))
yderr = tol * (abs (yroot(n)) + abs (yroot(n-1)))
if ((xdiff .le. xderr) .and.
& (ydiff .le. yderr)) then
if (atype(nroots) .eq. 'single') then
curve = c2 + 3.0 * c3 * xroot(n) + 6.0 * xroot(n)**2
if (curve .lt. 0.0) then
atype(nroots) = 'dbl max'
else
atype(nroots) = 'dbl min'
endif
elseif ((atype(nroots) .eq. 'dbl min') .or.
& (atype(nroots) .eq. 'dbl max')) then
atype(nroots) = 'triple'
elseif (atype(nroots) .eq. 'triple') then
atype(nroots) = 'quadrpl'
elseif (atype(nroots) .eq. 'complex') then
atype(nroots) = 'dbl com'
endif
else ! This root is unique.
nroots = nroots + 1
atype(nroots) = atype(n)
xroot(nroots) = xroot(n)
yroot(nroots) = yroot(n)
endif
230 continue
c.... Count real and imaginary roots.
do 250 n = 1, nroots
if (yroot(n) .eq. 0.0) then
nroot = nroot + 1
else
nrooti = nrooti + 1
endif
250 continue
if (nroots .lt. 4) then
do 260 n = 4, nroots + 1, -1
atype(n) = 'unknown'
260 continue
endif
c.... Store any duplicate roots in root locations nroots + 1, 4.
next = nroots
do 270 n = 1, nroots
if ((atype(n) .eq. 'dbl min') .or.
& (atype(n) .eq. 'dbl max') .or.
& (atype(n) .eq. 'dbl com')) then
next = next + 1
xroot(next) = xroot(n)
yroot(next) = yroot(n)
elseif (atype(n) .eq. 'triple') then
next = next + 1
xroot(next) = xroot(n)
yroot(next) = yroot(n)
next = next + 1
xroot(next) = xroot(n)
yroot(next) = yroot(n)
elseif (atype(n) .eq. 'quadrpl') then
next = next + 1
xroot(next) = xroot(n)
yroot(next) = yroot(n)
next = next + 1
xroot(next) = xroot(n)
yroot(next) = yroot(n)
next = next + 1
xroot(next) = xroot(n)
yroot(next) = yroot(n)
endif
270 continue
cbugc***DEBUG begins.
cbug 9920 format (/ 'aptquar results. nroot=',i2,' nrooti=',i3)
cbug write ( 3, 9920) nroot, nrooti
cbugc***DEBUG ends.
do 390 n = 1, nroots
frres = 0.0
frres = xroot(n)**4 + c3 * xroot(n)**3 + c2 * xroot(n)**2 +
& c1 * xroot(n) + c0 +
& yroot(n)**2 * (yroot(n)**2 - 6.0 * xroot(n)**2 -
& 3.0 * c3 * xroot(n) - c2)
frerr = tol * (xroot(n)**4 + abs (c3 * xroot(n)**3) +
& abs (c2) * xroot(n)**2 + abs (c1 * xroot(n)) +
& abs (c0) +
& yroot(n)**2 * (yroot(n)**2 + 6.0 * xroot(n)**2 +
& 3.0 * abs (c3 * xroot(n)) + abs (c2)))
fires = yroot(n) * (4.0 * xroot(n) *
& (xroot(n)**2 - yroot(n)**2) +
& c3 * (3.0 * xroot(n)**2 - yroot(n)**2) +
& 2.0 * c2 * xroot(n) + c1)
fierr = tol * (abs (yroot(n)) * (4.0 * abs (xroot(n)) *
& (xroot(n)**2 + yroot(n)**2) +
& abs (c3) * (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 if (kbad .eq. 1) then ! Try to improve with Newton's method.
cbug
cbugc.... Find the partial derivatives of the real and imaginary parts.
cbug
cbug xslop = c1 + 2.0 * c2 * xroot(n) + 3.0 * c3 * xroot(n)**2 +
cbug & 4.0 * xroot(n)**3 -
cbug & 3.0 * yroot(n)**2 * (c3 + 4.0 * xroot(n))
cbug
cbug xserr = tol * (abs (c1) + 2.0 * abs (c2 * xroot(n)) +
cbug & 3.0 * abs (c3) * xroot(n)**2 +
cbug & 4.0 * xroot(n)**3 +
cbug & 3.0 * yroot(n)**2 * (abs (c3) + 4.0 * abs (xroot(n))))
cbug
cbug yslop = 4.0 * yroot(n)**3 - 2.0 * yroot(n) *
cbug & (c2 + 3.0 * c3 * xroot(n) + 6.0 * xroot(n)**2)
cbug
cbug yserr = tol * (4.0 * abs (yroot(n)**3) +
cbug & 2.0 * abs (yroot(n)) *
cbug & (abs (c2) + 3.0 * abs (c3 * xroot(n)) +
cbug & 6.0 * xroot(n)**2))
cbug
cbug if (abs (xslop) .le. xserr) xslop = 0.0
cbug if (abs (yslop) .le. yserr) yslop = 0.0
cbug
cbug slsq = xslop**2 + yslop**2
cbug xdel = 0.0
cbug ydel = 0.0
cbug if (slsq .gt. 0.0) then
cbug xdel = (fires * yslop - frres * xslop) / slsq
cbug ydel = (frres * yslop + fires * xslop) / slsq
cbug endif
cbug xcorr = xroot(n) + xdel
cbug ycorr = yroot(n) + ydel
cbug endif ! Tried Newton's method.
cbugc***DEBUG ends.
cbugc***DEBUG begins.
cbug 9921 format (' Root ',i1,' r,i=',1p2e22.14, 2x,a8 /
cbug & ' P(z) r,i=',1p2e22.14 /
cbug & ' allowed r,i=',1p2e22.14 )
cbug 9922 format (' Slope r,i=',1p2e22.14 /
cbug & ' Delta r,i=',1p2e22.14 /
cbug & ' Iterate r,i=',1p2e22.14 )
cbug write ( 3, 9921) n, xroot(n), yroot(n), atype(n), frres, fires,
cbug & frerr, fierr
cbug if ((kbad .eq. 1) .and.
cbug & ((xdel .ne. 0.0) .or.
cbug & (ydel .ne. 0.0))) then
cbug write ( 3, 9922) xslop, yslop, xdel, ydel, xcorr, ycorr
cbug endif
cbugc***DEBUG ends.
390 continue
395 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 aptquar. (+1 line.)
end
UCRL-WEB-209832