stdlib_drotg Subroutine

     pure subroutine stdlib_drotg( a, b, c, s )
     !! The computation uses the formulas
     !! sigma = sgn(a)    if |a| >  |b|
     !! = sgn(b)    if |b| >= |a|
     !! r = sigma*sqrt( a**2 + b**2 )
     !! c = 1; s = 0      if r = 0
     !! c = a/r; s = b/r  if r != 0
     !! The subroutine also computes
     !! z = s    if |a| > |b|,
     !! = 1/c  if |b| >= |a| and c != 0
     !! = 1    if c = 0
     !! This allows c and s to be reconstructed from z as follows:
     !! If z = 1, set c = 0, s = 1.
     !! If |z| < 1, set c = sqrt(1 - z**2) and s = z.
     !! If |z| > 1, set c = 1/z and s = sqrt( 1 - c**2).
        ! -- reference blas level1 routine --
        ! -- reference blas is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
        ! Constants 
        integer, parameter :: wp = kind(1._dp)
        ! Scaling Constants 
        ! Scalar Arguments 
        real(dp), intent(inout) :: a, b
        real(dp), intent(out) :: c, s
        ! Local Scalars 
        real(dp) :: anorm, bnorm, scl, sigma, r, z
        anorm = abs(a)
        bnorm = abs(b)
        if( bnorm == zero ) then
           c = one
           s = zero
           b = zero
        else if( anorm == zero ) then
           c = zero
           s = one
           a = b
           b = one
        else
           scl = min( safmax, max( safmin, anorm, bnorm ) )
           if( anorm > bnorm ) then
              sigma = sign(one,a)
           else
              sigma = sign(one,b)
           end if
           r = sigma*( scl*sqrt((a/scl)**2 + (b/scl)**2) )
           c = a/r
           s = b/r
           if( anorm > bnorm ) then
              z = s
           else if( c /= zero ) then
              z = one/c
           else
              z = one
           end if
           a = r
           b = z
        end if
        return
     end subroutine stdlib_drotg