1: !> \brief \b DZNRM2
2: !
3: ! =========== DOCUMENTATION ===========
4: !
5: ! Online html documentation available at
6: ! http://www.netlib.org/lapack/explore-html/
7: !
8: ! Definition:
9: ! ===========
10: !
11: ! DOUBLE PRECISION FUNCTION DZNRM2(N,X,INCX)
12: !
13: ! .. Scalar Arguments ..
14: ! INTEGER INCX,N
15: ! ..
16: ! .. Array Arguments ..
17: ! DOUBLE COMPLEX X(*)
18: ! ..
19: !
20: !
21: !> \par Purpose:
22: ! =============
23: !>
24: !> \verbatim
25: !>
26: !> DZNRM2 returns the euclidean norm of a vector via the function
27: !> name, so that
28: !>
29: !> DZNRM2 := sqrt( x**H*x )
30: !> \endverbatim
31: !
32: ! Arguments:
33: ! ==========
34: !
35: !> \param[in] N
36: !> \verbatim
37: !> N is INTEGER
38: !> number of elements in input vector(s)
39: !> \endverbatim
40: !>
41: !> \param[in] X
42: !> \verbatim
43: !> X is COMPLEX*16 array, dimension (N)
44: !> complex vector with N elements
45: !> \endverbatim
46: !>
47: !> \param[in] INCX
48: !> \verbatim
49: !> INCX is INTEGER, storage spacing between elements of X
50: !> If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
51: !> If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
52: !> If INCX = 0, x isn't a vector so there is no need to call
53: !> this subroutine. If you call it anyway, it will count x(1)
54: !> in the vector norm N times.
55: !> \endverbatim
56: !
57: ! Authors:
58: ! ========
59: !
60: !> \author Edward Anderson, Lockheed Martin
61: !
62: !> \date August 2016
63: !
64: !> \ingroup single_blas_level1
65: !
66: !> \par Contributors:
67: ! ==================
68: !>
69: !> Weslley Pereira, University of Colorado Denver, USA
70: !
71: !> \par Further Details:
72: ! =====================
73: !>
74: !> \verbatim
75: !>
76: !> Anderson E. (2017)
77: !> Algorithm 978: Safe Scaling in the Level 1 BLAS
78: !> ACM Trans Math Softw 44:1--28
79: !> https://doi.org/10.1145/3061665
80: !>
81: !> Blue, James L. (1978)
82: !> A Portable Fortran Program to Find the Euclidean Norm of a Vector
83: !> ACM Trans Math Softw 4:15--23
84: !> https://doi.org/10.1145/355769.355771
85: !>
86: !> \endverbatim
87: !>
88: ! =====================================================================
89: function DZNRM2( n, x, incx )
90: integer, parameter :: wp = kind(1.d0)
91: real(wp) :: DZNRM2
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 scaling 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) - digits(0._wp)) * 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: DZNRM2 = 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: DZNRM2 = scl*sqrt( sumsq )
208: return
209: end function
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dznrm2.f90 CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas