pythag Function

private pure function pythag(a, b)

Compute the complex square root of a complex number without destructive overflow or underflow.

Finds sqrt(A**2+B**2) without overflow or destructive underflow

REVISION HISTORY (YYMMDD)

  • 811101 DATE WRITTEN
  • 890531 Changed all specific intrinsics to generic. (WRB)
  • 891214 Prologue converted to Version 4.0 format. (BAB)
  • 900402 Added TYPE section. (WRB)
  • Jacob Williams, 9/14/2022 : modernized this code

Arguments

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

Return Value real(kind=wp)


Called by

proc~~pythag~~CalledByGraph proc~pythag polyroots_module::pythag proc~comqr polyroots_module::comqr proc~comqr->proc~pythag proc~csroot polyroots_module::csroot proc~comqr->proc~csroot proc~csroot->proc~pythag proc~cpqr79 polyroots_module::cpqr79 proc~cpqr79->proc~comqr

Source Code

pure real(wp) function pythag(a, b)
    implicit none

    real(wp), intent(in) :: a, b

    real(wp) :: p, q, r, s, t

    p = max(abs(a), abs(b))
    q = min(abs(a), abs(b))

    if (q /= 0.0_wp) then
        do
            r = (q/p)**2
            t = 4.0_wp + r
            if (t == 4.0_wp) exit
            s = r/t
            p = p + 2.0_wp*p*s
            q = q*s
        end do
    end if
    pythag = p

end function pythag