Improve xSLASSQ by removing one operation like sqrt( sumsq )**2

!> If scale * sqrt( sumsq ) > tbig then
!>    we now require:   scale >= sqrt( TINY*EPS ) / sbig   on entry,
!> and if scale * sqrt( sumsq ) < tsml then
!>    we now require:   scale <= sqrt( HUGE ) / ssml       on entry,
!> where
!>    tbig -- upper threshold for values whose square is representable;
!>    sbig -- scaling constant for big numbers; \see la_constants.f90
!>    tsml -- lower threshold for values whose square is representable;
!>    ssml -- scaling constant for small numbers; \see la_constants.f90
!> and
!>    TINY*EPS -- tiniest representable number;
!>    HUGE     -- biggest representable number.

Author: weslleyspereira

SRC/classq.f90

@@ -39,13 +39,24 @@
 !>    ( scl**2 )*smsq = 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 and  scl  returns the value
-!>
-!>    scl = max( scale, abs( x( i ) ) ).
+!> assumed to be non-negative.
 !>
 !> scale and sumsq must be supplied in SCALE and SUMSQ and
 !> scl and smsq are overwritten on SCALE and SUMSQ respectively.
 !>
+!> If scale * sqrt( sumsq ) > tbig then
+!>    we require:   scale >= sqrt( TINY*EPS ) / sbig   on entry,
+!> and if scale * sqrt( sumsq ) < tsml then
+!>    we require:   scale <= sqrt( HUGE ) / ssml       on entry,
+!> where
+!>    tbig -- upper threshold for values whose square is representable;
+!>    sbig -- scaling constant for big numbers; \see la_constants.f90
+!>    tsml -- lower threshold for values whose square is representable;
+!>    ssml -- scaling constant for small numbers; \see la_constants.f90
+!> and
+!>    TINY*EPS -- tiniest representable number;
+!>    HUGE     -- biggest representable number.
+!>
 !> \endverbatim
 !
 !  Arguments:
@@ -198,12 +209,14 @@ subroutine CLASSQ( n, x, incx, scl, sumsq )
    if( sumsq > zero ) then
       ax = scl*sqrt( sumsq )
       if (ax > tbig) then
-         abig = abig + (ax*sbig)**2
+!        We assume scl >= sqrt( TINY*EPS ) / sbig
+         abig = abig + (scl*sbig)**2 * sumsq
          notbig = .false.
       else if (ax < tsml) then
-         if (notbig) asml = asml + (ax*ssml)**2
+!        We assume scl <= sqrt( HUGE ) / ssml
+         if (notbig) asml = asml + (scl*ssml)**2 * sumsq
       else
-         amed = amed + ax**2
+         amed = amed + scl**2 * sumsq
       end if
    end if
 !

SRC/dlassq.f90

@@ -39,13 +39,24 @@
 !>    ( scl**2 )*smsq = 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 and  scl  returns the value
-!>
-!>    scl = max( scale, abs( x( i ) ) ).
+!> assumed to be non-negative.
 !>
 !> scale and sumsq must be supplied in SCALE and SUMSQ and
 !> scl and smsq are overwritten on SCALE and SUMSQ respectively.
 !>
+!> If scale * sqrt( sumsq ) > tbig then
+!>    we require:   scale >= sqrt( TINY*EPS ) / sbig   on entry,
+!> and if scale * sqrt( sumsq ) < tsml then
+!>    we require:   scale <= sqrt( HUGE ) / ssml       on entry,
+!> where
+!>    tbig -- upper threshold for values whose square is representable;
+!>    sbig -- scaling constant for big numbers; \see la_constants.f90
+!>    tsml -- lower threshold for values whose square is representable;
+!>    ssml -- scaling constant for small numbers; \see la_constants.f90
+!> and
+!>    TINY*EPS -- tiniest representable number;
+!>    HUGE     -- biggest representable number.
+!>
 !> \endverbatim
 !
 !  Arguments:
@@ -189,12 +200,14 @@ subroutine DLASSQ( n, x, incx, scl, sumsq )
    if( sumsq > zero ) then
       ax = scl*sqrt( sumsq )
       if (ax > tbig) then
-         abig = abig + (ax*sbig)**2
+!        We assume scl >= sqrt( TINY*EPS ) / sbig
+         abig = abig + (scl*sbig)**2 * sumsq
          notbig = .false.
       else if (ax < tsml) then
-         if (notbig) asml = asml + (ax*ssml)**2
+!        We assume scl <= sqrt( HUGE ) / ssml
+         if (notbig) asml = asml + (scl*ssml)**2 * sumsq
       else
-         amed = amed + ax**2
+         amed = amed + scl**2 * sumsq
       end if
    end if
 !

SRC/slassq.f90

@@ -39,13 +39,24 @@
 !>    ( scl**2 )*smsq = 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 and  scl  returns the value
-!>
-!>    scl = max( scale, abs( x( i ) ) ).
+!> assumed to be non-negative.
 !>
 !> scale and sumsq must be supplied in SCALE and SUMSQ and
 !> scl and smsq are overwritten on SCALE and SUMSQ respectively.
 !>
+!> If scale * sqrt( sumsq ) > tbig then
+!>    we require:   scale >= sqrt( TINY*EPS ) / sbig   on entry,
+!> and if scale * sqrt( sumsq ) < tsml then
+!>    we require:   scale <= sqrt( HUGE ) / ssml       on entry,
+!> where
+!>    tbig -- upper threshold for values whose square is representable;
+!>    sbig -- scaling constant for big numbers; \see la_constants.f90
+!>    tsml -- lower threshold for values whose square is representable;
+!>    ssml -- scaling constant for small numbers; \see la_constants.f90
+!> and
+!>    TINY*EPS -- tiniest representable number;
+!>    HUGE     -- biggest representable number.
+!>
 !> \endverbatim
 !
 !  Arguments:
@@ -189,12 +200,14 @@ subroutine SLASSQ( n, x, incx, scl, sumsq )
    if( sumsq > zero ) then
       ax = scl*sqrt( sumsq )
       if (ax > tbig) then
-         abig = abig + (ax*sbig)**2
+!        We assume scl >= sqrt( TINY*EPS ) / sbig
+         abig = abig + (scl*sbig)**2 * sumsq
          notbig = .false.
       else if (ax < tsml) then
-         if (notbig) asml = asml + (ax*ssml)**2
+!        We assume scl <= sqrt( HUGE ) / ssml
+         if (notbig) asml = asml + (scl*ssml)**2 * sumsq
       else
-         amed = amed + ax**2
+         amed = amed + scl**2 * sumsq
       end if
    end if
 !

SRC/zlassq.f90

@@ -39,13 +39,24 @@
 !>    ( scl**2 )*smsq = 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 and  scl  returns the value
-!>
-!>    scl = max( scale, abs( x( i ) ) ).
+!> assumed to be non-negative.
 !>
 !> scale and sumsq must be supplied in SCALE and SUMSQ and
 !> scl and smsq are overwritten on SCALE and SUMSQ respectively.
 !>
+!> If scale * sqrt( sumsq ) > tbig then
+!>    we require:   scale >= sqrt( TINY*EPS ) / sbig   on entry,
+!> and if scale * sqrt( sumsq ) < tsml then
+!>    we require:   scale <= sqrt( HUGE ) / ssml       on entry,
+!> where
+!>    tbig -- upper threshold for values whose square is representable;
+!>    sbig -- scaling constant for big numbers; \see la_constants.f90
+!>    tsml -- lower threshold for values whose square is representable;
+!>    ssml -- scaling constant for small numbers; \see la_constants.f90
+!> and
+!>    TINY*EPS -- tiniest representable number;
+!>    HUGE     -- biggest representable number.
+!>
 !> \endverbatim
 !
 !  Arguments:
@@ -198,12 +209,14 @@ subroutine ZLASSQ( n, x, incx, scl, sumsq )
    if( sumsq > zero ) then
       ax = scl*sqrt( sumsq )
       if (ax > tbig) then
-         abig = abig + (ax*sbig)**2
+!        We assume scl >= sqrt( TINY*EPS ) / sbig
+         abig = abig + (scl*sbig)**2 * sumsq
          notbig = .false.
       else if (ax < tsml) then
-         if (notbig) asml = asml + (ax*ssml)**2
+!        We assume scl <= sqrt( HUGE ) / ssml
+         if (notbig) asml = asml + (scl*ssml)**2 * sumsq
       else
-         amed = amed + ax**2
+         amed = amed + scl**2 * sumsq
       end if
    end if
 !