Annotation of rpl/lapack/blas/dznrm2.f90, revision 1.1
1.1 ! bertrand 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>