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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
*DECK SNRM2
REAL FUNCTION SNRM2 (N, SX, INCX)
C***BEGIN PROLOGUE SNRM2
C***PURPOSE Compute the Euclidean length (L2 norm) of a vector.
C***LIBRARY SLATEC (BLAS)
C***CATEGORY D1A3B
C***TYPE SINGLE PRECISION (SNRM2-S, DNRM2-D, SCNRM2-C)
C***KEYWORDS BLAS, EUCLIDEAN LENGTH, EUCLIDEAN NORM, L2,
C LINEAR ALGEBRA, UNITARY, VECTOR
C***AUTHOR Lawson, C. L., (JPL)
C Hanson, R. J., (SNLA)
C Kincaid, D. R., (U. of Texas)
C Krogh, F. T., (JPL)
C***DESCRIPTION
C
C B L A S Subprogram
C Description of Parameters
C
C --Input--
C N number of elements in input vector(s)
C SX single precision vector with N elements
C INCX storage spacing between elements of SX
C
C --Output--
C SNRM2 single precision result (zero if N .LE. 0)
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 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
C***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
C Krogh, Basic linear algebra subprograms for Fortran
C usage, Algorithm No. 539, Transactions on Mathematical
C Software 5, 3 (September 1979), pp. 308-323.
C***ROUTINES CALLED (NONE)
C***REVISION HISTORY (YYMMDD)
C 791001 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 890831 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE SNRM2
INTEGER NEXT
REAL SX(*), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE
SAVE CUTLO, CUTHI, ZERO, ONE
DATA ZERO, ONE /0.0E0, 1.0E0/
C
DATA CUTLO, CUTHI /4.441E-16, 1.304E19/
C***FIRST EXECUTABLE STATEMENT SNRM2
IF (N .GT. 0) GO TO 10
SNRM2 = ZERO
GO TO 300
C
10 ASSIGN 30 TO NEXT
SUM = ZERO
NN = N * INCX
C
C BEGIN MAIN LOOP
C
I = 1
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.
C
ASSIGN 70 TO NEXT
GO TO 105
C
C PREPARE FOR PHASE 4.
C
100 I = J
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 PREPARE FOR PHASE 3.
C
75 SUM = (SUM * XMAX) * XMAX
C
C FOR REAL OR D.P. SET HITEST = CUTHI/N
C FOR COMPLEX SET HITEST = CUTHI/(2*N)
C
85 HITEST = CUTHI / N
C
C PHASE 3. SUM IS MID-RANGE. NO SCALING.
C
DO 95 J = I,NN,INCX
IF (ABS(SX(J)) .GE. HITEST) GO TO 100
95 SUM = SUM + SX(J)**2
SNRM2 = SQRT( 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
SNRM2 = XMAX * SQRT(SUM)
300 CONTINUE
RETURN
END
|
