Annotation of rpl/lapack/blas/dnrm2.f90, revision 1.1
1.1 ! bertrand 1: !> \brief \b DNRM2
! 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 DNRM2(N,X,INCX)
! 12: !
! 13: ! .. Scalar Arguments ..
! 14: ! INTEGER INCX,N
! 15: ! ..
! 16: ! .. Array Arguments ..
! 17: ! DOUBLE PRECISION X(*)
! 18: ! ..
! 19: !
! 20: !
! 21: !> \par Purpose:
! 22: ! =============
! 23: !>
! 24: !> \verbatim
! 25: !>
! 26: !> DNRM2 returns the euclidean norm of a vector via the function
! 27: !> name, so that
! 28: !>
! 29: !> DNRM2 := sqrt( x'*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 DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
! 44: !> \endverbatim
! 45: !>
! 46: !> \param[in] INCX
! 47: !> \verbatim
! 48: !> INCX is INTEGER, storage spacing between elements of X
! 49: !> If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
! 50: !> If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
! 51: !> If INCX = 0, x isn't a vector so there is no need to call
! 52: !> this subroutine. If you call it anyway, it will count x(1)
! 53: !> in the vector norm N times.
! 54: !> \endverbatim
! 55: !
! 56: ! Authors:
! 57: ! ========
! 58: !
! 59: !> \author Edward Anderson, Lockheed Martin
! 60: !
! 61: !> \date August 2016
! 62: !
! 63: !> \ingroup single_blas_level1
! 64: !
! 65: !> \par Contributors:
! 66: ! ==================
! 67: !>
! 68: !> Weslley Pereira, University of Colorado Denver, USA
! 69: !
! 70: !> \par Further Details:
! 71: ! =====================
! 72: !>
! 73: !> \verbatim
! 74: !>
! 75: !> Anderson E. (2017)
! 76: !> Algorithm 978: Safe Scaling in the Level 1 BLAS
! 77: !> ACM Trans Math Softw 44:1--28
! 78: !> https://doi.org/10.1145/3061665
! 79: !>
! 80: !> Blue, James L. (1978)
! 81: !> A Portable Fortran Program to Find the Euclidean Norm of a Vector
! 82: !> ACM Trans Math Softw 4:15--23
! 83: !> https://doi.org/10.1145/355769.355771
! 84: !>
! 85: !> \endverbatim
! 86: !>
! 87: ! =====================================================================
! 88: function DNRM2( n, x, incx )
! 89: integer, parameter :: wp = kind(1.d0)
! 90: real(wp) :: DNRM2
! 91: !
! 92: ! -- Reference BLAS level1 routine (version 3.9.1) --
! 93: ! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
! 94: ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 95: ! March 2021
! 96: !
! 97: ! .. Constants ..
! 98: real(wp), parameter :: zero = 0.0_wp
! 99: real(wp), parameter :: one = 1.0_wp
! 100: real(wp), parameter :: maxN = huge(0.0_wp)
! 101: ! ..
! 102: ! .. Blue's scaling constants ..
! 103: real(wp), parameter :: tsml = real(radix(0._wp), wp)**ceiling( &
! 104: (minexponent(0._wp) - 1) * 0.5_wp)
! 105: real(wp), parameter :: tbig = real(radix(0._wp), wp)**floor( &
! 106: (maxexponent(0._wp) - digits(0._wp) + 1) * 0.5_wp)
! 107: real(wp), parameter :: ssml = real(radix(0._wp), wp)**( - floor( &
! 108: (minexponent(0._wp) - digits(0._wp)) * 0.5_wp))
! 109: real(wp), parameter :: sbig = real(radix(0._wp), wp)**( - ceiling( &
! 110: (maxexponent(0._wp) + digits(0._wp) - 1) * 0.5_wp))
! 111: ! ..
! 112: ! .. Scalar Arguments ..
! 113: integer :: incx, n
! 114: ! ..
! 115: ! .. Array Arguments ..
! 116: real(wp) :: x(*)
! 117: ! ..
! 118: ! .. Local Scalars ..
! 119: integer :: i, ix
! 120: logical :: notbig
! 121: real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
! 122: !
! 123: ! Quick return if possible
! 124: !
! 125: DNRM2 = zero
! 126: if( n <= 0 ) return
! 127: !
! 128: scl = one
! 129: sumsq = zero
! 130: !
! 131: ! Compute the sum of squares in 3 accumulators:
! 132: ! abig -- sums of squares scaled down to avoid overflow
! 133: ! asml -- sums of squares scaled up to avoid underflow
! 134: ! amed -- sums of squares that do not require scaling
! 135: ! The thresholds and multipliers are
! 136: ! tbig -- values bigger than this are scaled down by sbig
! 137: ! tsml -- values smaller than this are scaled up by ssml
! 138: !
! 139: notbig = .true.
! 140: asml = zero
! 141: amed = zero
! 142: abig = zero
! 143: ix = 1
! 144: if( incx < 0 ) ix = 1 - (n-1)*incx
! 145: do i = 1, n
! 146: ax = abs(x(ix))
! 147: if (ax > tbig) then
! 148: abig = abig + (ax*sbig)**2
! 149: notbig = .false.
! 150: else if (ax < tsml) then
! 151: if (notbig) asml = asml + (ax*ssml)**2
! 152: else
! 153: amed = amed + ax**2
! 154: end if
! 155: ix = ix + incx
! 156: end do
! 157: !
! 158: ! Combine abig and amed or amed and asml if more than one
! 159: ! accumulator was used.
! 160: !
! 161: if (abig > zero) then
! 162: !
! 163: ! Combine abig and amed if abig > 0.
! 164: !
! 165: if ( (amed > zero) .or. (amed > maxN) .or. (amed /= amed) ) then
! 166: abig = abig + (amed*sbig)*sbig
! 167: end if
! 168: scl = one / sbig
! 169: sumsq = abig
! 170: else if (asml > zero) then
! 171: !
! 172: ! Combine amed and asml if asml > 0.
! 173: !
! 174: if ( (amed > zero) .or. (amed > maxN) .or. (amed /= amed) ) then
! 175: amed = sqrt(amed)
! 176: asml = sqrt(asml) / ssml
! 177: if (asml > amed) then
! 178: ymin = amed
! 179: ymax = asml
! 180: else
! 181: ymin = asml
! 182: ymax = amed
! 183: end if
! 184: scl = one
! 185: sumsq = ymax**2*( one + (ymin/ymax)**2 )
! 186: else
! 187: scl = one / ssml
! 188: sumsq = asml
! 189: end if
! 190: else
! 191: !
! 192: ! Otherwise all values are mid-range
! 193: !
! 194: scl = one
! 195: sumsq = amed
! 196: end if
! 197: DNRM2 = scl*sqrt( sumsq )
! 198: return
! 199: end function
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dnrm2.f90 CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas