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