Computes the roots of a cubic polynomial of the form
a(1) + a(2)*z + a(3)*z**2 + a(4)*z**3 = 0
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in), | dimension(4) | :: | a |
coefficients |
|
| real(kind=wp), | intent(out), | dimension(3) | :: | zr |
real components of roots |
|
| real(kind=wp), | intent(out), | dimension(3) | :: | zi |
imaginary components of roots |
subroutine dcbcrt(a, zr, zi) implicit none real(wp), dimension(4), intent(in) :: a !! coefficients real(wp), dimension(3), intent(out) :: zr !! real components of roots real(wp), dimension(3), intent(out) :: zi !! imaginary components of roots real(wp) :: arg, c, cf, d, p, p1, q, q1, & r, ra, rb, rq, rt, r1, s, sf, sq, sum, & t, tol, t1, w, w1, w2, x, x1, x2, x3, y, & y1, y2, y3 real(wp), parameter :: sqrt3 = sqrt(3.0_wp) if (a(4) == 0.0_wp) then ! quadratic equation call dqdcrt(a(1:3), zr(1:2), zi(1:2)) !there are only two roots, so just duplicate the second one: zr(3) = zr(2) zi(3) = zi(2) else if (a(1) == 0.0_wp) then ! quadratic zr(1) = 0.0_wp zi(1) = 0.0_wp call dqdcrt(a(2:4), zr(2:3), zi(2:3)) else p = a(3)/(3.0_wp*a(4)) q = a(2)/a(4) r = a(1)/a(4) tol = 4.0_wp*eps c = 0.0_wp t = a(2) - p*a(3) if (abs(t) > tol*abs(a(2))) c = t/a(4) t = 2.0_wp*p*p - q if (abs(t) <= tol*abs(q)) t = 0.0_wp d = r + p*t if (abs(d) <= tol*abs(r)) then !case when d = 0 zr(1) = -p zi(1) = 0.0_wp w = sqrt(abs(c)) if (c < 0.0_wp) then if (p /= 0.0_wp) then x = -(p + sign(w, p)) zr(3) = x zi(2) = 0.0_wp zi(3) = 0.0_wp t = 3.0_wp*a(1)/(a(3)*x) if (abs(p) > abs(t)) then zr(2) = zr(1) zr(1) = t else zr(2) = t end if else zr(2) = w zr(3) = -w zi(2) = 0.0_wp zi(3) = 0.0_wp end if else zr(2) = -p zr(3) = zr(2) zi(2) = w zi(3) = -w end if else !set sq = (a(4)/s)**2 * (c**3/27 + d**2/4) s = max(abs(a(1)), abs(a(2)), abs(a(3))) p1 = a(3)/(3.0_wp*s) q1 = a(2)/s r1 = a(1)/s t1 = q - 2.25_wp*p*p if (abs(t1) <= tol*abs(q)) t1 = 0.0_wp w = 0.25_wp*r1*r1 w1 = 0.5_wp*p1*r1*t w2 = q1*q1*t1/27.0_wp if (w1 < 0.0_wp) then if (w2 < 0.0_wp) then sq = w + (w1 + w2) else w = w + w2 sq = w + w1 end if else w = w + w1 sq = w + w2 end if if (abs(sq) <= tol*w) sq = 0.0_wp rq = abs(s/a(4))*sqrt(abs(sq)) if (sq < 0.0_wp) then !all roots are real arg = atan2(rq, -0.5_wp*d) cf = cos(arg/3.0_wp) sf = sin(arg/3.0_wp) rt = sqrt(-c/3.0_wp) y1 = 2.0_wp*rt*cf y2 = -rt*(cf + sqrt3*sf) y3 = -(d/y1)/y2 x1 = y1 - p x2 = y2 - p x3 = y3 - p if (abs(x1) > abs(x2)) then t = x1 x1 = x2 x2 = t end if if (abs(x2) > abs(x3)) then t = x2 x2 = x3 x3 = t if (abs(x1) > abs(x2)) then t = x1 x1 = x2 x2 = t end if end if w = x3 if (abs(x2) < 0.1_wp*abs(x3)) then call roundoff() else if (abs(x1) < 0.1_wp*abs(x2)) x1 = -(r/x3)/x2 zr(1) = x1 zr(2) = x2 zr(3) = x3 zi(1) = 0.0_wp zi(2) = 0.0_wp zi(3) = 0.0_wp end if else !real and complex roots ra = dcbrt(-0.5_wp*d - sign(rq, d)) rb = -c/(3.0_wp*ra) t = ra + rb w = -p x = -p if (abs(t) > tol*abs(ra)) then w = t - p x = -0.5_wp*t - p if (abs(x) <= tol*abs(p)) x = 0.0_wp end if t = abs(ra - rb) y = 0.5_wp*sqrt3*t if (t > tol*abs(ra)) then if (abs(x) < abs(y)) then s = abs(y) t = x/y else s = abs(x) t = y/x end if if (s < 0.1_wp*abs(w)) then call roundoff() else w1 = w/s sum = 1.0_wp + t*t if (w1*w1 < 1.0e-2_wp*sum) w = -((r/sum)/s)/s zr(1) = w zr(2) = x zr(3) = x zi(1) = 0.0_wp zi(2) = y zi(3) = -y end if else !at least two roots are equal zi(1) = 0.0_wp zi(2) = 0.0_wp zi(3) = 0.0_wp if (abs(x) < abs(w)) then if (abs(x) < 0.1_wp*abs(w)) then call roundoff() else zr(1) = x zr(2) = x zr(3) = w end if else if (abs(w) < 0.1_wp*abs(x)) w = -(r/x)/x zr(1) = w zr(2) = x zr(3) = x end if end if end if end if end if end if contains !************************************************************* subroutine roundoff() !! here `w` is much larger in magnitude than the other roots. !! as a result, the other roots may be exceedingly inaccurate !! because of roundoff error. to deal with this, a quadratic !! is formed whose roots are the same as the smaller roots of !! the cubic. this quadratic is then solved. !! !! this code was written by william l. davis (nswc). implicit none real(wp), dimension(3) :: aq aq(1) = a(1) aq(2) = a(2) + a(1)/w aq(3) = -a(4)*w call dqdcrt(aq, zr(1:2), zi(1:2)) zr(3) = w zi(3) = 0.0_wp if (zi(1) /= 0.0_wp) then zr(3) = zr(2) zi(3) = zi(2) zr(2) = zr(1) zi(2) = zi(1) zr(1) = w zi(1) = 0.0_wp end if end subroutine roundoff !************************************************************* end subroutine dcbcrt