Skip to content

Fix accumulation in LAPACK ?LASSQ (Reference-LAPACK PR 909) #4294

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
martin-frbg merged 1 commit into from
Nov 7, 2023
Merged

Fix accumulation in LAPACK ?LASSQ (Reference-LAPACK PR 909) #4294

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 40 additions & 41 deletions lapack-netlib/SRC/classq.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,28 +34,15 @@
!>
!> \verbatim
!>
!> CLASSQ returns the values scl and smsq such that
!> CLASSQ returns the values scale_out and sumsq_out such that
!>
!> ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,
!> (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
!> 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
!> 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 0 < 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.
!> scale_out and sumsq_out are overwritten on SCALE and SUMSQ respectively.
!>
!> \endverbatim
!
Expand All @@ -72,7 +59,7 @@
!> \verbatim
!> X is COMPLEX 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.
!> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
!> \endverbatim
!>
!> \param[in] INCX
Expand All @@ -82,24 +69,24 @@
!> 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)
!> this subroutine. If you call it anyway, it will count x(1)
!> in the vector norm N times.
!> \endverbatim
!>
!> \param[in,out] SCALE
!> \verbatim
!> SCALE is REAL
!> On entry, the value scale in the equation above.
!> On exit, SCALE is overwritten with scl , the scaling factor
!> 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.
!> \endverbatim
!>
!> \param[in,out] SUMSQ
!> \verbatim
!> SUMSQ is REAL
!> On entry, the value sumsq in the equation above.
!> On exit, SUMSQ is overwritten with smsq , the basic sum of
!> squares from which scl has been factored out.
!> 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.
!> \endverbatim
!
! Authors:
Expand Down Expand Up @@ -130,10 +117,10 @@
!>
!> \endverbatim
!
!> \ingroup OTHERauxiliary
!> \ingroup lassq
!
! =====================================================================
subroutine CLASSQ( n, x, incx, scl, sumsq )
subroutine CLASSQ( n, x, incx, scale, sumsq )
use LA_CONSTANTS, &
only: wp=>sp, zero=>szero, one=>sone, &
sbig=>ssbig, ssml=>sssml, tbig=>stbig, tsml=>stsml
Expand All @@ -145,7 +132,7 @@ subroutine CLASSQ( n, x, incx, scl, sumsq )
!
! .. Scalar Arguments ..
integer :: incx, n
real(wp) :: scl, sumsq
real(wp) :: scale, sumsq
! ..
! .. Array Arguments ..
complex(wp) :: x(*)
Expand All @@ -158,10 +145,10 @@ subroutine CLASSQ( n, x, incx, scl, sumsq )
!
! Quick return if possible
!
if( LA_ISNAN(scl) .or. LA_ISNAN(sumsq) ) return
if( sumsq == zero ) scl = one
if( scl == zero ) then
scl = one
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
Expand Down Expand Up @@ -207,15 +194,27 @@ subroutine CLASSQ( n, x, incx, scl, sumsq )
! Put the existing sum of squares into one of the accumulators
!
if( sumsq > zero ) then
ax = scl*sqrt( sumsq )
ax = scale*sqrt( sumsq )
if (ax > tbig) then
! We assume scl >= sqrt( TINY*EPS ) / sbig
abig = abig + (scl*sbig)**2 * sumsq
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
! We assume scl <= sqrt( HUGE ) / ssml
if (notbig) asml = asml + (scl*ssml)**2 * sumsq
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 + scl**2 * sumsq
amed = amed + scale * (scale * sumsq)
end if
end if
!
Expand All @@ -229,7 +228,7 @@ subroutine CLASSQ( n, x, incx, scl, sumsq )
if (amed > zero .or. LA_ISNAN(amed)) then
abig = abig + (amed*sbig)*sbig
end if
scl = one / sbig
scale = one / sbig
sumsq = abig
else if (asml > zero) then
!
Expand All @@ -245,17 +244,17 @@ subroutine CLASSQ( n, x, incx, scl, sumsq )
ymin = asml
ymax = amed
end if
scl = one
scale = one
sumsq = ymax**2*( one + (ymin/ymax)**2 )
else
scl = one / ssml
scale = one / ssml
sumsq = asml
end if
else
!
! Otherwise all values are mid-range or zero
!
scl = one
scale = one
sumsq = amed
end if
return
Expand Down
81 changes: 40 additions & 41 deletions lapack-netlib/SRC/dlassq.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,28 +34,15 @@
!>
!> \verbatim
!>
!> DLASSQ returns the values scl and smsq such that
!> DLASSQ returns the values scale_out and sumsq_out such that
!>
!> ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,
!> (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
!> 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
!> 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 0 < 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.
!> scale_out and sumsq_out are overwritten on SCALE and SUMSQ respectively.
!>
!> \endverbatim
!
Expand All @@ -72,7 +59,7 @@
!> \verbatim
!> 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.
!> x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
!> \endverbatim
!>
!> \param[in] INCX
Expand All @@ -82,24 +69,24 @@
!> 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)
!> this subroutine. If you call it anyway, it will count x(1)
!> in the vector norm N times.
!> \endverbatim
!>
!> \param[in,out] SCALE
!> \verbatim
!> SCALE is DOUBLE PRECISION
!> On entry, the value scale in the equation above.
!> On exit, SCALE is overwritten with scl , the scaling factor
!> 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.
!> \endverbatim
!>
!> \param[in,out] SUMSQ
!> \verbatim
!> SUMSQ is DOUBLE PRECISION
!> On entry, the value sumsq in the equation above.
!> On exit, SUMSQ is overwritten with smsq , the basic sum of
!> squares from which scl has been factored out.
!> 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.
!> \endverbatim
!
! Authors:
Expand Down Expand Up @@ -130,10 +117,10 @@
!>
!> \endverbatim
!
!> \ingroup OTHERauxiliary
!> \ingroup lassq
!
! =====================================================================
subroutine DLASSQ( n, x, incx, scl, sumsq )
subroutine DLASSQ( n, x, incx, scale, sumsq )
use LA_CONSTANTS, &
only: wp=>dp, zero=>dzero, one=>done, &
sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml
Expand All @@ -145,7 +132,7 @@ subroutine DLASSQ( n, x, incx, scl, sumsq )
!
! .. Scalar Arguments ..
integer :: incx, n
real(wp) :: scl, sumsq
real(wp) :: scale, sumsq
! ..
! .. Array Arguments ..
real(wp) :: x(*)
Expand All @@ -158,10 +145,10 @@ subroutine DLASSQ( n, x, incx, scl, sumsq )
!
! Quick return if possible
!
if( LA_ISNAN(scl) .or. LA_ISNAN(sumsq) ) return
if( sumsq == zero ) scl = one
if( scl == zero ) then
scl = one
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
Expand Down Expand Up @@ -198,15 +185,27 @@ subroutine DLASSQ( n, x, incx, scl, sumsq )
! Put the existing sum of squares into one of the accumulators
!
if( sumsq > zero ) then
ax = scl*sqrt( sumsq )
ax = scale*sqrt( sumsq )
if (ax > tbig) then
! We assume scl >= sqrt( TINY*EPS ) / sbig
abig = abig + (scl*sbig)**2 * sumsq
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
! We assume scl <= sqrt( HUGE ) / ssml
if (notbig) asml = asml + (scl*ssml)**2 * sumsq
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 + scl**2 * sumsq
amed = amed + scale * (scale * sumsq)
end if
end if
!
Expand All @@ -220,7 +219,7 @@ subroutine DLASSQ( n, x, incx, scl, sumsq )
if (amed > zero .or. LA_ISNAN(amed)) then
abig = abig + (amed*sbig)*sbig
end if
scl = one / sbig
scale = one / sbig
sumsq = abig
else if (asml > zero) then
!
Expand All @@ -236,17 +235,17 @@ subroutine DLASSQ( n, x, incx, scl, sumsq )
ymin = asml
ymax = amed
end if
scl = one
scale = one
sumsq = ymax**2*( one + (ymin/ymax)**2 )
else
scl = one / ssml
scale = one / ssml
sumsq = asml
end if
else
!
! Otherwise all values are mid-range or zero
!
scl = one
scale = one
sumsq = amed
end if
return
Expand Down
Loading