double precision function dznrm2( n, zx, incx)
logical imag, scale
integer next
double precision cutlo, cuthi, hitest, sum, xmax, absx, zero, one
double complex zx(1)
double precision dreal,dimag
double complex zdumr,zdumi
dreal(zdumr) = zdumr
dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
data zero, one /0.0d0, 1.0d0/
c
c unitary norm of the complex n-vector stored in zx() with storage
c increment incx .
c if n .le. 0 return with result = 0.
c if n .ge. 1 then incx must be .ge. 1
c
c c.l.lawson , 1978 jan 08
c
c four phase method using two built-in constants that are
c hopefully applicable to all machines.
c cutlo = maximum of sqrt(u/eps) over all known machines.
c cuthi = minimum of sqrt(v) over all known machines.
c where
c eps = smallest no. such that eps + 1. .gt. 1.
c u = smallest positive no. (underflow limit)
c v = largest no. (overflow limit)
c
c brief outline of algorithm..
c
c phase 1 scans zero components.
c move to phase 2 when a component is nonzero and .le. cutlo
c move to phase 3 when a component is .gt. cutlo
c move to phase 4 when a component is .ge. cuthi/m
c where m = n for x() real and m = 2*n for complex.
c
c values for cutlo and cuthi..
c from the environmental parameters listed in the imsl converter
c document the limiting values are as follows..
c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are
c univac and dec at 2**(-103)
c thus cutlo = 2**(-51) = 4.44089e-16
c cuthi, s.p. v = 2**127 for univac, honeywell, and dec.
c thus cuthi = 2**(63.5) = 1.30438e19
c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec.
c thus cutlo = 2**(-33.5) = 8.23181d-11
c cuthi, d.p. same as s.p. cuthi = 1.30438d19
c data cutlo, cuthi / 8.232d-11, 1.304d19 /
c data cutlo, cuthi / 4.441e-16, 1.304e19 /
data cutlo, cuthi / 8.232d-11, 1.304d19 /
c
if(n .gt. 0) go to 10
dznrm2 = zero
go to 300
c
10 assign 30 to next
sum = zero
nn = n * incx
c begin main loop
do 210 i=1,nn,incx
absx = dabs(dreal(zx(i)))
imag = .false.
go to next,(30, 50, 70, 90, 110)
30 if( absx .gt. cutlo) go to 85
assign 50 to next
scale = .false.
c
c phase 1. sum is zero
c
50 if( absx .eq. zero) go to 200
if( absx .gt. cutlo) go to 85
c
c prepare for phase 2.
assign 70 to next
go to 105
c
c prepare for phase 4.
c
100 assign 110 to next
sum = (sum / absx) / absx
105 scale = .true.
xmax = absx
go to 115
c
c phase 2. sum is small.
c scale to avoid destructive underflow.
c
70 if( absx .gt. cutlo ) go to 75
c
c common code for phases 2 and 4.
c in phase 4 sum is large. scale to avoid overflow.
c
110 if( absx .le. xmax ) go to 115
sum = one + sum * (xmax / absx)**2
xmax = absx
go to 200
c
115 sum = sum + (absx/xmax)**2
go to 200
c
c
c prepare for phase 3.
c
75 sum = (sum * xmax) * xmax
c
85 assign 90 to next
scale = .false.
c
c for real or d.p. set hitest = cuthi/n
c for complex set hitest = cuthi/(2*n)
c
hitest = cuthi/float( n )
c
c phase 3. sum is mid-range. no scaling.
c
90 if(absx .ge. hitest) go to 100
sum = sum + absx**2
200 continue
c control selection of real and imaginary parts.
c
if(imag) go to 210
absx = dabs(dimag(zx(i)))
imag = .true.
go to next,( 50, 70, 90, 110 )
c
210 continue
c
c end of main loop.
c compute square root and adjust for scaling.
c
dznrm2 = dsqrt(sum)
if(scale) dznrm2 = dznrm2 * xmax
300 continue
return
end