floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fix division by zero in zrotg 

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.
This commit is contained in:
Angelika Schwarz 2023-09-20 19:10:08 +02:00
parent bb2f1ec3b0
commit 6876ae0c3b
1 changed files with 5 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
interface/zrotg.c
View File

@ -61,16 +61,16 @@ void CNAME(void *VDA, void *VDB, FLOAT *C, void *VS) {
*(S1 + 0) = *(DB + 0); *(S1 + 0) = *(DB + 0);
*(S1 + 1) = *(DB + 1) *-1; *(S1 + 1) = *(DB + 1) *-1;
if (da_r == ZERO && da_i == ZERO) { if (da_r == ZERO && da_i == ZERO) {
*C = ZERO; *C = ZERO;
if (db_r == ZERO) { if (db_r == ZERO) {
(*DA) = fabsl(db_i); (*DA) = fabsl(db_i);
*S = *S1 /da_r; *S = *S1 /(*DA);
*(S+1) = *(S1+1) /da_r; *(S+1) = *(S1+1) /(*DA);
return; return;
} else if ( db_i == ZERO) { } else if ( db_i == ZERO) {
*DA = fabsl(db_r); *DA = fabsl(db_r);
*S = *S1 /da_r; *S = *S1 /(*DA);
*(S+1) = *(S1+1) /da_r; *(S+1) = *(S1+1) /(*DA);
return; return;
} else { } else {
long double g1 = MAX( fabsl(db_r), fabsl(db_i)); long double g1 = MAX( fabsl(db_r), fabsl(db_i));