![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to snrm2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / view2d / num|
Return to snrm2.f CVS log | Up to [Research Unix] / researchv10no / cmd / view2d / num |
1.1 root 1: real function snrm2 ( n, sx, incx)
2: integer next
3: real sx(1), cutlo, cuthi, hitest, sum, xmax, zero, one
4: data zero, one /0.0e0, 1.0e0/
5: c
6: c euclidean norm of the n-vector stored in sx() with storage
7: c increment incx .
8: c if n .le. 0 return with result = 0.
9: c if n .ge. 1 then incx must be .ge. 1
10: c
11: c c.l.lawson, 1978 jan 08
12: c
13: c four phase method using two built-in constants that are
14: c hopefully applicable to all machines.
15: c cutlo = maximum of sqrt(u/eps) over all known machines.
16: c cuthi = minimum of sqrt(v) over all known machines.
17: c where
18: c eps = smallest no. such that eps + 1. .gt. 1.
19: c u = smallest positive no. (underflow limit)
20: c v = largest no. (overflow limit)
21: c
22: c brief outline of algorithm..
23: c
24: c phase 1 scans zero components.
25: c move to phase 2 when a component is nonzero and .le. cutlo
26: c move to phase 3 when a component is .gt. cutlo
27: c move to phase 4 when a component is .ge. cuthi/m
28: c where m = n for x() real and m = 2*n for complex.
29: c
30: c values for cutlo and cuthi..
31: c from the environmental parameters listed in the imsl converter
32: c document the limiting values are as follows..
33: c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are
34: c univac and dec at 2**(-103)
35: c thus cutlo = 2**(-51) = 4.44089e-16
36: c cuthi, s.p. v = 2**127 for univac, honeywell, and dec.
37: c thus cuthi = 2**(63.5) = 1.30438e19
38: c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec.
39: c thus cutlo = 2**(-33.5) = 8.23181d-11
40: c cuthi, d.p. same as s.p. cuthi = 1.30438d19
41: c data cutlo, cuthi / 8.232d-11, 1.304d19 /
42: c data cutlo, cuthi / 4.441e-16, 1.304e19 /
43: data cutlo, cuthi / 4.441e-16, 1.304e19 /
44: c
45: if(n .gt. 0) go to 10
46: snrm2 = zero
47: go to 300
48: c
49: 10 assign 30 to next
50: sum = zero
51: nn = n * incx
52: c begin main loop
53: i = 1
54: 20 go to next,(30, 50, 70, 110)
55: 30 if( abs(sx(i)) .gt. cutlo) go to 85
56: assign 50 to next
57: xmax = zero
58: c
59: c phase 1. sum is zero
60: c
61: 50 if( sx(i) .eq. zero) go to 200
62: if( abs(sx(i)) .gt. cutlo) go to 85
63: c
64: c prepare for phase 2.
65: assign 70 to next
66: go to 105
67: c
68: c prepare for phase 4.
69: c
70: 100 i = j
71: assign 110 to next
72: sum = (sum / sx(i)) / sx(i)
73: 105 xmax = abs(sx(i))
74: go to 115
75: c
76: c phase 2. sum is small.
77: c scale to avoid destructive underflow.
78: c
79: 70 if( abs(sx(i)) .gt. cutlo ) go to 75
80: c
81: c common code for phases 2 and 4.
82: c in phase 4 sum is large. scale to avoid overflow.
83: c
84: 110 if( abs(sx(i)) .le. xmax ) go to 115
85: sum = one + sum * (xmax / sx(i))**2
86: xmax = abs(sx(i))
87: go to 200
88: c
89: 115 sum = sum + (sx(i)/xmax)**2
90: go to 200
91: c
92: c
93: c prepare for phase 3.
94: c
95: 75 sum = (sum * xmax) * xmax
96: c
97: c
98: c for real or d.p. set hitest = cuthi/n
99: c for complex set hitest = cuthi/(2*n)
100: c
101: 85 hitest = cuthi/float( n )
102: c
103: c phase 3. sum is mid-range. no scaling.
104: c
105: do 95 j =i,nn,incx
106: if(abs(sx(j)) .ge. hitest) go to 100
107: 95 sum = sum + sx(j)**2
108: snrm2 = sqrt( sum )
109: go to 300
110: c
111: 200 continue
112: i = i + incx
113: if ( i .le. nn ) go to 20
114: c
115: c end of main loop.
116: c
117: c compute square root and adjust for scaling.
118: c
119: snrm2 = xmax * sqrt(sum)
120: 300 continue
121: return
122: end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.