Skip to content

Commit c4f8085

Browse files
Applies fix following @angsch's suggestions
1 parent 78b3aab commit c4f8085

File tree

4 files changed

+28
-116
lines changed

4 files changed

+28
-116
lines changed

SRC/classq.f90

Lines changed: 7 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@
122122
! =====================================================================
123123
subroutine CLASSQ( n, x, incx, scale, sumsq )
124124
use LA_CONSTANTS, &
125-
only: wp=>sp, zero=>szero, one=>sone, safmin=>ssafmin, &
125+
only: wp=>sp, zero=>szero, one=>sone, &
126126
sbig=>ssbig, ssml=>sssml, tbig=>stbig, tsml=>stsml
127127
use LA_XISNAN
128128
!
@@ -140,11 +140,7 @@ subroutine CLASSQ( n, x, incx, scale, sumsq )
140140
! .. Local Scalars ..
141141
integer :: i, ix
142142
logical :: notbig
143-
real(wp) :: abig, amed, asml, ax, ymax, ymin, sqrtmin, sqrtmax
144-
! ..
145-
! .. Set constants ..
146-
sqrtmin = sqrt(safmin)
147-
sqrtmax = one / sqrtmin
143+
real(wp) :: abig, amed, asml, ax, ymax, ymin
148144
! ..
149145
!
150146
! Quick return if possible
@@ -201,42 +197,24 @@ subroutine CLASSQ( n, x, incx, scale, sumsq )
201197
ax = scale*sqrt( sumsq )
202198
if (ax > tbig) then
203199
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
200+
scale = scale * sbig
201+
abig = abig + scale * (scale * sumsq)
212202
else
213203
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
214204
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
215205
end if
216206
else if (ax < tsml) then
217207
if (notbig) then
218208
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
209+
scale = scale * ssml
210+
asml = asml + scale * (scale * sumsq)
227211
else
228212
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
229213
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
230214
end if
231215
end if
232216
else
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
217+
amed = amed + scale * (scale * sumsq)
240218
end if
241219
end if
242220
!

SRC/dlassq.f90

Lines changed: 7 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@
122122
! =====================================================================
123123
subroutine DLASSQ( n, x, incx, scale, sumsq )
124124
use LA_CONSTANTS, &
125-
only: wp=>dp, zero=>dzero, one=>done, safmin=>dsafmin, &
125+
only: wp=>dp, zero=>dzero, one=>done, &
126126
sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml
127127
use LA_XISNAN
128128
!
@@ -140,11 +140,7 @@ subroutine DLASSQ( n, x, incx, scale, sumsq )
140140
! .. Local Scalars ..
141141
integer :: i, ix
142142
logical :: notbig
143-
real(wp) :: abig, amed, asml, ax, ymax, ymin, sqrtmin, sqrtmax
144-
! ..
145-
! .. Set constants ..
146-
sqrtmin = sqrt(safmin)
147-
sqrtmax = one / sqrtmin
143+
real(wp) :: abig, amed, asml, ax, ymax, ymin
148144
! ..
149145
!
150146
! Quick return if possible
@@ -192,42 +188,24 @@ subroutine DLASSQ( n, x, incx, scale, sumsq )
192188
ax = scale*sqrt( sumsq )
193189
if (ax > tbig) then
194190
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
191+
scale = scale * sbig
192+
abig = abig + scale * (scale * sumsq)
203193
else
204194
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
205195
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
206196
end if
207197
else if (ax < tsml) then
208198
if (notbig) then
209199
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
200+
scale = scale * ssml
201+
asml = asml + scale * (scale * sumsq)
218202
else
219203
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
220204
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
221205
end if
222206
end if
223207
else
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
208+
amed = amed + scale * (scale * sumsq)
231209
end if
232210
end if
233211
!

SRC/slassq.f90

Lines changed: 7 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@
122122
! =====================================================================
123123
subroutine SLASSQ( n, x, incx, scale, sumsq )
124124
use LA_CONSTANTS, &
125-
only: wp=>sp, zero=>szero, one=>sone, safmin=>ssafmin, &
125+
only: wp=>sp, zero=>szero, one=>sone, &
126126
sbig=>ssbig, ssml=>sssml, tbig=>stbig, tsml=>stsml
127127
use LA_XISNAN
128128
!
@@ -140,11 +140,7 @@ subroutine SLASSQ( n, x, incx, scale, sumsq )
140140
! .. Local Scalars ..
141141
integer :: i, ix
142142
logical :: notbig
143-
real(wp) :: abig, amed, asml, ax, ymax, ymin, sqrtmin, sqrtmax
144-
! ..
145-
! .. Set constants ..
146-
sqrtmin = sqrt(safmin)
147-
sqrtmax = one / sqrtmin
143+
real(wp) :: abig, amed, asml, ax, ymax, ymin
148144
! ..
149145
!
150146
! Quick return if possible
@@ -192,42 +188,24 @@ subroutine SLASSQ( n, x, incx, scale, sumsq )
192188
ax = scale*sqrt( sumsq )
193189
if (ax > tbig) then
194190
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
191+
scale = scale * sbig
192+
abig = abig + scale * (scale * sumsq)
203193
else
204194
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
205195
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
206196
end if
207197
else if (ax < tsml) then
208198
if (notbig) then
209199
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
200+
scale = scale * ssml
201+
asml = asml + scale * (scale * sumsq)
218202
else
219203
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
220204
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
221205
end if
222206
end if
223207
else
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
208+
amed = amed + scale * (scale * sumsq)
231209
end if
232210
end if
233211
!

SRC/zlassq.f90

Lines changed: 7 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@
122122
! =====================================================================
123123
subroutine ZLASSQ( n, x, incx, scale, sumsq )
124124
use LA_CONSTANTS, &
125-
only: wp=>dp, zero=>dzero, one=>done, safmin=>dsafmin, &
125+
only: wp=>dp, zero=>dzero, one=>done, &
126126
sbig=>dsbig, ssml=>dssml, tbig=>dtbig, tsml=>dtsml
127127
use LA_XISNAN
128128
!
@@ -140,11 +140,7 @@ subroutine ZLASSQ( n, x, incx, scale, sumsq )
140140
! .. Local Scalars ..
141141
integer :: i, ix
142142
logical :: notbig
143-
real(wp) :: abig, amed, asml, ax, ymax, ymin, sqrtmin, sqrtmax
144-
! ..
145-
! .. Set constants ..
146-
sqrtmin = sqrt(safmin)
147-
sqrtmax = one / sqrtmin
143+
real(wp) :: abig, amed, asml, ax, ymax, ymin
148144
! ..
149145
!
150146
! Quick return if possible
@@ -201,42 +197,24 @@ subroutine ZLASSQ( n, x, incx, scale, sumsq )
201197
ax = scale*sqrt( sumsq )
202198
if (ax > tbig) then
203199
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
200+
scale = scale * sbig
201+
abig = abig + scale * (scale * sumsq)
212202
else
213203
! sumsq > tbig^2 => (sbig * (sbig * sumsq)) is representable
214204
abig = abig + scale * (scale * (sbig * (sbig * sumsq)))
215205
end if
216206
else if (ax < tsml) then
217207
if (notbig) then
218208
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
209+
scale = scale * ssml
210+
asml = asml + scale * (scale * sumsq)
227211
else
228212
! sumsq < tsml^2 => (ssml * (ssml * sumsq)) is representable
229213
asml = asml + scale * (scale * (ssml * (ssml * sumsq)))
230214
end if
231215
end if
232216
else
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
217+
amed = amed + scale * (scale * sumsq)
240218
end if
241219
end if
242220
!

0 commit comments

Comments
 (0)