dlassq()

subroutine dlassq	(	integer	n,
		real(wp), dimension(*)	x,
		integer	incx,
		real(wp)	scale,
		real(wp)	sumsq
	)

DLASSQ updates a sum of squares represented in scaled form.

Download DLASSQ + dependencies [TGZ] [ZIP] [TXT]

Purpose:

 DLASSQ returns the values scale_out and sumsq_out such that

    (scale_out**2)*sumsq_out = x( 1 )**2 +...+ x( n )**2 + (scale**2)*sumsq,

 where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is
 assumed to be non-negative.

 scale and sumsq must be supplied in SCALE and SUMSQ and
 scale_out and sumsq_out are overwritten on SCALE and SUMSQ respectively.

Parameters

[in]	N	N is INTEGER The number of elements to be used from the vector x.
[in]	X	X is DOUBLE PRECISION array, dimension (1+(N-1)*abs(INCX)) The vector for which a scaled sum of squares is computed. x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
[in]	INCX	INCX is INTEGER The increment between successive values of the vector x. If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n If INCX = 0, x isn't a vector so there is no need to call this subroutine. If you call it anyway, it will count x(1) in the vector norm N times.
[in,out]	SCALE	SCALE is DOUBLE PRECISION On entry, the value scale in the equation above. On exit, SCALE is overwritten by scale_out, the scaling factor for the sum of squares.
[in,out]	SUMSQ	SUMSQ is DOUBLE PRECISION On entry, the value sumsq in the equation above. On exit, SUMSQ is overwritten by sumsq_out, the basic sum of squares from which scale_out has been factored out.

Author

Edward Anderson, Lockheed Martin

Contributors:

Weslley Pereira, University of Colorado Denver, USA Nick Papior, Technical University of Denmark, DK

Further Details:

  Anderson E. (2017)
  Algorithm 978: Safe Scaling in the Level 1 BLAS
  ACM Trans Math Softw 44:1--28
  https://doi.org/10.1145/3061665

  Blue, James L. (1978)
  A Portable Fortran Program to Find the Euclidean Norm of a Vector
  ACM Trans Math Softw 4:15--23
  https://doi.org/10.1145/355769.355771

Definition at line 123 of file dlassq.f90.

   use la_constants, &
      only: wp=>dp, zero=>dzero, one=>done, &
            sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml
   use la_xisnan
!
!  -- LAPACK auxiliary routine --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!
!  .. Scalar Arguments ..
   integer :: incx, n
   real(wp) :: scale, sumsq
!  ..
!  .. Array Arguments ..
   real(wp) :: x(*)
!  ..
!  .. Local Scalars ..
   integer :: i, ix
   logical :: notbig
   real(wp) :: abig, amed, asml, ax, ymax, ymin
!  ..
!
!  Quick return if possible
!
   if( la_isnan(scale) .or. la_isnan(sumsq) ) return
   if( sumsq == zero ) scale = one
   if( scale == zero ) then
      scale = one
      sumsq = zero
   end if
   if (n <= 0) then
      return
   end if
!
!  Compute the sum of squares in 3 accumulators:
!     abig -- sums of squares scaled down to avoid overflow
!     asml -- sums of squares scaled up to avoid underflow
!     amed -- sums of squares that do not require scaling
!  The thresholds and multipliers are
!     tbig -- values bigger than this are scaled down by sbig
!     tsml -- values smaller than this are scaled up by ssml
!
   notbig = .true.
   asml = zero
   amed = zero
   abig = zero
   ix = 1
   if( incx < 0 ) ix = 1 - (n-1)*incx
   do i = 1, n
      ax = abs(x(ix))
      if (ax > tbig) then
         abig = abig + (ax*sbig)**2
         notbig = .false.
      else if (ax < tsml) then
         if (notbig) asml = asml + (ax*ssml)**2
      else
         amed = amed + ax**2
      end if
      ix = ix + incx
   end do
!
!  Put the existing sum of squares into one of the accumulators
!
   if( sumsq > zero ) then
      ax = scale*sqrt( sumsq )
      if (ax > tbig) then
         if (scale > one) then
            scale = scale * sbig
            abig = abig + scale * (scale * sumsq)
         else
            ! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
            abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
         end if
      else if (ax < tsml) then
         if (notbig) then
            if (scale < one) then
               scale = scale * ssml
               asml = asml + scale * (scale * sumsq)
            else
               ! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
               asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
            end if
         end if
      else
         amed = amed + scale * (scale * sumsq)
      end if
   end if
!
!  Combine abig and amed or amed and asml if more than one
!  accumulator was used.
!
   if (abig > zero) then
!
!     Combine abig and amed if abig > 0.
!
      if (amed > zero .or. la_isnan(amed)) then
         abig = abig + (amed*sbig)*sbig
      end if
      scale = one / sbig
      sumsq = abig
   else if (asml > zero) then
!
!     Combine amed and asml if asml > 0.
!
      if (amed > zero .or. la_isnan(amed)) then
         amed = sqrt(amed)
         asml = sqrt(asml) / ssml
         if (asml > amed) then
            ymin = amed
            ymax = asml
         else
            ymin = asml
            ymax = amed
         end if
         scale = one
         sumsq = ymax**2*( one + (ymin/ymax)**2 )
      else
         scale = one / ssml
         sumsq = asml
      end if
   else
!
!     Otherwise all values are mid-range or zero
!
      scale = one
      sumsq = amed
   end if
   return

Here is the caller graph for this function: