90   integer, parameter :: wp = kind(1.e0)
   91   real(wp) :: SCNRM2
   92
   93
   94
   95
   96
   97
   98
   99   real(wp), parameter :: zero = 0.0_wp
  100   real(wp), parameter :: one  = 1.0_wp
  101   real(wp), parameter :: maxN = huge(0.0_wp)
  102
  103
  104   real(wp), parameter :: tsml = real(radix(0._wp), wp)**ceiling( &
  105       (minexponent(0._wp) - 1) * 0.5_wp)
  106   real(wp), parameter :: tbig = real(radix(0._wp), wp)**floor( &
  107       (maxexponent(0._wp) - digits(0._wp) + 1) * 0.5_wp)
  108   real(wp), parameter :: ssml = real(radix(0._wp), wp)**( - floor( &
  109       (minexponent(0._wp) - digits(0._wp)) * 0.5_wp))
  110   real(wp), parameter :: sbig = real(radix(0._wp), wp)**( - ceiling( &
  111       (maxexponent(0._wp) + digits(0._wp) - 1) * 0.5_wp))
  112
  113
  114   integer :: incx, n
  115
  116
  117   complex(wp) :: x(*)
  118
  119
  120   integer :: i, ix
  121   logical :: notbig
  122   real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
  123
  124
  125
  127   if( n <= 0 ) return
  128
  129   scl = one
  130   sumsq = zero
  131
  132
  133
  134
  135
  136
  137
  138
  139
  140   notbig = .true.
  141   asml = zero
  142   amed = zero
  143   abig = zero
  144   ix = 1
  145   if( incx < 0 ) ix = 1 - (n-1)*incx
  146   do i = 1, n
  147      ax = abs(real(x(ix)))
  148      if (ax > tbig) then
  149         abig = abig + (ax*sbig)**2
  150         notbig = .false.
  151      else if (ax < tsml) then
  152         if (notbig) asml = asml + (ax*ssml)**2
  153      else
  154         amed = amed + ax**2
  155      end if
  156      ax = abs(aimag(x(ix)))
  157      if (ax > tbig) then
  158         abig = abig + (ax*sbig)**2
  159         notbig = .false.
  160      else if (ax < tsml) then
  161         if (notbig) asml = asml + (ax*ssml)**2
  162      else
  163         amed = amed + ax**2
  164      end if
  165      ix = ix + incx
  166   end do
  167
  168
  169
  170
  171   if (abig > zero) then
  172
  173
  174
  175      if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
  176         abig = abig + (amed*sbig)*sbig
  177      end if
  178      scl = one / sbig
  179      sumsq = abig
  180   else if (asml > zero) then
  181
  182
  183
  184      if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
  185         amed = sqrt(amed)
  186         asml = sqrt(asml) / ssml
  187         if (asml > amed) then
  188            ymin = amed
  189            ymax = asml
  190         else
  191            ymin = asml
  192            ymax = amed
  193         end if
  194         scl = one
  195         sumsq = ymax**2*( one + (ymin/ymax)**2 )
  196      else
  197         scl = one / ssml
  198         sumsq = asml
  199      end if
  200   else
  201
  202
  203
  204      scl = one
  205      sumsq = amed
  206   end if
  207   scnrm2 = scl*sqrt( sumsq )
 
  208   return
real(wp) function scnrm2(n, x, incx)
SCNRM2