?LASQ3 could cause NaNs appearance in ?GESVD and ?GESDD

[akobotov]

Hi,

We recently noticed that some of matrices could cause ?GESVD to produce NaNs on output. One example is a matrix with all elements being equal to 1. After bidiagonalization it naturally has a lot of denormalized numbers in it. But not every sequence of such numbers fails, and appearance of a failing sequence affected by matrix sizes, compiler and it’s options used to compile BLAS, and so on. However the rootcause is clear.

NaNs appear in ?LASQ3 on call tree ?GESVD->?BDSQR->?LASQ1->?LASQ2->?LASQ3, because operation on line 271 of ?LASQ3 which is T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) ) could cause T to be equal to ZERO for very close denormalized values in Z after taking the half. Followed operations with use of the T which is equal to ZERO causes NaN appearance: S = something/T and next step in computing S contains S/T.

It looks like the preceding condition Z( NN-5 ).GT.Z( NN-3 )*TOL2 is weak, and doesn’t take into account possibility of T to become ZERO.
Easy fix could be just replacing:
IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 ) THEN
T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
With:
T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) )
IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2 .AND. T.NE.ZERO ) THEN

Thanks,
Alexander

[rodney]

Alexander,

Thanks for finding this. I will put your solution in the svn and it will be in the next release.

Regards,

Rodney