rpl/lapack/blas/dnrm2.f90 - view

File: [local] / rpl / lapack / blas / dnrm2.f90
Revision 1.1: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:55:28 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD

Ajout de fichiers de lapack 3.11

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