Split test condition in LU computation - non-denormal for computation, exact zero for reporting singularity

This commit is contained in:
Martin Kroeker 2023-03-29 22:14:21 +02:00 committed by GitHub
parent 23f2c4ca5b
commit 6c431239da
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 37 additions and 28 deletions

View File

@ -100,16 +100,21 @@ blasint CNAME(blas_arg_t *args, BLASLONG *range_m, BLASLONG *range_n, FLOAT *sa,
jp--; jp--;
temp1 = *(b + jp); temp1 = *(b + jp);
//if (temp1 != ZERO) { if (temp1 != ZERO) {
#if defined(DOUBLE)
if (fabs(temp1) >= DBL_MIN ) { if (fabs(temp1) >= DBL_MIN ) {
temp1 = dp1 / temp1; #else
if (fabs(temp1) >= FLT_MIN ) {
#endif
temp1 = dp1 / temp1;
if (jp != j) { if (jp != j) {
SWAP_K(j + 1, 0, 0, ZERO, a + j, lda, a + jp, lda, NULL, 0); SWAP_K(j + 1, 0, 0, ZERO, a + j, lda, a + jp, lda, NULL, 0);
} }
if (j + 1 < m) { if (j + 1 < m) {
SCAL_K(m - j - 1, 0, 0, temp1, b + j + 1, 1, NULL, 0, NULL, 0); SCAL_K(m - j - 1, 0, 0, temp1, b + j + 1, 1, NULL, 0, NULL, 0);
} }
}
} else { } else {
if (!info) info = j + 1; if (!info) info = j + 1;
} }

View File

@ -106,30 +106,34 @@ blasint CNAME(blas_arg_t *args, BLASLONG *range_m, BLASLONG *range_n, FLOAT *sa,
temp1 = *(b + jp * 2 + 0); temp1 = *(b + jp * 2 + 0);
temp2 = *(b + jp * 2 + 1); temp2 = *(b + jp * 2 + 1);
// if ((temp1 != ZERO) || (temp2 != ZERO)) { if ((temp1 != ZERO) || (temp2 != ZERO)) {
#if defined(DOUBLE)
if ((fabs(temp1) >= DBL_MIN) || (fabs(temp2) >= DBL_MIN)) { if ((fabs(temp1) >= DBL_MIN) || (fabs(temp2) >= DBL_MIN)) {
#else
if (jp != j) { if ((fabs(temp1) >= FLT_MIN) || (fabs(temp2) >= FLT_MIN)) {
SWAP_K(j + 1, 0, 0, ZERO, ZERO, a + j * 2, lda, #endif
if (jp != j) {
SWAP_K(j + 1, 0, 0, ZERO, ZERO, a + j * 2, lda,
a + jp * 2, lda, NULL, 0); a + jp * 2, lda, NULL, 0);
} }
if (fabs(temp1) >= fabs(temp2)){ if (fabs(temp1) >= fabs(temp2)){
ratio = temp2 / temp1; ratio = temp2 / temp1;
den = dp1 /(temp1 * ( 1 + ratio * ratio)); den = dp1 /(temp1 * ( 1 + ratio * ratio));
temp3 = den; temp3 = den;
temp4 = -ratio * den; temp4 = -ratio * den;
} else { } else {
ratio = temp1 / temp2; ratio = temp1 / temp2;
den = dp1 /(temp2 * ( 1 + ratio * ratio)); den = dp1 /(temp2 * ( 1 + ratio * ratio));
temp3 = ratio * den; temp3 = ratio * den;
temp4 = -den; temp4 = -den;
} }
if (j + 1 < m) { if (j + 1 < m) {
SCAL_K(m - j - 1, 0, 0, temp3, temp4, SCAL_K(m - j - 1, 0, 0, temp3, temp4,
b + (j + 1) * 2, 1, NULL, 0, NULL, 0); b + (j + 1) * 2, 1, NULL, 0, NULL, 0);
} }
}
} else { } else {
if (!info) info = j + 1; if (!info) info = j + 1;
} }