dcsqrt Subroutine

private subroutine dcsqrt(z, w)

w = sqrt(z) for the double precision complex number z

z and w are interpreted as double precision complex numbers. it is assumed that z(1) and z(2) are the real and imaginary parts of the complex number z, and that w(1) and w(2) are the real and imaginary parts of w.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: z(2)
real(kind=wp), intent(out) :: w(2)

Calls

proc~~dcsqrt~~CallsGraph proc~dcsqrt polyroots_module::dcsqrt proc~dcpabs polyroots_module::dcpabs proc~dcsqrt->proc~dcpabs

Called by

proc~~dcsqrt~~CalledByGraph proc~dcsqrt polyroots_module::dcsqrt proc~dqtcrt polyroots_module::dqtcrt proc~dqtcrt->proc~dcsqrt

Source Code

subroutine dcsqrt(z,w)

    real(wp),intent(in) :: z(2)
    real(wp),intent(out) :: w(2)

    real(wp) :: x , y , r

    x = z(1)
    y = z(2)
    if ( x<0 ) then
       if ( y/=0.0_wp ) then
          r = dcpabs(x,y)
          w(2) = sqrt(0.5_wp*(r-x))
          w(2) = sign(w(2),y)
          w(1) = 0.5_wp*y/w(2)
       else
          w(1) = 0.0_wp
          w(2) = sqrt(abs(x))
       endif
    elseif ( x==0.0_wp ) then
       if ( y/=0.0_wp ) then
          w(1) = sqrt(0.5_wp*abs(y))
          w(2) = sign(w(1),y)
       else
          w(1) = 0.0_wp
          w(2) = 0.0_wp
       endif
    elseif ( y/=0.0_wp ) then
       r = dcpabs(x,y)
       w(1) = sqrt(0.5_wp*(r+x))
       w(2) = 0.5_wp*y/w(1)
    else
       w(1) = sqrt(x)
       w(2) = 0.0_wp
    endif

 end subroutine dcsqrt