dqdcrt Subroutine

public pure subroutine dqdcrt(a, zr, zi)

Computes the roots of a quadratic polynomial of the form a(1) + a(2)*z + a(3)*z**2 = 0

See also

  • A. H. Morris, "NSWC Library of Mathematical Subroutines", Naval Surface Warfare Center, NSWCDD/TR-92/425, January 1993

See also

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in), dimension(3) :: a

coefficients
real(kind=wp), intent(out), dimension(2) :: zr

real components of roots
real(kind=wp), intent(out), dimension(2) :: zi

imaginary components of roots

Called by

proc~~dqdcrt~~CalledByGraph proc~dqdcrt polyroots_module::dqdcrt proc~dcbcrt polyroots_module::dcbcrt proc~dcbcrt->proc~dqdcrt proc~dqtcrt polyroots_module::dqtcrt proc~dqtcrt->proc~dcbcrt

Source Code

pure subroutine dqdcrt(a, zr, zi)

    implicit none

    real(wp), dimension(3), intent(in) :: a     !! coefficients
    real(wp), dimension(2), intent(out) :: zr   !! real components of roots
    real(wp), dimension(2), intent(out) :: zi   !! imaginary components of roots

    real(wp) :: d, r, w

    if (a(3) == 0.0_wp) then

        !it is really a linear equation:

        if (a(2) == 0.0_wp) then  !degenerate case, just return zeros

            zr = 0.0_wp
            zi = 0.0_wp

        else

            !there is only one root (real), so just duplicate it:

            zr(1) = -a(1)/a(2)
            zr(2) = zr(1)

            zi = 0.0_wp

        end if

    else

        if (a(1) == 0.0_wp) then

            zr(1) = 0.0_wp
            zi(1) = 0.0_wp
            zr(2) = -a(2)/a(3)
            zi(2) = 0.0_wp

        else

            d = a(2)*a(2) - 4.0_wp*a(1)*a(3)

            if (abs(d) <= 2.0_wp*eps*a(2)*a(2)) then

                !equal real roots

                zr(1) = -0.5_wp*a(2)/a(3)
                zr(2) = zr(1)
                zi(1) = 0.0_wp
                zi(2) = 0.0_wp

            else

                r = sqrt(abs(d))
                if (d < 0.0_wp) then

                    !complex roots

                    zr(1) = -0.5_wp*a(2)/a(3)
                    zr(2) = zr(1)
                    zi(1) = abs(0.5_wp*r/a(3))
                    zi(2) = -zi(1)

                else

                    !distinct real roots

                    zi(1) = 0.0_wp
                    zi(2) = 0.0_wp

                    if (a(2) /= 0.0_wp) then

                        w = -(a(2) + sign(r, a(2)))
                        zr(1) = 2.0_wp*a(1)/w
                        zr(2) = 0.5_wp*w/a(3)

                    else

                        zr(1) = abs(0.5_wp*r/a(3))
                        zr(2) = -zr(1)

                    end if

                end if

            end if

        end if

    end if

end subroutine dqdcrt