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

Annotation of /trunk/blas/dnrm2.f

Parent Directory Parent Directory | Revision Log Revision Log


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

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