LAPACK  3.10.0 LAPACK: Linear Algebra PACKage
scnrm2.f90
Go to the documentation of this file.
1
2 !
3 ! =========== DOCUMENTATION ===========
4 !
5 ! Online html documentation available at
6 ! http://www.netlib.org/lapack/explore-html/
7 !
8 ! Definition:
9 ! ===========
10 !
11 ! REAL FUNCTION SCNRM2(N,X,INCX)
12 !
13 ! .. Scalar Arguments ..
14 ! INTEGER INCX,N
15 ! ..
16 ! .. Array Arguments ..
17 ! COMPLEX X(*)
18 ! ..
19 !
20 !
22 ! =============
31 !
32 ! Arguments:
33 ! ==========
34 !
56 !
57 ! Authors:
58 ! ========
59 !
61 !
63 !
65 !
67 ! ==================
70 !
72 ! =====================
88 ! =====================================================================
89 function scnrm2( n, x, incx )
90  integer, parameter :: wp = kind(1.e0)
91  real(wp) :: scnrm2
92 !
93 ! -- Reference BLAS level1 routine (version 3.9.1) --
94 ! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
95 ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
96 ! March 2021
97 !
98 ! .. Constants ..
99  real(wp), parameter :: zero = 0.0_wp
100  real(wp), parameter :: one = 1.0_wp
101  real(wp), parameter :: maxn = huge(0.0_wp)
102 ! ..
103 ! .. Blue's ccaling constants ..
104  real(wp), parameter :: tsml = real(radix(0._wp), wp)**ceiling( &
105  (minexponent(0._wp) - 1) * 0.5_wp)
106  real(wp), parameter :: tbig = real(radix(0._wp), wp)**floor( &
107  (maxexponent(0._wp) - digits(0._wp) + 1) * 0.5_wp)
108  real(wp), parameter :: ssml = real(radix(0._wp), wp)**( - floor( &
109  (minexponent(0._wp) - 1) * 0.5_wp))
110  real(wp), parameter :: sbig = real(radix(0._wp), wp)**( - ceiling( &
111  (maxexponent(0._wp) - digits(0._wp) + 1) * 0.5_wp))
112 ! ..
113 ! .. Scalar Arguments ..
114  integer :: incx, n
115 ! ..
116 ! .. Array Arguments ..
117  complex(wp) :: x(*)
118 ! ..
119 ! .. Local Scalars ..
120  integer :: i, ix
121  logical :: notbig
122  real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
123 !
124 ! Quick return if possible
125 !
126  scnrm2 = zero
127  if( n <= 0 ) return
128 !
129  scl = one
130  sumsq = zero
131 !
132 ! Compute the sum of squares in 3 accumulators:
133 ! abig -- sums of squares scaled down to avoid overflow
134 ! asml -- sums of squares scaled up to avoid underflow
135 ! amed -- sums of squares that do not require scaling
136 ! The thresholds and multipliers are
137 ! tbig -- values bigger than this are scaled down by sbig
138 ! tsml -- values smaller than this are scaled up by ssml
139 !
140  notbig = .true.
141  asml = zero
142  amed = zero
143  abig = zero
144  ix = 1
145  if( incx < 0 ) ix = 1 - (n-1)*incx
146  do i = 1, n
147  ax = abs(real(x(ix)))
148  if (ax > tbig) then
149  abig = abig + (ax*sbig)**2
150  notbig = .false.
151  else if (ax < tsml) then
152  if (notbig) asml = asml + (ax*ssml)**2
153  else
154  amed = amed + ax**2
155  end if
156  ax = abs(aimag(x(ix)))
157  if (ax > tbig) then
158  abig = abig + (ax*sbig)**2
159  notbig = .false.
160  else if (ax < tsml) then
161  if (notbig) asml = asml + (ax*ssml)**2
162  else
163  amed = amed + ax**2
164  end if
165  ix = ix + incx
166  end do
167 !
168 ! Combine abig and amed or amed and asml if more than one
169 ! accumulator was used.
170 !
171  if (abig > zero) then
172 !
173 ! Combine abig and amed if abig > 0.
174 !
175  if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
176  abig = abig + (amed*sbig)*sbig
177  end if
178  scl = one / sbig
179  sumsq = abig
180  else if (asml > zero) then
181 !
182 ! Combine amed and asml if asml > 0.
183 !
184  if ( (amed > zero) .or. (amed > maxn) .or. (amed /= amed) ) then
185  amed = sqrt(amed)
186  asml = sqrt(asml) / ssml
187  if (asml > amed) then
188  ymin = amed
189  ymax = asml
190  else
191  ymin = asml
192  ymax = amed
193  end if
194  scl = one
195  sumsq = ymax**2*( one + (ymin/ymax)**2 )
196  else
197  scl = one / ssml
198  sumsq = asml
199  end if
200  else
201 !
202 ! Otherwise all values are mid-range
203 !
204  scl = one
205  sumsq = amed
206  end if
207  scnrm2 = scl*sqrt( sumsq )
208  return
209 end function
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition: scnrm2.f90:90