From 8bef08dd7d3c0c1ca7e150a7cbbb768b98d9152b Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Mon, 7 Sep 2026 11:53:21 +0200 Subject: [PATCH] Guard against a NaN pivot (Reference-LAPACK PR 1378) --- lapack-netlib/SRC/chptrf.f | 27 ++++++++++++++++----------- lapack-netlib/SRC/csptrf.f | 24 ++++++++++++++---------- lapack-netlib/SRC/dsptrf.f | 24 ++++++++++++++---------- lapack-netlib/SRC/ssptrf.f | 24 ++++++++++++++---------- lapack-netlib/SRC/zhptrf.f | 27 ++++++++++++++++----------- lapack-netlib/SRC/zsptrf.f | 24 ++++++++++++++---------- 6 files changed, 88 insertions(+), 62 deletions(-) diff --git a/lapack-netlib/SRC/chptrf.f b/lapack-netlib/SRC/chptrf.f index d9bfffb06e..3d584fd513 100644 --- a/lapack-netlib/SRC/chptrf.f +++ b/lapack-netlib/SRC/chptrf.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download CHPTRF + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -107,7 +105,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complexOTHERcomputational +*> \ingroup hptrf * *> \par Further Details: * ===================== @@ -156,6 +154,7 @@ *> * ===================================================================== SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -187,10 +186,10 @@ SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO ) COMPLEX D12, D21, T, WK, WKM1, WKP1, ZDUM * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, SISNAN INTEGER ICAMAX REAL SLAPY2 - EXTERNAL LSAME, ICAMAX, SLAPY2 + EXTERNAL LSAME, ICAMAX, SLAPY2, SISNAN * .. * .. External Subroutines .. EXTERNAL CHPR, CSSCAL, CSWAP, XERBLA @@ -257,9 +256,11 @@ SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ SISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -460,9 +461,11 @@ SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ SISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -524,7 +527,8 @@ SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO ) * submatrix A(k:n,k:n) * IF( KP.LT.N ) - $ CALL CSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ), + $ CALL CSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, + $ AP( KPC+1 ), $ 1 ) KX = KNC + KP - KK DO 80 J = KK + 1, KP - 1 @@ -592,7 +596,8 @@ SUBROUTINE CHPTRF( UPLO, N, AP, IPIV, INFO ) * where L(k) and L(k+1) are the k-th and (k+1)-th * columns of L * - D = SLAPY2( REAL( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ), + D = SLAPY2( + $ REAL( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ), $ AIMAG( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ) ) D11 = REAL( AP( K+1+K*( 2*N-K-1 ) / 2 ) ) / D D22 = REAL( AP( K+( K-1 )*( 2*N-K ) / 2 ) ) / D diff --git a/lapack-netlib/SRC/csptrf.f b/lapack-netlib/SRC/csptrf.f index 702ed0c518..64f8e678b6 100644 --- a/lapack-netlib/SRC/csptrf.f +++ b/lapack-netlib/SRC/csptrf.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download CSPTRF + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -108,7 +106,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complexOTHERcomputational +*> \ingroup hptrf * *> \par Further Details: * ===================== @@ -155,6 +153,7 @@ *> * ===================================================================== SUBROUTINE CSPTRF( UPLO, N, AP, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -187,9 +186,9 @@ SUBROUTINE CSPTRF( UPLO, N, AP, IPIV, INFO ) COMPLEX D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, ZDUM * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, SISNAN INTEGER ICAMAX - EXTERNAL LSAME, ICAMAX + EXTERNAL LSAME, ICAMAX, SISNAN * .. * .. External Subroutines .. EXTERNAL CSCAL, CSPR, CSWAP, XERBLA @@ -256,9 +255,11 @@ SUBROUTINE CSPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ SISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -443,9 +444,11 @@ SUBROUTINE CSPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ SISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -506,7 +509,8 @@ SUBROUTINE CSPTRF( UPLO, N, AP, IPIV, INFO ) * submatrix A(k:n,k:n) * IF( KP.LT.N ) - $ CALL CSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ), + $ CALL CSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, + $ AP( KPC+1 ), $ 1 ) KX = KNC + KP - KK DO 80 J = KK + 1, KP - 1 diff --git a/lapack-netlib/SRC/dsptrf.f b/lapack-netlib/SRC/dsptrf.f index c81d8bcbf8..acd0290435 100644 --- a/lapack-netlib/SRC/dsptrf.f +++ b/lapack-netlib/SRC/dsptrf.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download DSPTRF + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -107,7 +105,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup doubleOTHERcomputational +*> \ingroup hptrf * *> \par Further Details: * ===================== @@ -156,6 +154,7 @@ *> * ===================================================================== SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -186,9 +185,9 @@ SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO ) $ ROWMAX, T, WK, WKM1, WKP1 * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, DISNAN INTEGER IDAMAX - EXTERNAL LSAME, IDAMAX + EXTERNAL LSAME, IDAMAX, DISNAN * .. * .. External Subroutines .. EXTERNAL DSCAL, DSPR, DSWAP, XERBLA @@ -249,9 +248,11 @@ SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ DISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -437,9 +438,11 @@ SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ DISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -500,7 +503,8 @@ SUBROUTINE DSPTRF( UPLO, N, AP, IPIV, INFO ) * submatrix A(k:n,k:n) * IF( KP.LT.N ) - $ CALL DSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ), + $ CALL DSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, + $ AP( KPC+1 ), $ 1 ) KX = KNC + KP - KK DO 80 J = KK + 1, KP - 1 diff --git a/lapack-netlib/SRC/ssptrf.f b/lapack-netlib/SRC/ssptrf.f index 4a25e3d073..3c4456a148 100644 --- a/lapack-netlib/SRC/ssptrf.f +++ b/lapack-netlib/SRC/ssptrf.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download SSPTRF + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -107,7 +105,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup realOTHERcomputational +*> \ingroup hptrf * *> \par Further Details: * ===================== @@ -154,6 +152,7 @@ *> * ===================================================================== SUBROUTINE SSPTRF( UPLO, N, AP, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -184,9 +183,9 @@ SUBROUTINE SSPTRF( UPLO, N, AP, IPIV, INFO ) $ ROWMAX, T, WK, WKM1, WKP1 * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, SISNAN INTEGER ISAMAX - EXTERNAL LSAME, ISAMAX + EXTERNAL LSAME, ISAMAX, SISNAN * .. * .. External Subroutines .. EXTERNAL SSCAL, SSPR, SSWAP, XERBLA @@ -247,9 +246,11 @@ SUBROUTINE SSPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ SISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -435,9 +436,11 @@ SUBROUTINE SSPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ SISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -498,7 +501,8 @@ SUBROUTINE SSPTRF( UPLO, N, AP, IPIV, INFO ) * submatrix A(k:n,k:n) * IF( KP.LT.N ) - $ CALL SSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ), + $ CALL SSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, + $ AP( KPC+1 ), $ 1 ) KX = KNC + KP - KK DO 80 J = KK + 1, KP - 1 diff --git a/lapack-netlib/SRC/zhptrf.f b/lapack-netlib/SRC/zhptrf.f index 9c736f3352..6558a1bb80 100644 --- a/lapack-netlib/SRC/zhptrf.f +++ b/lapack-netlib/SRC/zhptrf.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download ZHPTRF + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -107,7 +105,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16OTHERcomputational +*> \ingroup hptrf * *> \par Further Details: * ===================== @@ -156,6 +154,7 @@ * * ===================================================================== SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -187,10 +186,10 @@ SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO ) COMPLEX*16 D12, D21, T, WK, WKM1, WKP1, ZDUM * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, DISNAN INTEGER IZAMAX DOUBLE PRECISION DLAPY2 - EXTERNAL LSAME, IZAMAX, DLAPY2 + EXTERNAL LSAME, IZAMAX, DLAPY2, DISNAN * .. * .. External Subroutines .. EXTERNAL XERBLA, ZDSCAL, ZHPR, ZSWAP @@ -257,9 +256,11 @@ SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ DISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -460,9 +461,11 @@ SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ DISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -524,7 +527,8 @@ SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO ) * submatrix A(k:n,k:n) * IF( KP.LT.N ) - $ CALL ZSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ), + $ CALL ZSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, + $ AP( KPC+1 ), $ 1 ) KX = KNC + KP - KK DO 80 J = KK + 1, KP - 1 @@ -592,7 +596,8 @@ SUBROUTINE ZHPTRF( UPLO, N, AP, IPIV, INFO ) * where L(k) and L(k+1) are the k-th and (k+1)-th * columns of L * - D = DLAPY2( DBLE( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ), + D = DLAPY2( + $ DBLE( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ), $ DIMAG( AP( K+1+( K-1 )*( 2*N-K ) / 2 ) ) ) D11 = DBLE( AP( K+1+K*( 2*N-K-1 ) / 2 ) ) / D D22 = DBLE( AP( K+( K-1 )*( 2*N-K ) / 2 ) ) / D diff --git a/lapack-netlib/SRC/zsptrf.f b/lapack-netlib/SRC/zsptrf.f index bad30b7a42..4a5a7a50b8 100644 --- a/lapack-netlib/SRC/zsptrf.f +++ b/lapack-netlib/SRC/zsptrf.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download ZSPTRF + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -108,7 +106,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16OTHERcomputational +*> \ingroup hptrf * *> \par Further Details: * ===================== @@ -155,6 +153,7 @@ *> * ===================================================================== SUBROUTINE ZSPTRF( UPLO, N, AP, IPIV, INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -187,9 +186,9 @@ SUBROUTINE ZSPTRF( UPLO, N, AP, IPIV, INFO ) COMPLEX*16 D11, D12, D21, D22, R1, T, WK, WKM1, WKP1, ZDUM * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL LSAME, DISNAN INTEGER IZAMAX - EXTERNAL LSAME, IZAMAX + EXTERNAL LSAME, IZAMAX, DISNAN * .. * .. External Subroutines .. EXTERNAL XERBLA, ZSCAL, ZSPR, ZSWAP @@ -256,9 +255,11 @@ SUBROUTINE ZSPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ DISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -443,9 +444,11 @@ SUBROUTINE ZSPTRF( UPLO, N, AP, IPIV, INFO ) COLMAX = ZERO END IF * - IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN + IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. + $ DISNAN(ABSAKK) ) THEN * -* Column K is zero: set INFO and continue +* Column K is zero or underflow, or contains a NaN: +* set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K @@ -506,7 +509,8 @@ SUBROUTINE ZSPTRF( UPLO, N, AP, IPIV, INFO ) * submatrix A(k:n,k:n) * IF( KP.LT.N ) - $ CALL ZSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, AP( KPC+1 ), + $ CALL ZSWAP( N-KP, AP( KNC+KP-KK+1 ), 1, + $ AP( KPC+1 ), $ 1 ) KX = KNC + KP - KK DO 80 J = KK + 1, KP - 1