Skip to content

Commit 02ea2d3

Browse files
langoulangou
authored
Merge pull request #909 from weslleyspereira/fix-xlassq
Fix issue #908 related to accumulation in xLASSQ
2 parents 461cc2c + c4f8085 commit 02ea2d3

File tree

4 files changed

+68
-72
lines changed

4 files changed

+68
-72
lines changed

‎SRC/classq.f90

Lines changed: 17 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -44,19 +44,6 @@
4444
!> scale and sumsq must be supplied in SCALE and SUMSQ and
4545
!> scale_out and sumsq_out are overwritten on SCALE and SUMSQ respectively.
4646
!>
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-
!>
6047
!> \endverbatim
6148
!
6249
! Arguments:
@@ -209,13 +196,25 @@ subroutine CLASSQ( n, x, incx, scale, sumsq )
209196
if( sumsq > zero ) then
210197
ax = scale*sqrt( sumsq )
211198
if (ax > tbig) then
212-
! We assume scale >= sqrt( TINY*EPS ) / sbig
213-
abig = abig + (scale*sbig)**2 * sumsq
199+
if (scale > one) then
200+
scale = scale * sbig
201+
abig = abig + scale * (scale * sumsq)
202+
else
203+
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
204+
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
205+
end if
214206
else if (ax < tsml) then
215-
! We assume scale <= sqrt( HUGE ) / ssml
216-
if (notbig) asml = asml + (scale*ssml)**2 * sumsq
207+
if (notbig) then
208+
if (scale < one) then
209+
scale = scale * ssml
210+
asml = asml + scale * (scale * sumsq)
211+
else
212+
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
213+
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
214+
end if
215+
end if
217216
else
218-
amed = amed + scale**2 * sumsq
217+
amed = amed + scale * (scale * sumsq)
219218
end if
220219
end if
221220
!

‎SRC/dlassq.f90

Lines changed: 17 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -44,19 +44,6 @@
4444
!> scale and sumsq must be supplied in SCALE and SUMSQ and
4545
!> scale_out and sumsq_out are overwritten on SCALE and SUMSQ respectively.
4646
!>
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-
!>
6047
!> \endverbatim
6148
!
6249
! Arguments:
@@ -200,13 +187,25 @@ subroutine DLASSQ( n, x, incx, scale, sumsq )
200187
if( sumsq > zero ) then
201188
ax = scale*sqrt( sumsq )
202189
if (ax > tbig) then
203-
! We assume scale >= sqrt( TINY*EPS ) / sbig
204-
abig = abig + (scale*sbig)**2 * sumsq
190+
if (scale > one) then
191+
scale = scale * sbig
192+
abig = abig + scale * (scale * sumsq)
193+
else
194+
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
195+
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
196+
end if
205197
else if (ax < tsml) then
206-
! We assume scale <= sqrt( HUGE ) / ssml
207-
if (notbig) asml = asml + (scale*ssml)**2 * sumsq
198+
if (notbig) then
199+
if (scale < one) then
200+
scale = scale * ssml
201+
asml = asml + scale * (scale * sumsq)
202+
else
203+
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
204+
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
205+
end if
206+
end if
208207
else
209-
amed = amed + scale**2 * sumsq
208+
amed = amed + scale * (scale * sumsq)
210209
end if
211210
end if
212211
!

‎SRC/slassq.f90

Lines changed: 17 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -44,19 +44,6 @@
4444
!> scale and sumsq must be supplied in SCALE and SUMSQ and
4545
!> scale_out and sumsq_out are overwritten on SCALE and SUMSQ respectively.
4646
!>
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-
!>
6047
!> \endverbatim
6148
!
6249
! Arguments:
@@ -200,13 +187,25 @@ subroutine SLASSQ( n, x, incx, scale, sumsq )
200187
if( sumsq > zero ) then
201188
ax = scale*sqrt( sumsq )
202189
if (ax > tbig) then
203-
! We assume scale >= sqrt( TINY*EPS ) / sbig
204-
abig = abig + (scale*sbig)**2 * sumsq
190+
if (scale > one) then
191+
scale = scale * sbig
192+
abig = abig + scale * (scale * sumsq)
193+
else
194+
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
195+
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
196+
end if
205197
else if (ax < tsml) then
206-
! We assume scale <= sqrt( HUGE ) / ssml
207-
if (notbig) asml = asml + (scale*ssml)**2 * sumsq
198+
if (notbig) then
199+
if (scale < one) then
200+
scale = scale * ssml
201+
asml = asml + scale * (scale * sumsq)
202+
else
203+
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
204+
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
205+
end if
206+
end if
208207
else
209-
amed = amed + scale**2 * sumsq
208+
amed = amed + scale * (scale * sumsq)
210209
end if
211210
end if
212211
!

‎SRC/zlassq.f90

Lines changed: 17 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -44,19 +44,6 @@
4444
!> scale and sumsq must be supplied in SCALE and SUMSQ and
4545
!> scale_out and sumsq_out are overwritten on SCALE and SUMSQ respectively.
4646
!>
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-
!>
6047
!> \endverbatim
6148
!
6249
! Arguments:
@@ -209,13 +196,25 @@ subroutine ZLASSQ( n, x, incx, scale, sumsq )
209196
if( sumsq > zero ) then
210197
ax = scale*sqrt( sumsq )
211198
if (ax > tbig) then
212-
! We assume scale >= sqrt( TINY*EPS ) / sbig
213-
abig = abig + (scale*sbig)**2 * sumsq
199+
if (scale > one) then
200+
scale = scale * sbig
201+
abig = abig + scale * (scale * sumsq)
202+
else
203+
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
204+
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
205+
end if
214206
else if (ax < tsml) then
215-
! We assume scale <= sqrt( HUGE ) / ssml
216-
if (notbig) asml = asml + (scale*ssml)**2 * sumsq
207+
if (notbig) then
208+
if (scale < one) then
209+
scale = scale * ssml
210+
asml = asml + scale * (scale * sumsq)
211+
else
212+
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
213+
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
214+
end if
215+
end if
217216
else
218-
amed = amed + scale**2 * sumsq
217+
amed = amed + scale * (scale * sumsq)
219218
end if
220219
end if
221220
!

0 commit comments

Comments
 (0)