Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix division by zero in [z]rotg #4235

Merged
merged 2 commits into from
Oct 9, 2023
Merged

Conversation

angsch
Copy link
Contributor

@angsch angsch commented Sep 20, 2023

Fix division by zero.

With that if condition, the division by da_r corresponds to a division by zero.

if (da_r == ZERO && da_i == ZERO) {
*C = ZERO;
if (db_r == ZERO) {
(*DA) = fabsl(db_i);
*S = *S1 /da_r;
*(S+1) = *(S1+1) /da_r;

Fixed to such that S = conjg(db)/|db|

Remove redundant branch.

With the algorithm change in reference LAPACK, the case da == ZERO is a special branch.

long double ada = fabsl(da);

} else if (ada == ZERO) {
*C = ZERO;
*S = ONE;
*DA = *DB;
*DB = ONE;
} else {

Consequently, there is no longer the need to check.
if (da != ZERO) {

The commit removes the branch.

Add missing check for zero.

} else {
z = ONE / c;

does not have a check c != 0 in the long double code. The float/double version does have this check. The new code is essentially identical to the float/double version

The cases
[ c  s ] * [ 0      ] = [ |db_i| ]
[-s  c ]   [ i*db_i ]   [  0     ]
and
[ c  s ] * [ 0      ] = [ |db_r| ]
[-s  c ]   [ db_r   ]   [  0     ]
computed s incorrectly. To flip the entries of vector,
s should be conjg(db)/|db| and not conjg(db) / da,
where da == 0.0.
* The check da != ZERO is no longer necessary since there
  is a special case ada == ZERO, where ada = |da|.
* Add the missing check c != ZERO before the division.

Note that with these two changes the long double code
follows the float/double version of the code.
@angsch
Copy link
Contributor Author

angsch commented Sep 20, 2023

@vladimir-ch You have commented on rotg before and may want to have a look

@martin-frbg martin-frbg added this to the 0.3.25 milestone Oct 9, 2023
@martin-frbg martin-frbg merged commit 4a0f863 into OpenMathLib:develop Oct 9, 2023
58 checks passed
@vladimir-ch
Copy link
Contributor

Thanks for the notification @angsch. I've checked and our Gonum BLAS already has this fix.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants