Annotation of rpl/lapack/lapack/zlassq.f90, revision 1.1
1.1 ! bertrand 1: !> \brief \b ZLASSQ updates a sum of squares represented in scaled form.
! 2: !
! 3: ! =========== DOCUMENTATION ===========
! 4: !
! 5: ! Online html documentation available at
! 6: ! http://www.netlib.org/lapack/explore-html/
! 7: !
! 8: !> \htmlonly
! 9: !> Download ZLASSQ + dependencies
! 10: !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlassq.f90">
! 11: !> [TGZ]</a>
! 12: !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlassq.f90">
! 13: !> [ZIP]</a>
! 14: !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlassq.f90">
! 15: !> [TXT]</a>
! 16: !> \endhtmlonly
! 17: !
! 18: ! Definition:
! 19: ! ===========
! 20: !
! 21: ! SUBROUTINE ZLASSQ( N, X, INCX, SCALE, SUMSQ )
! 22: !
! 23: ! .. Scalar Arguments ..
! 24: ! INTEGER INCX, N
! 25: ! DOUBLE PRECISION SCALE, SUMSQ
! 26: ! ..
! 27: ! .. Array Arguments ..
! 28: ! DOUBLE COMPLEX X( * )
! 29: ! ..
! 30: !
! 31: !
! 32: !> \par Purpose:
! 33: ! =============
! 34: !>
! 35: !> \verbatim
! 36: !>
! 37: !> ZLASSQ returns the values scl and smsq such that
! 38: !>
! 39: !> ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,
! 40: !>
! 41: !> where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is
! 42: !> assumed to be non-negative.
! 43: !>
! 44: !> scale and sumsq must be supplied in SCALE and SUMSQ and
! 45: !> scl and smsq are overwritten on SCALE and SUMSQ respectively.
! 46: !>
! 47: !> If scale * sqrt( sumsq ) > tbig then
! 48: !> we require: scale >= sqrt( TINY*EPS ) / sbig on entry,
! 49: !> and if 0 < scale * sqrt( sumsq ) < tsml then
! 50: !> we require: scale <= sqrt( HUGE ) / ssml on entry,
! 51: !> where
! 52: !> tbig -- upper threshold for values whose square is representable;
! 53: !> sbig -- scaling constant for big numbers; \see la_constants.f90
! 54: !> tsml -- lower threshold for values whose square is representable;
! 55: !> ssml -- scaling constant for small numbers; \see la_constants.f90
! 56: !> and
! 57: !> TINY*EPS -- tiniest representable number;
! 58: !> HUGE -- biggest representable number.
! 59: !>
! 60: !> \endverbatim
! 61: !
! 62: ! Arguments:
! 63: ! ==========
! 64: !
! 65: !> \param[in] N
! 66: !> \verbatim
! 67: !> N is INTEGER
! 68: !> The number of elements to be used from the vector x.
! 69: !> \endverbatim
! 70: !>
! 71: !> \param[in] X
! 72: !> \verbatim
! 73: !> X is DOUBLE COMPLEX array, dimension (1+(N-1)*abs(INCX))
! 74: !> The vector for which a scaled sum of squares is computed.
! 75: !> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
! 76: !> \endverbatim
! 77: !>
! 78: !> \param[in] INCX
! 79: !> \verbatim
! 80: !> INCX is INTEGER
! 81: !> The increment between successive values of the vector x.
! 82: !> If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
! 83: !> If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
! 84: !> If INCX = 0, x isn't a vector so there is no need to call
! 85: !> this subroutine. If you call it anyway, it will count x(1)
! 86: !> in the vector norm N times.
! 87: !> \endverbatim
! 88: !>
! 89: !> \param[in,out] SCALE
! 90: !> \verbatim
! 91: !> SCALE is DOUBLE PRECISION
! 92: !> On entry, the value scale in the equation above.
! 93: !> On exit, SCALE is overwritten with scl , the scaling factor
! 94: !> for the sum of squares.
! 95: !> \endverbatim
! 96: !>
! 97: !> \param[in,out] SUMSQ
! 98: !> \verbatim
! 99: !> SUMSQ is DOUBLE PRECISION
! 100: !> On entry, the value sumsq in the equation above.
! 101: !> On exit, SUMSQ is overwritten with smsq , the basic sum of
! 102: !> squares from which scl has been factored out.
! 103: !> \endverbatim
! 104: !
! 105: ! Authors:
! 106: ! ========
! 107: !
! 108: !> \author Edward Anderson, Lockheed Martin
! 109: !
! 110: !> \par Contributors:
! 111: ! ==================
! 112: !>
! 113: !> Weslley Pereira, University of Colorado Denver, USA
! 114: !> Nick Papior, Technical University of Denmark, DK
! 115: !
! 116: !> \par Further Details:
! 117: ! =====================
! 118: !>
! 119: !> \verbatim
! 120: !>
! 121: !> Anderson E. (2017)
! 122: !> Algorithm 978: Safe Scaling in the Level 1 BLAS
! 123: !> ACM Trans Math Softw 44:1--28
! 124: !> https://doi.org/10.1145/3061665
! 125: !>
! 126: !> Blue, James L. (1978)
! 127: !> A Portable Fortran Program to Find the Euclidean Norm of a Vector
! 128: !> ACM Trans Math Softw 4:15--23
! 129: !> https://doi.org/10.1145/355769.355771
! 130: !>
! 131: !> \endverbatim
! 132: !
! 133: !> \ingroup OTHERauxiliary
! 134: !
! 135: ! =====================================================================
! 136: subroutine ZLASSQ( n, x, incx, scl, sumsq )
! 137: use LA_CONSTANTS, &
! 138: only: wp=>dp, zero=>dzero, one=>done, &
! 139: sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml
! 140: use LA_XISNAN
! 141: !
! 142: ! -- LAPACK auxiliary routine --
! 143: ! -- LAPACK is a software package provided by Univ. of Tennessee, --
! 144: ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 145: !
! 146: ! .. Scalar Arguments ..
! 147: integer :: incx, n
! 148: real(wp) :: scl, sumsq
! 149: ! ..
! 150: ! .. Array Arguments ..
! 151: complex(wp) :: x(*)
! 152: ! ..
! 153: ! .. Local Scalars ..
! 154: integer :: i, ix
! 155: logical :: notbig
! 156: real(wp) :: abig, amed, asml, ax, ymax, ymin
! 157: ! ..
! 158: !
! 159: ! Quick return if possible
! 160: !
! 161: if( LA_ISNAN(scl) .or. LA_ISNAN(sumsq) ) return
! 162: if( sumsq == zero ) scl = one
! 163: if( scl == zero ) then
! 164: scl = one
! 165: sumsq = zero
! 166: end if
! 167: if (n <= 0) then
! 168: return
! 169: end if
! 170: !
! 171: ! Compute the sum of squares in 3 accumulators:
! 172: ! abig -- sums of squares scaled down to avoid overflow
! 173: ! asml -- sums of squares scaled up to avoid underflow
! 174: ! amed -- sums of squares that do not require scaling
! 175: ! The thresholds and multipliers are
! 176: ! tbig -- values bigger than this are scaled down by sbig
! 177: ! tsml -- values smaller than this are scaled up by ssml
! 178: !
! 179: notbig = .true.
! 180: asml = zero
! 181: amed = zero
! 182: abig = zero
! 183: ix = 1
! 184: if( incx < 0 ) ix = 1 - (n-1)*incx
! 185: do i = 1, n
! 186: ax = abs(real(x(ix)))
! 187: if (ax > tbig) then
! 188: abig = abig + (ax*sbig)**2
! 189: notbig = .false.
! 190: else if (ax < tsml) then
! 191: if (notbig) asml = asml + (ax*ssml)**2
! 192: else
! 193: amed = amed + ax**2
! 194: end if
! 195: ax = abs(aimag(x(ix)))
! 196: if (ax > tbig) then
! 197: abig = abig + (ax*sbig)**2
! 198: notbig = .false.
! 199: else if (ax < tsml) then
! 200: if (notbig) asml = asml + (ax*ssml)**2
! 201: else
! 202: amed = amed + ax**2
! 203: end if
! 204: ix = ix + incx
! 205: end do
! 206: !
! 207: ! Put the existing sum of squares into one of the accumulators
! 208: !
! 209: if( sumsq > zero ) then
! 210: ax = scl*sqrt( sumsq )
! 211: if (ax > tbig) then
! 212: ! We assume scl >= sqrt( TINY*EPS ) / sbig
! 213: abig = abig + (scl*sbig)**2 * sumsq
! 214: else if (ax < tsml) then
! 215: ! We assume scl <= sqrt( HUGE ) / ssml
! 216: if (notbig) asml = asml + (scl*ssml)**2 * sumsq
! 217: else
! 218: amed = amed + scl**2 * sumsq
! 219: end if
! 220: end if
! 221: !
! 222: ! Combine abig and amed or amed and asml if more than one
! 223: ! accumulator was used.
! 224: !
! 225: if (abig > zero) then
! 226: !
! 227: ! Combine abig and amed if abig > 0.
! 228: !
! 229: if (amed > zero .or. LA_ISNAN(amed)) then
! 230: abig = abig + (amed*sbig)*sbig
! 231: end if
! 232: scl = one / sbig
! 233: sumsq = abig
! 234: else if (asml > zero) then
! 235: !
! 236: ! Combine amed and asml if asml > 0.
! 237: !
! 238: if (amed > zero .or. LA_ISNAN(amed)) then
! 239: amed = sqrt(amed)
! 240: asml = sqrt(asml) / ssml
! 241: if (asml > amed) then
! 242: ymin = amed
! 243: ymax = asml
! 244: else
! 245: ymin = asml
! 246: ymax = amed
! 247: end if
! 248: scl = one
! 249: sumsq = ymax**2*( one + (ymin/ymax)**2 )
! 250: else
! 251: scl = one / ssml
! 252: sumsq = asml
! 253: end if
! 254: else
! 255: !
! 256: ! Otherwise all values are mid-range or zero
! 257: !
! 258: scl = one
! 259: sumsq = amed
! 260: end if
! 261: return
! 262: end subroutine
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlassq.f90 CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack