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