1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
real function snrm2 ( n, sx, incx)
integer i, incx, ix, j, n, next
real sx(1), cutlo, cuthi, hitest, sum, xmax, zero, one
data zero, one /0.0e0, 1.0e0/
c
c euclidean norm of the n-vector stored in sx() 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 modified to correct problem with negative increment, 8/21/90.
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 /
c data cutlo, cuthi / 4.441e-16, 1.304e19 /
c... from Ed Anderson (for Cray)
c data cutlo / 0300315520236314774737b /
c data cuthi / 0500004000000000000000b /
c... from Ed Anderson (for Sun4)
data cutlo / 0.44408921E-15 /
data cuthi / 0.18446743E+20 /
c
if(n .gt. 0) go to 10
snrm2 = zero
go to 300
c
10 assign 30 to next
sum = zero
i = 1
if(incx.lt.0)i = (-n+1)*incx + 1
ix = 1
c begin main loop
20 go to next,(30, 50, 70, 110)
30 if( abs(sx(i)) .gt. cutlo) go to 85
assign 50 to next
xmax = zero
c
c phase 1. sum is zero
c
50 if( sx(i) .eq. zero) go to 200
if( abs(sx(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 continue
assign 110 to next
sum = (sum / sx(i)) / sx(i)
105 xmax = abs(sx(i))
go to 115
c
c phase 2. sum is small.
c scale to avoid destructive underflow.
c
70 if( abs(sx(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( abs(sx(i)) .le. xmax ) go to 115
sum = one + sum * (xmax / sx(i))**2
xmax = abs(sx(i))
go to 200
c
115 sum = sum + (sx(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 = ix, n
if(abs(sx(i)) .ge. hitest) go to 100
sum = sum + sx(i)**2
i = i + incx
95 continue
snrm2 = sqrt( sum )
go to 300
c
200 continue
ix = ix + 1
i = i + incx
if( ix .le. n ) go to 20
c
c end of main loop.
c
c compute square root and adjust for scaling.
c
snrm2 = xmax * sqrt(sum)
300 continue
return
end
|
