trunk/blas/dnrm2.f

C  dnrm2.f
C  is freely available from netlib. It is not subject to any GNU License
C  set by the authors of the ASCEND math programming system.
C  $Date: 1996/04/30 18:14:07 $ $Revision: 1.1.1.1 $
C
      double precision function dnrm2 ( n, dx, incx)
      integer          next
      double precision   dx(*), cutlo, cuthi, hitest, sum, xmax,zero,one
      data   zero, one /0.0d0, 1.0d0/
c
c     euclidean norm of the n-vector stored in dx() 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  dsqrt(u/eps)  over all known machines.
c         cuthi = minimum of  dsqrt(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
         dnrm2  = zero
         go to 300
c
   10 assign 30 to next
      sum = zero
      nn = n * incx
c                                                 begin main loop
      i = 1
   20    go to next,(30, 50, 70, 110)
   30 if( dabs(dx(i)) .gt. cutlo) go to 85
      assign 50 to next
      xmax = zero
c
c                        phase 1.  sum is zero
c
   50 if( dx(i) .eq. zero) go to 200
      if( dabs(dx(i)) .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 i = j
      assign 110 to next
      sum = (sum / dx(i)) / dx(i)
  105 xmax = dabs(dx(i))
      go to 115
c
c                   phase 2.  sum is small.
c                             scale to avoid destructive underflow.
c
   70 if( dabs(dx(i)) .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( dabs(dx(i)) .le. xmax ) go to 115
         sum = one + sum * (xmax / dx(i))**2
         xmax = dabs(dx(i))
         go to 200
c
  115 sum = sum + (dx(i)/xmax)**2
      go to 200
c
c
c                  prepare for phase 3.
c
   75 sum = (sum * xmax) * xmax
c
c
c     for real or d.p. set hitest = cuthi/n
c     for complex      set hitest = cuthi/(2*n)
c
   85 hitest = cuthi/float( n )
c
c                   phase 3.  sum is mid-range.  no scaling.
c
      do 95 j =i,nn,incx
      if(dabs(dx(j)) .ge. hitest) go to 100
   95    sum = sum + dx(j)**2
      dnrm2 = dsqrt( sum )
      go to 300
c
  200 continue
      i = i + incx
      if ( i .le. nn ) go to 20
c
c              end of main loop.
c
c              compute square root and adjust for scaling.
c
      dnrm2 = xmax * dsqrt(sum)
  300 continue
      return
      end

1	aw0a	1	C dnrm2.f
2			C is freely available from netlib. It is not subject to any GNU License
3			C set by the authors of the ASCEND math programming system.
4			C $Date: 1996/04/30 18:14:07 $ $Revision: 1.1.1.1 $
5			C
6			double precision function dnrm2 ( n, dx, incx)
7			integer next
8			double precision dx(*), cutlo, cuthi, hitest, sum, xmax,zero,one
9			data zero, one /0.0d0, 1.0d0/
10			c
11			c euclidean norm of the n-vector stored in dx() with storage
12			c increment incx .
13			c if n .le. 0 return with result = 0.
14			c if n .ge. 1 then incx must be .ge. 1
15			c
16			c c.l.lawson, 1978 jan 08
17			c
18			c four phase method using two built-in constants that are
19			c hopefully applicable to all machines.
20			c cutlo = maximum of dsqrt(u/eps) over all known machines.
21			c cuthi = minimum of dsqrt(v) over all known machines.
22			c where
23			c eps = smallest no. such that eps + 1. .gt. 1.
24			c u = smallest positive no. (underflow limit)
25			c v = largest no. (overflow limit)
26			c
27			c brief outline of algorithm..
28			c
29			c phase 1 scans zero components.
30			c move to phase 2 when a component is nonzero and .le. cutlo
31			c move to phase 3 when a component is .gt. cutlo
32			c move to phase 4 when a component is .ge. cuthi/m
33			c where m = n for x() real and m = 2*n for complex.
34			c
35			c values for cutlo and cuthi..
36			c from the environmental parameters listed in the imsl converter
37			c document the limiting values are as follows..
38			c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are
39			c univac and dec at 2**(-103)
40			c thus cutlo = 2**(-51) = 4.44089e-16
41			c cuthi, s.p. v = 2**127 for univac, honeywell, and dec.
42			c thus cuthi = 2**(63.5) = 1.30438e19
43			c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec.
44			c thus cutlo = 2**(-33.5) = 8.23181d-11
45			c cuthi, d.p. same as s.p. cuthi = 1.30438d19
46			c data cutlo, cuthi / 8.232d-11, 1.304d19 /
47			c data cutlo, cuthi / 4.441e-16, 1.304e19 /
48			data cutlo, cuthi / 8.232d-11, 1.304d19 /
49			c
50			if(n .gt. 0) go to 10
51			dnrm2 = zero
52			go to 300
53			c
54			10 assign 30 to next
55			sum = zero
56			nn = n * incx
57			c begin main loop
58			i = 1
59			20 go to next,(30, 50, 70, 110)
60			30 if( dabs(dx(i)) .gt. cutlo) go to 85
61			assign 50 to next
62			xmax = zero
63			c
64			c phase 1. sum is zero
65			c
66			50 if( dx(i) .eq. zero) go to 200
67			if( dabs(dx(i)) .gt. cutlo) go to 85
68			c
69			c prepare for phase 2.
70			assign 70 to next
71			go to 105
72			c
73			c prepare for phase 4.
74			c
75			100 i = j
76			assign 110 to next
77			sum = (sum / dx(i)) / dx(i)
78			105 xmax = dabs(dx(i))
79			go to 115
80			c
81			c phase 2. sum is small.
82			c scale to avoid destructive underflow.
83			c
84			70 if( dabs(dx(i)) .gt. cutlo ) go to 75
85			c
86			c common code for phases 2 and 4.
87			c in phase 4 sum is large. scale to avoid overflow.
88			c
89			110 if( dabs(dx(i)) .le. xmax ) go to 115
90			sum = one + sum * (xmax / dx(i))**2
91			xmax = dabs(dx(i))
92			go to 200
93			c
94			115 sum = sum + (dx(i)/xmax)**2
95			go to 200
96			c
97			c
98			c prepare for phase 3.
99			c
100			75 sum = (sum * xmax) * xmax
101			c
102			c
103			c for real or d.p. set hitest = cuthi/n
104			c for complex set hitest = cuthi/(2*n)
105			c
106			85 hitest = cuthi/float( n )
107			c
108			c phase 3. sum is mid-range. no scaling.
109			c
110			do 95 j =i,nn,incx
111			if(dabs(dx(j)) .ge. hitest) go to 100
112			95 sum = sum + dx(j)**2
113			dnrm2 = dsqrt( sum )
114			go to 300
115			c
116			200 continue
117			i = i + incx
118			if ( i .le. nn ) go to 20
119			c
120			c end of main loop.
121			c
122			c compute square root and adjust for scaling.
123			c
124			dnrm2 = xmax * dsqrt(sum)
125			300 continue
126			return
127			end
128