Skip to content

Commit 78b3aab

Browse files
weslleyspereiraweslleyspereira
committed
fix issue #908 
1 parent ae2f14f commit 78b3aab

File tree

4 files changed

+164
-80
lines changed

4 files changed

+164
-80
lines changed

‎SRC/classq.f90

Lines changed: 41 additions & 20 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:
@@ -135,7 +122,7 @@
135122
! =====================================================================
136123
subroutine CLASSQ( n, x, incx, scale, sumsq )
137124
use LA_CONSTANTS, &
138-
only: wp=>sp, zero=>szero, one=>sone, &
125+
only: wp=>sp, zero=>szero, one=>sone, safmin=>ssafmin, &
139126
sbig=>ssbig, ssml=>sssml, tbig=>stbig, tsml=>stsml
140127
use LA_XISNAN
141128
!
@@ -153,7 +140,11 @@ subroutine CLASSQ( n, x, incx, scale, sumsq )
153140
! .. Local Scalars ..
154141
integer :: i, ix
155142
logical :: notbig
156-
real(wp) :: abig, amed, asml, ax, ymax, ymin
143+
real(wp) :: abig, amed, asml, ax, ymax, ymin, sqrtmin, sqrtmax
144+
! ..
145+
! .. Set constants ..
146+
sqrtmin = sqrt(safmin)
147+
sqrtmax = one / sqrtmin
157148
! ..
158149
!
159150
! Quick return if possible
@@ -209,13 +200,43 @@ subroutine CLASSQ( n, x, incx, scale, sumsq )
209200
if( sumsq > zero ) then
210201
ax = scale*sqrt( sumsq )
211202
if (ax > tbig) then
212-
! We assume scale >= sqrt( TINY*EPS ) / sbig
213-
abig = abig + (scale*sbig)**2 * sumsq
203+
if (scale > one) then
204+
scale = scale * sbig ! sbig < scale <= sbig * max
205+
if (scale > sqrtmin) then
206+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
207+
abig = abig + (scale * scale) * sumsq
208+
else
209+
! Do not square scale, as it may underflow
210+
abig = abig + scale * (scale * sumsq)
211+
end if
212+
else
213+
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
214+
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
215+
end if
214216
else if (ax < tsml) then
215-
! We assume scale <= sqrt( HUGE ) / ssml
216-
if (notbig) asml = asml + (scale*ssml)**2 * sumsq
217+
if (notbig) then
218+
if (scale < one) then
219+
scale = scale * ssml ! ssml * min <= scale < ssml
220+
if (scale < sqrtmax) then
221+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
222+
asml = asml + (scale * scale) * sumsq
223+
else
224+
! Do not square scale, as it may overflow
225+
asml = asml + scale * (scale * sumsq)
226+
end if
227+
else
228+
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
229+
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
230+
end if
231+
end if
217232
else
218-
amed = amed + scale**2 * sumsq
233+
if (scale > sqrtmin .and. scale < sqrtmax) then
234+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
235+
amed = amed + (scale * scale) * sumsq
236+
else
237+
! Do not square scale, as it may overflow
238+
amed = amed + scale * (scale * sumsq)
239+
end if
219240
end if
220241
end if
221242
!

‎SRC/dlassq.f90

Lines changed: 41 additions & 20 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:
@@ -135,7 +122,7 @@
135122
! =====================================================================
136123
subroutine DLASSQ( n, x, incx, scale, sumsq )
137124
use LA_CONSTANTS, &
138-
only: wp=>dp, zero=>dzero, one=>done, &
125+
only: wp=>dp, zero=>dzero, one=>done, safmin=>dsafmin, &
139126
sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml
140127
use LA_XISNAN
141128
!
@@ -153,7 +140,11 @@ subroutine DLASSQ( n, x, incx, scale, sumsq )
153140
! .. Local Scalars ..
154141
integer :: i, ix
155142
logical :: notbig
156-
real(wp) :: abig, amed, asml, ax, ymax, ymin
143+
real(wp) :: abig, amed, asml, ax, ymax, ymin, sqrtmin, sqrtmax
144+
! ..
145+
! .. Set constants ..
146+
sqrtmin = sqrt(safmin)
147+
sqrtmax = one / sqrtmin
157148
! ..
158149
!
159150
! Quick return if possible
@@ -200,13 +191,43 @@ subroutine DLASSQ( n, x, incx, scale, sumsq )
200191
if( sumsq > zero ) then
201192
ax = scale*sqrt( sumsq )
202193
if (ax > tbig) then
203-
! We assume scale >= sqrt( TINY*EPS ) / sbig
204-
abig = abig + (scale*sbig)**2 * sumsq
194+
if (scale > one) then
195+
scale = scale * sbig ! sbig < scale <= sbig * max
196+
if (scale > sqrtmin) then
197+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
198+
abig = abig + (scale * scale) * sumsq
199+
else
200+
! Do not square scale, as it may underflow
201+
abig = abig + scale * (scale * sumsq)
202+
end if
203+
else
204+
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
205+
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
206+
end if
205207
else if (ax < tsml) then
206-
! We assume scale <= sqrt( HUGE ) / ssml
207-
if (notbig) asml = asml + (scale*ssml)**2 * sumsq
208+
if (notbig) then
209+
if (scale < one) then
210+
scale = scale * ssml ! ssml * min <= scale < ssml
211+
if (scale < sqrtmax) then
212+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
213+
asml = asml + (scale * scale) * sumsq
214+
else
215+
! Do not square scale, as it may overflow
216+
asml = asml + scale * (scale * sumsq)
217+
end if
218+
else
219+
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
220+
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
221+
end if
222+
end if
208223
else
209-
amed = amed + scale**2 * sumsq
224+
if (scale > sqrtmin .and. scale < sqrtmax) then
225+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
226+
amed = amed + (scale * scale) * sumsq
227+
else
228+
! Do not square scale, as it may overflow
229+
amed = amed + scale * (scale * sumsq)
230+
end if
210231
end if
211232
end if
212233
!

‎SRC/slassq.f90

Lines changed: 41 additions & 20 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:
@@ -135,7 +122,7 @@
135122
! =====================================================================
136123
subroutine SLASSQ( n, x, incx, scale, sumsq )
137124
use LA_CONSTANTS, &
138-
only: wp=>sp, zero=>szero, one=>sone, &
125+
only: wp=>sp, zero=>szero, one=>sone, safmin=>ssafmin, &
139126
sbig=>ssbig, ssml=>sssml, tbig=>stbig, tsml=>stsml
140127
use LA_XISNAN
141128
!
@@ -153,7 +140,11 @@ subroutine SLASSQ( n, x, incx, scale, sumsq )
153140
! .. Local Scalars ..
154141
integer :: i, ix
155142
logical :: notbig
156-
real(wp) :: abig, amed, asml, ax, ymax, ymin
143+
real(wp) :: abig, amed, asml, ax, ymax, ymin, sqrtmin, sqrtmax
144+
! ..
145+
! .. Set constants ..
146+
sqrtmin = sqrt(safmin)
147+
sqrtmax = one / sqrtmin
157148
! ..
158149
!
159150
! Quick return if possible
@@ -200,13 +191,43 @@ subroutine SLASSQ( n, x, incx, scale, sumsq )
200191
if( sumsq > zero ) then
201192
ax = scale*sqrt( sumsq )
202193
if (ax > tbig) then
203-
! We assume scale >= sqrt( TINY*EPS ) / sbig
204-
abig = abig + (scale*sbig)**2 * sumsq
194+
if (scale > one) then
195+
scale = scale * sbig ! sbig < scale <= sbig * max
196+
if (scale > sqrtmin) then
197+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
198+
abig = abig + (scale * scale) * sumsq
199+
else
200+
! Do not square scale, as it may underflow
201+
abig = abig + scale * (scale * sumsq)
202+
end if
203+
else
204+
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
205+
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
206+
end if
205207
else if (ax < tsml) then
206-
! We assume scale <= sqrt( HUGE ) / ssml
207-
if (notbig) asml = asml + (scale*ssml)**2 * sumsq
208+
if (notbig) then
209+
if (scale < one) then
210+
scale = scale * ssml ! ssml * min <= scale < ssml
211+
if (scale < sqrtmax) then
212+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
213+
asml = asml + (scale * scale) * sumsq
214+
else
215+
! Do not square scale, as it may overflow
216+
asml = asml + scale * (scale * sumsq)
217+
end if
218+
else
219+
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
220+
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
221+
end if
222+
end if
208223
else
209-
amed = amed + scale**2 * sumsq
224+
if (scale > sqrtmin .and. scale < sqrtmax) then
225+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
226+
amed = amed + (scale * scale) * sumsq
227+
else
228+
! Do not square scale, as it may overflow
229+
amed = amed + scale * (scale * sumsq)
230+
end if
210231
end if
211232
end if
212233
!

‎SRC/zlassq.f90

Lines changed: 41 additions & 20 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:
@@ -135,7 +122,7 @@
135122
! =====================================================================
136123
subroutine ZLASSQ( n, x, incx, scale, sumsq )
137124
use LA_CONSTANTS, &
138-
only: wp=>dp, zero=>dzero, one=>done, &
125+
only: wp=>dp, zero=>dzero, one=>done, safmin=>dsafmin, &
139126
sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml
140127
use LA_XISNAN
141128
!
@@ -153,7 +140,11 @@ subroutine ZLASSQ( n, x, incx, scale, sumsq )
153140
! .. Local Scalars ..
154141
integer :: i, ix
155142
logical :: notbig
156-
real(wp) :: abig, amed, asml, ax, ymax, ymin
143+
real(wp) :: abig, amed, asml, ax, ymax, ymin, sqrtmin, sqrtmax
144+
! ..
145+
! .. Set constants ..
146+
sqrtmin = sqrt(safmin)
147+
sqrtmax = one / sqrtmin
157148
! ..
158149
!
159150
! Quick return if possible
@@ -209,13 +200,43 @@ subroutine ZLASSQ( n, x, incx, scale, sumsq )
209200
if( sumsq > zero ) then
210201
ax = scale*sqrt( sumsq )
211202
if (ax > tbig) then
212-
! We assume scale >= sqrt( TINY*EPS ) / sbig
213-
abig = abig + (scale*sbig)**2 * sumsq
203+
if (scale > one) then
204+
scale = scale * sbig ! sbig < scale <= sbig * max
205+
if (scale > sqrtmin) then
206+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
207+
abig = abig + (scale * scale) * sumsq
208+
else
209+
! Do not square scale, as it may underflow
210+
abig = abig + scale * (scale * sumsq)
211+
end if
212+
else
213+
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
214+
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
215+
end if
214216
else if (ax < tsml) then
215-
! We assume scale <= sqrt( HUGE ) / ssml
216-
if (notbig) asml = asml + (scale*ssml)**2 * sumsq
217+
if (notbig) then
218+
if (scale < one) then
219+
scale = scale * ssml ! ssml * min <= scale < ssml
220+
if (scale < sqrtmax) then
221+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
222+
asml = asml + (scale * scale) * sumsq
223+
else
224+
! Do not square scale, as it may overflow
225+
asml = asml + scale * (scale * sumsq)
226+
end if
227+
else
228+
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
229+
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
230+
end if
231+
end if
217232
else
218-
amed = amed + scale**2 * sumsq
233+
if (scale > sqrtmin .and. scale < sqrtmax) then
234+
! sqrtmin < scale < sqrtmax, so it is safe to square scale
235+
amed = amed + (scale * scale) * sumsq
236+
else
237+
! Do not square scale, as it may overflow
238+
amed = amed + scale * (scale * sumsq)
239+
end if
219240
end if
220241
end if
221242
!

0 commit comments

Comments
 (0)