stdlib_crotg Subroutine

     pure subroutine stdlib_crotg( a, b, c, s )
     !! The computation uses the formulas
     !! |x| = sqrt( Re(x)**2 + Im(x)**2 )
     !! sgn(x) = x / |x|  if x /= 0
     !! = 1        if x  = 0
     !! c = |a| / sqrt(|a|**2 + |b|**2)
     !! s = sgn(a) * conjg(b) / sqrt(|a|**2 + |b|**2)
     !! When a and b are real and r /= 0, the formulas simplify to
     !! r = sgn(a)*sqrt(|a|**2 + |b|**2)
     !! c = a / r
     !! s = b / r
     !! the same as in SROTG when |a| > |b|.  When |b| >= |a|, the
     !! sign of c and s will be different from those computed by SROTG
     !! if the signs of a and b are not the same.
        ! -- 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..--
        ! Scaling Constants 
        ! Scalar Arguments 
        real(sp), intent(out) :: c
        complex(sp), intent(inout) :: a
        complex(sp), intent(in) :: b
        complex(sp), intent(out) :: s
        ! Local Scalars 
        real(sp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w
        complex(sp) :: f, fs, g, gs, r, t
        ! Intrinsic Functions 
        intrinsic :: abs,aimag,conjg,max,min,real,sqrt
        ! Statement Functions 
        real(sp) :: abssq
        ! Statement Function Definitions 
        abssq( t ) = real( t,KIND=sp)**2 + aimag( t )**2
        ! Executable Statements 
        f = a
        g = b
        if( g == czero ) then
           c = one
           s = czero
           r = f
        else if( f == czero ) then
           c = zero
           g1 = max( abs(real(g,KIND=sp)), abs(aimag(g)) )
           if( g1 > rtmin .and. g1 < rtmax ) then
              ! use unscaled algorithm
              g2 = abssq( g )
              d = sqrt( g2 )
              s = conjg( g ) / d
              r = d
           else
              ! use scaled algorithm
              u = min( safmax, max( safmin, g1 ) )
              uu = one / u
              gs = g*uu
              g2 = abssq( gs )
              d = sqrt( g2 )
              s = conjg( gs ) / d
              r = d*u
           end if
        else
           f1 = max( abs(real(f,KIND=sp)), abs(aimag(f)) )
           g1 = max( abs(real(g,KIND=sp)), abs(aimag(g)) )
     if( f1 > rtmin .and. f1 < rtmax .and.          g1 > rtmin .and. g1 < rtmax ) then
              ! use unscaled algorithm
              f2 = abssq( f )
              g2 = abssq( g )
              h2 = f2 + g2
              if( f2 > rtmin .and. h2 < rtmax ) then
                 d = sqrt( f2*h2 )
              else
                 d = sqrt( f2 )*sqrt( h2 )
              end if
              p = 1 / d
              c = f2*p
              s = conjg( g )*( f*p )
              r = f*( h2*p )
           else
              ! use scaled algorithm
              u = min( safmax, max( safmin, f1, g1 ) )
              uu = one / u
              gs = g*uu
              g2 = abssq( gs )
              if( f1*uu < rtmin ) then
                 ! f is not well-scaled when scaled by g1.
                 ! use a different scaling for f.
                 v = min( safmax, max( safmin, f1 ) )
                 vv = one / v
                 w = v * uu
                 fs = f*vv
                 f2 = abssq( fs )
                 h2 = f2*w**2 + g2
              else
                 ! otherwise use the same scaling for f and g.
                 w = one
                 fs = f*uu
                 f2 = abssq( fs )
                 h2 = f2 + g2
              end if
              if( f2 > rtmin .and. h2 < rtmax ) then
                 d = sqrt( f2*h2 )
              else
                 d = sqrt( f2 )*sqrt( h2 )
              end if
              p = 1 / d
              c = ( f2*p )*w
              s = conjg( gs )*( fs*p )
              r = ( fs*( h2*p ) )*u
           end if
        end if
        a = r
        return
     end subroutine stdlib_crotg