Revert "Merge pull request #290" #570
Conversation
Codecov Report
@@ Coverage Diff @@
## master #570 +/- ##
===========================================
- Coverage 82.37% 0.00% -82.38%
===========================================
Files 1894 1892 -2
Lines 190679 198735 +8056
===========================================
- Hits 157065 0 -157065
- Misses 33614 198735 +165121
Continue to review full report at Codecov.
|
|
Running tests with:
My system:
Results:
Wed May 26 15:16:25 2021
./tester --verbose 1 --type s,d,c,z --align 32 --dim 4000 --norm fro lange
LAPACK++ version 2020.10.02, id 0eb8400
input: ./tester --verbose 1 --type 's,d,c,z' --align 32 --dim 4000 --norm fro lange
LAPACK++ Ref.
type norm m n align error time (s) time (s) status
A m= 4000, n= 4000, lda= 4000
norm_tst = 2.28695190e+03
norm_ref = 2.28695190e+03
test matrix A: rand, cond(S) = NA
s fro 4000 4000 32 0.00e+00 0.04597 0.05621 pass
A m= 4000, n= 4000, lda= 4000
norm_tst = 2.30912303e+03
norm_ref = 2.30912303e+03
d fro 4000 4000 32 0.00e+00 0.04700 0.05795 pass
A m= 4000, n= 4000, lda= 4000
norm_tst = 3.18711377e+03
norm_ref = 3.18711377e+03
c fro 4000 4000 32 0.00e+00 0.09132 0.1071 pass
A m= 4000, n= 4000, lda= 4000
norm_tst = 3.26590201e+03
norm_ref = 3.26590201e+03
z fro 4000 4000 32 0.00e+00 0.09260 0.1143 pass
All tests passed for lange.
pass
All routines passed.
Elapsed 2.58 sec
Wed May 26 15:16:27 2021Using double precision (type d) as reference, this gives a relative error for single precision (type s) of 9.6015e-03.
Wed May 26 16:34:24 2021
./tester --verbose 1 --type s,d,c,z --align 32 --dim 4000 --norm fro lange
LAPACK++ version 2020.10.02, id 0eb8400
input: ./tester --verbose 1 --type 's,d,c,z' --align 32 --dim 4000 --norm fro lange
LAPACK++ Ref.
type norm m n align error time (s) time (s) status
A m= 4000, n= 4000, lda= 4000
norm_tst = 2.30936792e+03
norm_ref = 2.30936792e+03
test matrix A: rand, cond(S) = NA
s fro 4000 4000 32 0.00e+00 0.02936 0.03555 pass
A m= 4000, n= 4000, lda= 4000
norm_tst = 2.30912303e+03
norm_ref = 2.30912303e+03
d fro 4000 4000 32 0.00e+00 0.02205 0.03970 pass
A m= 4000, n= 4000, lda= 4000
norm_tst = 3.26619067e+03
norm_ref = 3.26619067e+03
c fro 4000 4000 32 0.00e+00 0.04240 0.05933 pass
A m= 4000, n= 4000, lda= 4000
norm_tst = 3.26590201e+03
norm_ref = 3.26590201e+03
z fro 4000 4000 32 0.00e+00 0.04327 0.06538 pass
All tests passed for lange.
pass
All routines passed.
Elapsed 2.34 sec
Wed May 26 16:34:26 2021Using double precision (type d) as reference, this gives a relative error for single precision (type s) of 1.0605e-04.
Wed May 26 16:44:58 2021
./tester --verbose 1 --type s,d,c,z --align 32 --dim 4000 --norm fro lange
LAPACK++ version 2020.10.02, id 0eb8400
input: ./tester --verbose 1 --type 's,d,c,z' --align 32 --dim 4000 --norm fro lange
LAPACK++ Ref.
type norm m n align error time (s) time (s) status
A m= 4000, n= 4000, lda= 4000
norm_tst = 2.30937085e+03
norm_ref = 2.30937085e+03
test matrix A: rand, cond(S) = NA
s fro 4000 4000 32 0.00e+00 0.02963 0.03561 pass
A m= 4000, n= 4000, lda= 4000
norm_tst = 2.30912303e+03
norm_ref = 2.30912303e+03
d fro 4000 4000 32 0.00e+00 0.02254 0.03499 pass
A m= 4000, n= 4000, lda= 4000
norm_tst = 3.26619287e+03
norm_ref = 3.26619287e+03
c fro 4000 4000 32 0.00e+00 0.04273 0.06186 pass
A m= 4000, n= 4000, lda= 4000
norm_tst = 3.26590201e+03
norm_ref = 3.26590201e+03
z fro 4000 4000 32 0.00e+00 0.04479 0.06828 pass
All tests passed for lange.
pass
All routines passed.
Elapsed 2.40 sec
Wed May 26 16:45:01 2021Using double precision (type d) as reference, this gives a relative error for single precision (type s) of 1.0732e-04. Summary:
|
|
Hi. I did more tests on my laptop. LAPACK compiled with gfortran Source: // Copyright 2021 Weslley S Pereira
#include <cmath>
#include <limits>
#include <iostream>
#include <bitset>
using Integer = std::int32_t;
extern "C" {
double dlange_(
char* norm, Integer* m, Integer* n,
const double* A, Integer* lda,
double* work, std::size_t norm_len);
}
template<typename T>
void print_binary( const T& a ){
const char* a_char = reinterpret_cast<const char*>(&a);
for (size_t i = 0; i < sizeof(T); ++i)
{
std::bitset<8> x(a_char[sizeof(T)-i-1]);
std::cout << x;
}
std::cout << std::endl;
}
int main() {
using Real = double;
auto normid = 'F';
const Integer N = 3793;
Integer m, n;
Integer *lda = &m;
auto work = std::numeric_limits<Real>::quiet_NaN();
Real norm, expectedNorm, a, *A;
a = pow(2, 1-1023);
std::cout << "a = " << a << " = (memory) ";
print_binary( a );
m = N;
n = N;
A = new Real[m*n];
for (Integer i = 0; i < m*n; ++i) A[i] = a;
expectedNorm = N*a;
std::cout << "expectedNorm = " << expectedNorm << " = (memory) ";
print_binary( expectedNorm );
norm = dlange_(&normid, &m, &n, A, lda, &work, 1);
std::cout << "dlange = " << norm << " = (memory) ";
print_binary( norm );
if(!std::isfinite(norm) || norm != expectedNorm) {
std::fprintf(
stderr, "\nexpected norm %8.2e, got %e\n", expectedNorm, norm);
}
delete[] A;
return 0;
}Results:
Conclusion: Using both the recent master branch bd6add2 and master at c3b03d8 (before #290), we obtain the exact norm. The modifications from this PR (which reverts 8d23489) introduce errors in the computation. I will investigate what is happening. |
!> 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.
weslleyspereira
force-pushed
the
fix-bug-in-xCOMBSSQ
branch
from
bb6d5cf
to
e7d572c
Compare
|
I found one operation with precision loss in the xSLASSQ subroutines. My last commit (bb6d5cf) replaces !
! Put the existing sum of squares into one of the accumulators
!
if( sumsq > zero ) then
ax = scl*sqrt( sumsq )
if (ax > tbig) then
abig = abig + (ax*sbig)**2
notbig = .false.
else if (ax < tsml) then
if (notbig) asml = asml + (ax*ssml)**2
else
amed = amed + ax**2
end if
end ifby !
! Put the existing sum of squares into one of the accumulators
!
if( sumsq > zero ) then
ax = scl*sqrt( sumsq )
if (ax > tbig) then
! We assume scl >= sqrt( TINY*EPS ) / sbig
abig = abig + (scl*sbig)**2 * sumsq
notbig = .false.
else if (ax < tsml) then
! We assume scl <= sqrt( HUGE ) / ssml
if (notbig) asml = asml + (scl*ssml)**2 * sumsq
else
amed = amed + scl**2 * sumsq
end if
end ifSo, I keep the original algorithm but don't compute !> 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.This is reasonable. It basically says that:
In floating-point arithmetic with base
The output pair Current branch at e7d572c: |
|
I can also propose an algorithm that avoids the requirements over I prefer to just ask for |
|
This PR closes #424 |
There was a problem hiding this comment.
Nicely done @weslleyspereira.
(1) I like the fact that we remove xCOMBSSQ, it is better to have a good xLASSQ than to have a bad xLASSQ fixed by xCOMBSSQ on top, now that Ed Anderson improved xLASSQ, it is possible to remove xCOMBSSQ, perfect.
(2) good job for fixing a small inaccuracy problem in xLASSQ, a few more multiplications, but this is worse it for the improved accuracy for this range of numbers.
Sounds great to me.
SRC/dlassq.f90
Outdated
| @@ -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. | |||
There was a problem hiding this comment.
Not an issue but I think that setting notbig here is unnecessary, it's not read again.
There was a problem hiding this comment.
Indeed. Thanks for noticing that!
| @@ -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 | |||
There was a problem hiding this comment.
What happens if scl=1 and sumsq=inf (assuming inf is a valid input for sumsq)? In the previous version, ax would be inf and so would be abig. In the new version, sbig**2 underflows to 0, and times inf gives abig = NaN. I think this is unexpected and undesired.
I've made these changes in our Go port and this is just one of many test cases that started to fail and that is easy to reason about. I haven't looked into the rest, I wanted to share this quickly before the PR is merged.
There was a problem hiding this comment.
What happens if
scl=1andsumsq=inf(assuming inf is a valid input for sumsq)? In the previous version, ax would be inf and so would be abig. In the new version,sbig**2underflows to 0, and times inf gives abig = NaN. I think this is unexpected and undesired.
scl=1 and sumsq=inf would be an invalid input for the proposed change since, in this case, 1 < sqrt( TINY*EPS ) / sbig. If we want to deal if this case, we should add another if-tense inside the scope of if (ax > tbig) then.
Do you think this is an important use case for LASSQ?
There was a problem hiding this comment.
Do you think this is an important use case for LASSQ?
Not at all. I've re-read your comments above about why this change is what it is and now it all clicked. IIUC, the tradeoff is a slightly limited input range for better accuracy. I will adjust our tests accordingly.
There was a problem hiding this comment.
Just for fun. If we want to support all input ranges, the code would be something like:
!
! Put the existing sum of squares into one of the accumulators
!
if( sumsq > zero ) then
ax = scl*sqrt( sumsq )
if (ax > tbig) then
if ( scl >= sqrt( TINY*EPS ) / sbig ) then
abig = abig + (scl*sbig)**2 * sumsq
else if ( scl < tsml ) then
abig = abig + (scl*ssml)**2 * (sbig/ssml)**2 * sumsq
else
abig = abig + (scl**2 * sbig) * (sbig * sumsq)
end if
else if (ax < tsml) then
if (notbig) then
if ( scl <= sqrt( HUGE ) / ssml ) then
...Maybe a little smaller if sqrt( TINY*EPS ) / sbig <= sqrt( HUGE ) / ssml. I didn't check. :)
There was a problem hiding this comment.
The upside of this would be that it would avoid the need to document all the constants and thresholds in the documentation for DLASSQ but I'm not in a position to decide whether this should be done or not. Feel free to re-approve and merge.
| !> 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, |
There was a problem hiding this comment.
Maybe it's nitpicking but documented like this isn't scale=1, sumsq=0 excluded?
Remove Dcombssq and rely on Dlassq to combine scaled sums. See the upstream change in Reference-LAPACK/lapack#570
Remove Dcombssq and rely on Dlassq to combine scaled sums. See the upstream change in Reference-LAPACK/lapack#570
This reverts commit 8d23489, reversing changes made to c3b03d8.
Description
The discussion in #569 put some doubts on the range of valid inputs for (s,d)COMBSSQ. See #569 for some examples where (s,d)COMBSSQ fails. These routines were introduced in #290 to improve the accuracy of the xLANxx subroutines.
This PR removes (s,d)COMBSSQ from the library, which brings us back to the old strategy: use xLASSQ to combine the scaled sums.
Edited: This PR closes #424
Checklist