From 5dbf71b11b4901dd61b5d820d1c937a2ec4d0fcf Mon Sep 17 00:00:00 2001 From: james Date: Thu, 11 Oct 2012 18:19:31 +0000 Subject: [PATCH] * added check for NaN after the norm of each block of H is computed * if NaN is detected, returns with INFO=-6 (H is the 6th parameter) * fixes problem found by Alexander Kobotov at Intel where a NaN in H can cause an infinite loop: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=3928 --- SRC/chsein.f | 10 +++++++--- SRC/dhsein.f | 10 +++++++--- SRC/shsein.f | 10 +++++++--- SRC/zhsein.f | 10 +++++++--- 4 files changed, 28 insertions(+), 12 deletions(-) diff --git a/SRC/chsein.f b/SRC/chsein.f index fa87a73..a97b700 100644 --- a/SRC/chsein.f +++ b/SRC/chsein.f @@ -104,6 +104,7 @@ *> \verbatim *> H is COMPLEX array, dimension (LDH,N) *> The upper Hessenberg matrix H. +*> If a NaN is detected in H, the routine will return with INFO=-6. *> \endverbatim *> *> \param[in] LDH @@ -276,9 +277,9 @@ COMPLEX CDUM, WK * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, SISNAN REAL CLANHS, SLAMCH - EXTERNAL LSAME, CLANHS, SLAMCH + EXTERNAL LSAME, CLANHS, SLAMCH, SISNAN * .. * .. External Subroutines .. EXTERNAL CLAEIN, XERBLA @@ -399,7 +400,10 @@ * has not ben computed before. * HNORM = CLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK ) - IF( HNORM.GT.RZERO ) THEN + IF( SISNAN( HNORM ) ) THEN + INFO = -6 + RETURN + ELSE IF( (HNORM.GT.RZERO) ) THEN EPS3 = HNORM*ULP ELSE EPS3 = SMLNUM diff --git a/SRC/dhsein.f b/SRC/dhsein.f index d239a0e..5b1bad6 100644 --- a/SRC/dhsein.f +++ b/SRC/dhsein.f @@ -108,6 +108,7 @@ *> \verbatim *> H is DOUBLE PRECISION array, dimension (LDH,N) *> The upper Hessenberg matrix H. +*> If a NaN is detected in H, the routine will return with INFO=-6. *> \endverbatim *> *> \param[in] LDH @@ -291,9 +292,9 @@ $ WKR * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, DISNAN DOUBLE PRECISION DLAMCH, DLANHS - EXTERNAL LSAME, DLAMCH, DLANHS + EXTERNAL LSAME, DLAMCH, DLANHS, DISNAN * .. * .. External Subroutines .. EXTERNAL DLAEIN, XERBLA @@ -423,7 +424,10 @@ * has not ben computed before. * HNORM = DLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, WORK ) - IF( HNORM.GT.ZERO ) THEN + IF( DISNAN( HNORM ) ) THEN + INFO = -6 + RETURN + ELSE IF( HNORM.GT.ZERO ) THEN EPS3 = HNORM*ULP ELSE EPS3 = SMLNUM diff --git a/SRC/shsein.f b/SRC/shsein.f index 0c69e2c..322c5e7 100644 --- a/SRC/shsein.f +++ b/SRC/shsein.f @@ -108,6 +108,7 @@ *> \verbatim *> H is REAL array, dimension (LDH,N) *> The upper Hessenberg matrix H. +*> If a NaN is detected in H, the routine will return with INFO=-6. *> \endverbatim *> *> \param[in] LDH @@ -291,9 +292,9 @@ $ WKR * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, SISNAN REAL SLAMCH, SLANHS - EXTERNAL LSAME, SLAMCH, SLANHS + EXTERNAL LSAME, SLAMCH, SLANHS, SISNAN * .. * .. External Subroutines .. EXTERNAL SLAEIN, XERBLA @@ -423,7 +424,10 @@ * has not ben computed before. * HNORM = SLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, WORK ) - IF( HNORM.GT.ZERO ) THEN + IF( SISNAN( HNORM ) ) THEN + INFO = -6 + RETURN + ELSE IF( HNORM.GT.ZERO ) THEN EPS3 = HNORM*ULP ELSE EPS3 = SMLNUM diff --git a/SRC/zhsein.f b/SRC/zhsein.f index 6c7173f..059f366 100644 --- a/SRC/zhsein.f +++ b/SRC/zhsein.f @@ -104,6 +104,7 @@ *> \verbatim *> H is COMPLEX*16 array, dimension (LDH,N) *> The upper Hessenberg matrix H. +*> If a NaN is detected in H, the routine will return with INFO=-6. *> \endverbatim *> *> \param[in] LDH @@ -276,9 +277,9 @@ COMPLEX*16 CDUM, WK * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, DISNAN DOUBLE PRECISION DLAMCH, ZLANHS - EXTERNAL LSAME, DLAMCH, ZLANHS + EXTERNAL LSAME, DLAMCH, ZLANHS, DISNAN * .. * .. External Subroutines .. EXTERNAL XERBLA, ZLAEIN @@ -399,7 +400,10 @@ * has not ben computed before. * HNORM = ZLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK ) - IF( HNORM.GT.RZERO ) THEN + IF( DISNAN( HNORM ) ) THEN + INFO = -6 + RETURN + ELSE IF( HNORM.GT.RZERO ) THEN EPS3 = HNORM*ULP ELSE EPS3 = SMLNUM -- 2.7.4