From 09abbd017b2c50b86bc0ee9ad971bc8bddbad74f Mon Sep 17 00:00:00 2001 From: Weslley S Pereira Date: Thu, 2 Nov 2023 11:02:19 -0600 Subject: [PATCH 1/2] GECON returns 1 if RCOND is NaN --- SRC/cgecon.f | 13 +++++++++++-- SRC/dgecon.f | 13 +++++++++++-- SRC/sgecon.f | 13 +++++++++++-- SRC/zgecon.f | 13 +++++++++++-- 4 files changed, 44 insertions(+), 8 deletions(-) diff --git a/SRC/cgecon.f b/SRC/cgecon.f index af3a50c6a1..ed705f796a 100644 --- a/SRC/cgecon.f +++ b/SRC/cgecon.f @@ -106,7 +106,7 @@ *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value -*> =-5: if ANORM is NAN or negative. +*> = 1: RCOND = NaN. *> \endverbatim * * Authors: @@ -183,7 +183,7 @@ SUBROUTINE CGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 - ELSE IF( ANORM.LT.ZERO .OR. SISNAN( ANORM ) ) THEN + ELSE IF( ANORM.LT.ZERO ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN @@ -199,6 +199,10 @@ SUBROUTINE CGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, RETURN ELSE IF( ANORM.EQ.ZERO ) THEN RETURN + ELSE IF( SISNAN( ANORM ) ) THEN + RCOND = ANORM + INFO = 1 + RETURN END IF * SMLNUM = SLAMCH( 'Safe minimum' ) @@ -258,6 +262,11 @@ SUBROUTINE CGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, * IF( AINVNM.NE.ZERO ) $ RCOND = ( ONE / AINVNM ) / ANORM +* +* Check for NaNs +* + IF( SISNAN( RCOND ) ) + $ INFO = 1 * 20 CONTINUE RETURN diff --git a/SRC/dgecon.f b/SRC/dgecon.f index 5e6dacebf1..40361a318e 100644 --- a/SRC/dgecon.f +++ b/SRC/dgecon.f @@ -106,7 +106,7 @@ *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value -*> =-5: if ANORM is NAN or negative. +*> = 1: RCOND = NaN. *> \endverbatim * * Authors: @@ -176,7 +176,7 @@ SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 - ELSE IF( ANORM.LT.ZERO .OR. DISNAN( ANORM ) ) THEN + ELSE IF( ANORM.LT.ZERO ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN @@ -192,6 +192,10 @@ SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, RETURN ELSE IF( ANORM.EQ.ZERO ) THEN RETURN + ELSE IF( DISNAN( ANORM ) ) THEN + RCOND = ANORM + INFO = 1 + RETURN END IF * SMLNUM = DLAMCH( 'Safe minimum' ) @@ -250,6 +254,11 @@ SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, * IF( AINVNM.NE.ZERO ) $ RCOND = ( ONE / AINVNM ) / ANORM +* +* Check for NaNs +* + IF( DISNAN( RCOND ) ) + $ INFO = 1 * 20 CONTINUE RETURN diff --git a/SRC/sgecon.f b/SRC/sgecon.f index 64b96bf1f9..49c581a9b1 100644 --- a/SRC/sgecon.f +++ b/SRC/sgecon.f @@ -106,7 +106,7 @@ *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value -*> =-5: if ANORM is NAN or negative. +*> = 1: RCOND = NaN. *> \endverbatim * * Authors: @@ -176,7 +176,7 @@ SUBROUTINE SGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 - ELSE IF( ANORM.LT.ZERO .OR. SISNAN( ANORM ) ) THEN + ELSE IF( ANORM.LT.ZERO ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN @@ -192,6 +192,10 @@ SUBROUTINE SGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, RETURN ELSE IF( ANORM.EQ.ZERO ) THEN RETURN + ELSE IF( SISNAN( ANORM ) ) THEN + RCOND = ANORM + INFO = 1 + RETURN END IF * SMLNUM = SLAMCH( 'Safe minimum' ) @@ -250,6 +254,11 @@ SUBROUTINE SGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, * IF( AINVNM.NE.ZERO ) $ RCOND = ( ONE / AINVNM ) / ANORM +* +* Check for NaNs +* + IF( SISNAN( RCOND ) ) + $ INFO = 1 * 20 CONTINUE RETURN diff --git a/SRC/zgecon.f b/SRC/zgecon.f index 03fa09421f..0f8bba18e0 100644 --- a/SRC/zgecon.f +++ b/SRC/zgecon.f @@ -106,7 +106,7 @@ *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value -*> =-5: if ANORM is NAN or negative. +*> = 1: RCOND = NaN. *> \endverbatim * * Authors: @@ -183,7 +183,7 @@ SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 - ELSE IF( ANORM.LT.ZERO .OR. DISNAN( ANORM ) ) THEN + ELSE IF( ANORM.LT.ZERO ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN @@ -199,6 +199,10 @@ SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, RETURN ELSE IF( ANORM.EQ.ZERO ) THEN RETURN + ELSE IF( DISNAN( ANORM ) ) THEN + RCOND = ANORM + INFO = 1 + RETURN END IF * SMLNUM = DLAMCH( 'Safe minimum' ) @@ -258,6 +262,11 @@ SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, * IF( AINVNM.NE.ZERO ) $ RCOND = ( ONE / AINVNM ) / ANORM +* +* Check for NaNs +* + IF( DISNAN( RCOND ) ) + $ INFO = 1 * 20 CONTINUE RETURN From 68eafc69a850c5b8ee19c164de743739318b1a17 Mon Sep 17 00:00:00 2001 From: Weslley S Pereira Date: Wed, 8 Nov 2023 13:11:30 -0700 Subject: [PATCH 2/2] Makes NaN an illegal value for ANORM on input of GECON. Let NaNs propagate in the computation of RCOND --- SRC/cgecon.f | 32 ++++++++++++++++++++++++-------- SRC/dgecon.f | 32 ++++++++++++++++++++++++-------- SRC/sgecon.f | 32 ++++++++++++++++++++++++-------- SRC/zgecon.f | 32 ++++++++++++++++++++++++-------- 4 files changed, 96 insertions(+), 32 deletions(-) diff --git a/SRC/cgecon.f b/SRC/cgecon.f index ed705f796a..e018b18bb8 100644 --- a/SRC/cgecon.f +++ b/SRC/cgecon.f @@ -105,8 +105,15 @@ *> \verbatim *> INFO is INTEGER *> = 0: successful exit -*> < 0: if INFO = -i, the i-th argument had an illegal value -*> = 1: RCOND = NaN. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> NaNs are illegal values for ANORM, and they propagate to +*> the output parameter RCOND. +*> Infinity is illegal for ANORM, and it propagates to the output +*> parameter RCOND as 0. +*> = 1: if RCOND = NaN, or +*> RCOND = Inf, or +*> the computed norm of the inverse of A is 0. +*> In the latter, RCOND = 0 is returned. *> \endverbatim * * Authors: @@ -147,7 +154,7 @@ SUBROUTINE CGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, LOGICAL ONENRM CHARACTER NORMIN INTEGER IX, KASE, KASE1 - REAL AINVNM, SCALE, SL, SMLNUM, SU + REAL AINVNM, SCALE, SL, SMLNUM, SU, HUGEVAL COMPLEX ZDUM * .. * .. Local Arrays .. @@ -172,6 +179,8 @@ SUBROUTINE CGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) * .. * .. Executable Statements .. +* + HUGEVAL = SLAMCH( 'Overflow' ) * * Test the input parameters. * @@ -201,7 +210,10 @@ SUBROUTINE CGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, RETURN ELSE IF( SISNAN( ANORM ) ) THEN RCOND = ANORM - INFO = 1 + INFO = -5 + RETURN + ELSE IF( ANORM.GT.HUGEVAL ) THEN + INFO = -5 RETURN END IF * @@ -260,12 +272,16 @@ SUBROUTINE CGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, * * Compute the estimate of the reciprocal condition number. * - IF( AINVNM.NE.ZERO ) - $ RCOND = ( ONE / AINVNM ) / ANORM + IF( AINVNM.NE.ZERO ) THEN + RCOND = ( ONE / AINVNM ) / ANORM + ELSE + INFO = 1 + RETURN + END IF * -* Check for NaNs +* Check for NaNs and Infs * - IF( SISNAN( RCOND ) ) + IF( SISNAN( RCOND ) .OR. RCOND.GT.HUGEVAL ) $ INFO = 1 * 20 CONTINUE diff --git a/SRC/dgecon.f b/SRC/dgecon.f index 40361a318e..a543eb03a3 100644 --- a/SRC/dgecon.f +++ b/SRC/dgecon.f @@ -105,8 +105,15 @@ *> \verbatim *> INFO is INTEGER *> = 0: successful exit -*> < 0: if INFO = -i, the i-th argument had an illegal value -*> = 1: RCOND = NaN. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> NaNs are illegal values for ANORM, and they propagate to +*> the output parameter RCOND. +*> Infinity is illegal for ANORM, and it propagates to the output +*> parameter RCOND as 0. +*> = 1: if RCOND = NaN, or +*> RCOND = Inf, or +*> the computed norm of the inverse of A is 0. +*> In the latter, RCOND = 0 is returned. *> \endverbatim * * Authors: @@ -147,7 +154,7 @@ SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, LOGICAL ONENRM CHARACTER NORMIN INTEGER IX, KASE, KASE1 - DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU + DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU, HUGEVAL * .. * .. Local Arrays .. INTEGER ISAVE( 3 ) @@ -165,6 +172,8 @@ SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INTRINSIC ABS, MAX * .. * .. Executable Statements .. +* + HUGEVAL = DLAMCH( 'Overflow' ) * * Test the input parameters. * @@ -194,7 +203,10 @@ SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, RETURN ELSE IF( DISNAN( ANORM ) ) THEN RCOND = ANORM - INFO = 1 + INFO = -5 + RETURN + ELSE IF( ANORM.GT.HUGEVAL ) THEN + INFO = -5 RETURN END IF * @@ -252,12 +264,16 @@ SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, * * Compute the estimate of the reciprocal condition number. * - IF( AINVNM.NE.ZERO ) - $ RCOND = ( ONE / AINVNM ) / ANORM + IF( AINVNM.NE.ZERO ) THEN + RCOND = ( ONE / AINVNM ) / ANORM + ELSE + INFO = 1 + RETURN + END IF * -* Check for NaNs +* Check for NaNs and Infs * - IF( DISNAN( RCOND ) ) + IF( DISNAN( RCOND ) .OR. RCOND.GT.HUGEVAL ) $ INFO = 1 * 20 CONTINUE diff --git a/SRC/sgecon.f b/SRC/sgecon.f index 49c581a9b1..82f463ebb1 100644 --- a/SRC/sgecon.f +++ b/SRC/sgecon.f @@ -105,8 +105,15 @@ *> \verbatim *> INFO is INTEGER *> = 0: successful exit -*> < 0: if INFO = -i, the i-th argument had an illegal value -*> = 1: RCOND = NaN. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> NaNs are illegal values for ANORM, and they propagate to +*> the output parameter RCOND. +*> Infinity is illegal for ANORM, and it propagates to the output +*> parameter RCOND as 0. +*> = 1: if RCOND = NaN, or +*> RCOND = Inf, or +*> the computed norm of the inverse of A is 0. +*> In the latter, RCOND = 0 is returned. *> \endverbatim * * Authors: @@ -147,7 +154,7 @@ SUBROUTINE SGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, LOGICAL ONENRM CHARACTER NORMIN INTEGER IX, KASE, KASE1 - REAL AINVNM, SCALE, SL, SMLNUM, SU + REAL AINVNM, SCALE, SL, SMLNUM, SU, HUGEVAL * .. * .. Local Arrays .. INTEGER ISAVE( 3 ) @@ -165,6 +172,8 @@ SUBROUTINE SGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INTRINSIC ABS, MAX * .. * .. Executable Statements .. +* + HUGEVAL = SLAMCH( 'Overflow' ) * * Test the input parameters. * @@ -194,7 +203,10 @@ SUBROUTINE SGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, RETURN ELSE IF( SISNAN( ANORM ) ) THEN RCOND = ANORM - INFO = 1 + INFO = -5 + RETURN + ELSE IF( ANORM.GT.HUGEVAL ) THEN + INFO = -5 RETURN END IF * @@ -252,12 +264,16 @@ SUBROUTINE SGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, * * Compute the estimate of the reciprocal condition number. * - IF( AINVNM.NE.ZERO ) - $ RCOND = ( ONE / AINVNM ) / ANORM + IF( AINVNM.NE.ZERO ) THEN + RCOND = ( ONE / AINVNM ) / ANORM + ELSE + INFO = 1 + RETURN + END IF * -* Check for NaNs +* Check for NaNs and Infs * - IF( SISNAN( RCOND ) ) + IF( SISNAN( RCOND ) .OR. RCOND.GT.HUGEVAL ) $ INFO = 1 * 20 CONTINUE diff --git a/SRC/zgecon.f b/SRC/zgecon.f index 0f8bba18e0..ef567d7c2a 100644 --- a/SRC/zgecon.f +++ b/SRC/zgecon.f @@ -105,8 +105,15 @@ *> \verbatim *> INFO is INTEGER *> = 0: successful exit -*> < 0: if INFO = -i, the i-th argument had an illegal value -*> = 1: RCOND = NaN. +*> < 0: if INFO = -i, the i-th argument had an illegal value. +*> NaNs are illegal values for ANORM, and they propagate to +*> the output parameter RCOND. +*> Infinity is illegal for ANORM, and it propagates to the output +*> parameter RCOND as 0. +*> = 1: if RCOND = NaN, or +*> RCOND = Inf, or +*> the computed norm of the inverse of A is 0. +*> In the latter, RCOND = 0 is returned. *> \endverbatim * * Authors: @@ -147,7 +154,7 @@ SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, LOGICAL ONENRM CHARACTER NORMIN INTEGER IX, KASE, KASE1 - DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU + DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU, HUGEVAL COMPLEX*16 ZDUM * .. * .. Local Arrays .. @@ -172,6 +179,8 @@ SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) * .. * .. Executable Statements .. +* + HUGEVAL = DLAMCH( 'Overflow' ) * * Test the input parameters. * @@ -201,7 +210,10 @@ SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, RETURN ELSE IF( DISNAN( ANORM ) ) THEN RCOND = ANORM - INFO = 1 + INFO = -5 + RETURN + ELSE IF( ANORM.GT.HUGEVAL ) THEN + INFO = -5 RETURN END IF * @@ -260,12 +272,16 @@ SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, * * Compute the estimate of the reciprocal condition number. * - IF( AINVNM.NE.ZERO ) - $ RCOND = ( ONE / AINVNM ) / ANORM + IF( AINVNM.NE.ZERO ) THEN + RCOND = ( ONE / AINVNM ) / ANORM + ELSE + INFO = 1 + RETURN + END IF * -* Check for NaNs +* Check for NaNs and Infs * - IF( DISNAN( RCOND ) ) + IF( DISNAN( RCOND ) .OR. RCOND.GT.HUGEVAL ) $ INFO = 1 * 20 CONTINUE