Break out of potentially infinite rescaling loops in xLARGV, xLARTG and xLARTGP

[martin-frbg]

avoids the hang seen in #411 in the same trivial manner as used for xLARFG in #196 but does not address the related topic of overflows in the caller that led to the appearance of Inf arguments

[thijssteel]

PR looks fine to me. The only thing I would recommend changing is whether this fixes #411. Although this PR will remove the deadlock, it will result in inf values. That issue should remain open untill the overflow has been fixed.

[martin-frbg]

PR comment changed to reflect this, thanks.

[codecov]

Codecov Report

Merging #419 into master will increase coverage by 0.06%.
The diff coverage is 50.00%.

@@            Coverage Diff             @@
##           master     #419      +/-   ##
==========================================
+ Coverage   83.24%   83.31%   +0.06%     
==========================================
  Files        1808     1808              
  Lines      170197   170873     +676     
==========================================
+ Hits       141681   142362     +681     
+ Misses      28516    28511       -5     

Impacted Files	Coverage Δ
SRC/clargv.f	63.51% <0.00%> (ø)
SRC/dlartgp.f	46.00% <0.00%> (ø)
SRC/slartgp.f	46.00% <0.00%> (ø)
SRC/zlargv.f	52.70% <0.00%> (ø)
SRC/clartg.f	96.77% <100.00%> (ø)
SRC/dlartg.f	100.00% <100.00%> (ø)
SRC/slartg.f	100.00% <100.00%> (ø)
SRC/zlartg.f	96.77% <100.00%> (ø)
SRC/dbbcsd.f	94.87% <0.00%> (+1.90%)	⬆️
SRC/cbbcsd.f	94.89% <0.00%> (+1.92%)	⬆️
... and 6 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 8e1e16c...0065263. Read the comment docs.

[thijssteel]

@langou if we're releasing soon, i'd include this PR too. Judging by the amount of times i've personally run into this issue, there are probably a lot of people that suffer from it without even realising.

[langou]

Hi @thijssteel, thanks for helping out in what to release and what not. Much appreciated.
Hi @martin-frbg, thanks for the PR, looks good to me.
Let us wait a few days to see if we have more feedback, and then I'll merge the PR.

[langou]

Quick review of ZLARTG. Throwing some comments here.

1. Not related, but not a big fan of:
   LOG( SAFMIN / EPS ) / LOG( DLAMCH( 'B' ) ) / TWO
   Maybe we could add a couple of parenthesis around:
   ( LOG( SAFMIN / EPS ) / LOG( DLAMCH( 'B' ) ) ) / TWO
   I do not like when we rely on FORTAN using operations from left to right. (Which it does.) I would feel better by putting parenthesis around each / sign. (Personal preference.)

2. There is a similar routine in the BLAS called ZROTG. The comment of ZLARTG reads: "This is a faster version of the BLAS1 routine ZROTG, except for the following differences: (1) F and G are unchanged on return; (2) If G=0, then CS=1 and SN=0; (3) If F=0, then CS=0 and SN is chosen so that R is real." And then we can read "3-5-96 - Modified with a new algorithm by W. Kahan and J. Demmel." And "This version has a few statements commented out for thread safety (machine parameters are computed on each entry). 10 feb 03, SJH". Where SJH is Sven Hammarling.

3. I think the goal of the scaling is to be able to SQUARE numbers without overflow or underflow. However DLAPY2 also takes care of this too. I am not clear on how the scaling of ZLARTG and ZROTG interact together. Maybe reviewing the whole code one day would be needed. Maybe we can call a non-scaling version of DLAPY2 from ZLARTG.

4. Assuming IEEE arithmetic, SAFMIN in DOUBLE is 2^(-1022) which is about 2.2251E-308.
   EPSILON in DOUBLE is 2^(-53) which is about 1.1102E-16.
   SAFMIN2 in DOUBLE is 2^(-484) which is about 2.0021e-146.
   SAFMX2 in DOUBLE is 2^(+484) which is about 4.9948e+145.
   The largest number in DOUBLE is about 2^(1024). The smallest positive number in DOUBLE is 2^(-1022).

5. Assuming IEEE arithmetic, SAFMIN in SINGLE is 2^(-126) which is about 1.1755E-38.
   EPSILON in SINGLE is 2^(-24) which is about 5.9605E-08.
   SAFMIN2 in SINGLE is 2^(-51) which is about 2.2518E+15.
   SAFMX2 in SINGLE is 2^(+51) which is about 2.2518E+15.
   The largest number in SINGLE is about 2^(128). The smallest positive number in DOUBLE is 2^(-126).

6. It could be good to have our test suite testing all this in ways we understand. Maybe a unit test for this particular set of routines is good.

[thijssteel]

@langou I've spent quite a lot of time looking into the scaling. (my research code spends up to 30% of its time in that routine, so it really matters for me). I don't fully understand it, but i can say that the scaling does matter. Disabling it significantly increases the backward error (almost factor 2 in some cases).

[langou]

Hi @thijssteel, Super. Good to know. Thanks for letting us know. By a factor of 2 in backward error, do you mean from from 1e-12 to 2e-12 for example? (Or from 1e-12 to 1e-6?) I am assuming the former. On our end, we are looking at various papers right now to try to understand better. Here is some of the literature on scaling we are looking at. Julien.

• J. L. Blue. 1978. A portable fortran program to find the euclidean norm of a vector. ACM Trans. Math. Softw. 4, 1 (Mar. 1978), 15--23.
• E. Anderson and M. Fahey. 1997. Performance Improvements to LAPACK for the Cray Scientific Library. UT-CS-97-359 (LAPACK Working Note 126). University of Tennessee, Knoxville.
• D. Bindel, J. Demmel, W. Kahan, and O. Marques. 2002, On computing Givens rotations reliably and efficiently, ACM Trans. Math. Softw. 28, 2, 206--238.
• E. Anderson. 2017. Algorithm 978: Safe Scaling in the Level 1 BLAS. ACM Transactions on Mathematical Software. Vol. 44, No. 1, Article 12.

[thijssteel]

Thanks for those references.

It was indeed the first case. In a full run of the double complex QR algorithm (single shifted) the relative backward error went from 3.1E-15 to 4.3E-15 in the worst case, sometimes there is no difference. Small note though, this is an edited version of ZLARTG, because i needed another normalization of the phase. The same problems might not be present in the current implementation.

I haven't had time to read all the references, but if you're worried about the performance impact of the do loop to scale, perhaps it can be replaced with direct calculations using EXPONENT.

[langou]

Hi @thijssteel,

Not worried at all of the performance impact of this DO-LOOP. I think what worries us with this loop is part of the general concerns: [ ranked from major concern, to less major concern ] (1) the potential infinite loop in the presence of NaN, or Infinity, and (2) the creation of NaN while scaling, and (3) the non-propagation of NaN. The performance impact of this loop does not worry me for the least. The iteration COUNT at 20 is fine with me, though it seems that for any input (SINGLE or DOUBLE) the max number for COUNT is at most 3. I am glad that there is hard limit (20 is a fine number) in case things go eery. I would like to understand all this better.

Related: Ed Anderson sent us some codes recently with patches that goes in the same vein.

Note: I do not understand why scaling improves accuracy. This is a tat weird to me. We scale with powers of 2 on purpose so as to keep the base-2 computation the same. But I am not surprised to read that base-2 computation are not exactly the same after a scaling. A few threshold here and there might be absolute and not relative. Might hit denorm faster, etc.

Small note though, this is an edited version of ZLARTG, because i needed another normalization of the phase. The same problems might not be present in the current implementation.

Which normalization do you need? We already have LAPACK/ZLARTG and BLAS/ZROTG. A third one? Argh. Yes, we have a couple of routines to generate Householder reflections for a similar reason. I can see we need a few to generate Givens rotations.

[thijssteel]

In the current routine, the phase of R is taken as the phase of F. I needed it to be G, this also meant that C needed to be complex and S could be real. And then ZROT also needed to be changed... Probably not worth including in the library as its a very very specific case.

I'll try to find out the exact parameters that resulted in the increased error.

[langou]

Thanks for explaining. No need to explain for now more on my end. J.

[langou]

Comment: @weslleyspereira is working on another PR with some codes from Ed Anderson. (Spoiler Alert: There will be a lots of Fortran 90.) We would like to see how the community reacts to these Fortran 90 codes. Ed provided us a lots of routines (all the family the norm computation routines: xLANxx for example, but more). Ed has a new xLARGV so @weslleyspereira will prepare a PR based on this. And we will see what people think about the F90 thing.

[langou]

Hi all,

I think I am going to merge this pull request.

Last round of opinion, polling:
In favor: thumbs up.
Against: thumbs down.
Closing polling in a day.

Cheers,
Julien.

[martin-frbg]

Have to abstain as it is "my" PR, but this trivial fix has been distributed with OpenBLAS since June...

[weslleyspereira]

I am OK with merging the current state of this PR. I just want to let some math registered here. It is about @langou's comment:

Assuming IEEE arithmetic, SAFMIN in DOUBLE is 2^(-1022) which is about 2.2251E-308.
EPSILON in DOUBLE is 2^(-53) which is about 1.1102E-16.
SAFMIN2 in DOUBLE is 2^(-484) which is about 2.0021e-146.
SAFMX2 in DOUBLE is 2^(+484) which is about 4.9948e+145.
The largest number in DOUBLE is about 2^(1024). The smallest positive number in DOUBLE is 2^(-1022).

Assuming IEEE arithmetic, SAFMIN in SINGLE is 2^(-126) which is about 1.1755E-38.
EPSILON in SINGLE is 2^(-24) which is about 5.9605E-08.
SAFMIN2 in SINGLE is 2^(-51) which is about 2.2518E+15.
SAFMX2 in SINGLE is 2^(+51) which is about 2.2518E+15.
The largest number in SINGLE is about 2^(128). The smallest positive number in DOUBLE is 2^(-126).

In the standard IEEE arithmetic:

• If we scale twice, the largest float ~2^(128) turns into ~2^(26) which is behind SAFMX2.
• Using doubles, and again using 2 iterations, the largest double ~2^(1024) turns into ~2^(56) which is also behind SAFMX2.

which means we need only 2 iterations in the loops.

In a generic IEEE float type:

• the largest floating-point number is: HUGE = 2^(e_max) = 2^(3 - e_min), where e_max and e_min are the largest and lowest exponents, respectively.
• and SAFMIN2 = 2^((e_min-1+t)/2), where t is the number of digits in the mantissa.
• then, after doing some math, the loop would need:
  ( 5+t-e_min ) / ( 1-t-e_min )
  iterations to take the number HUGE behind SAFMX2.

I am OK with merging the current state of this PR.