DROTMG returns unexpected values

[vladimir-ch]

In some cases, DROTMG returns values that, when used in DROTM, do not zero out the second component of the given vector.

This can be observed for example with the following input values for DROTMG:

D1 = 1600000000
D2 = 800000000
X1 = 8
Y1 = 7

After DROTM, Y1 is currently mapped to -24052.672896090524.

The issue we see is in the rescaling part. The setup phase sets the non-unit elements of H. When rescaling is done, the unit elements are set explicitly before the scaling. Rescaling for D1 scales H11,H12 and sets FLAG=-1 but then rescaling for D2 overwrites H12 with 1 which is wrong.

In Gonum we fixed the issue by adopting in gonum/gonum#420 the implementation proposed in:

Tim Hopkins. 1997. Restructuring the BLAS level 1 routine for computing the modified
givens transformation. SIGNUM Newsl. 32, 4 (October 1997), 2-14.
DOI=http://dx.doi.org/10.1145/289251.289253

OpenBLAS later fixed their rotmg implementation with the same approach in OpenMathLib/OpenBLAS#1480.
There was also some discussion in OpenMathLib/OpenBLAS#1452 where the initial report about strange values returned by DROTMG was done.

Also, as I comment in that thread, we would like to clarify whether the value of GAM could be increased and to what value. From what I read in the above paper, the current value of 4096 was chosen very conservatively when BLAS Level 1 was created, that is before IEEE floating point. A larger value of GAM would reduce the necessity for rescaling (flag = -1) which is good because having units in H saves operations in DROTM which is the whole point of DROTMG in my understanding. However, it is not clear what an appropriate value would be. What are your views on this?

Lastly, it would worth revising the documentation for DROTMG and thus improve slightly its status of an underdog. It took us some time to fully understand how to interpret the values the routine takes and returns. The most mysterious was the square root of D1 and D2, some projects that provide interface to BLAS even drop the square root from their documentation.

[jacobbogers]

Hi, can anyone give the status on this ussue, is it resolved?

[jacobbogers]

What is the status on this issue?

[langou]

This might be related to #344