92   integer, 
parameter :: wp = kind(1.e0)
 
   99   real(wp), 
parameter :: zero = 0.0_wp
 
  100   real(wp), 
parameter :: one  = 1.0_wp
 
  103   real(wp), 
parameter :: safmin = real(radix(0._wp),wp)**max( &
 
  104      minexponent(0._wp)-1, &
 
  105      1-maxexponent(0._wp) &
 
  107   real(wp), 
parameter :: safmax = real(radix(0._wp),wp)**max( &
 
  108      1-minexponent(0._wp), &
 
  109      maxexponent(0._wp)-1 &
 
  113   real(wp) :: a, b, c, s
 
  116   real(wp) :: anorm, bnorm, scl, sigma, r, z
 
  120   if( bnorm == zero ) 
then 
  124   else if( anorm == zero ) 
then 
  130      scl = min( safmax, max( safmin, anorm, bnorm ) )
 
  131      if( anorm > bnorm ) 
then 
  136      r = sigma*( scl*sqrt((a/scl)**2 + (b/scl)**2) )
 
  139      if( anorm > bnorm ) 
then 
  141      else if( c /= zero ) 
then