mat_wsqrt Subroutine

subroutine mat_wsqrt(x_real,x_imag,y_real,y_imag)

! ident_21="@(#)M_LA::mat_wsqrt(3fp): y = sqrt(x) with y_real .ge. 0.0 and sign(y_imag) .eq. sign(x_imag)"

doubleprecision,intent(in)  :: x_real
doubleprecision,intent(in)  :: x_imag
doubleprecision,intent(out) :: y_real
doubleprecision,intent(out) :: y_imag
doubleprecision             :: s
doubleprecision             :: tr
doubleprecision             :: ti
!
   tr = x_real
   ti = x_imag
   s = dsqrt(0.5d0*(mat_pythag(tr,ti) + dabs(tr)))
   if (tr .ge. 0.0d0) y_real = mat_flop(s)
   if (ti .lt. 0.0d0) s = -s
   if (tr .le. 0.0d0) y_imag = mat_flop(s)
   if (tr .lt. 0.0d0) y_real = mat_flop(0.5d0*(ti/y_imag))
   if (tr .gt. 0.0d0) y_imag = mat_flop(0.5d0*(ti/y_real))
end subroutine mat_wsqrt