/[ascend]/trunk/blas/dnrm2.f
ViewVC logotype

Contents of /trunk/blas/dnrm2.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (15 years, 5 months ago) by aw0a
File size: 3986 byte(s)
Setting up web subdirectory in repository
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

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22