**** Fix bug 124 ****

**** (I also fixed the documentation of the six subroutines ****

Bug reported by Victor_K on the forum.

See forum http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4643

Bug fix by Julien Langou. Involved as well Nick Higham, Sven Hammarling and Craig Lucas.

** Bug **

xPSTF2 and xPSTRF returns NaNs when the matrix in input has N negative diagonal entries.

** Note **

A matrix with N negative diagonal entries is a highly non valid input, but, yet no reason to create some NaNs, better to have a cleaner solution.

** Description **

xPSTF2 and xPSTRF first scan the diagonal of the matrix to find the greatest element on the diagonal.

During this pass on the diagonal, xPSTF2/xPSTRF checks as well for ZERO or NAN on the diagonal and aborts if it finds a ZERO or a NAN on the diagonal.

*
*     Compute stopping value
*
      PVT = 1
      AJJ = A( PVT, PVT )
      DO I = 2, N
         IF( A( I, I ).GT.AJJ ) THEN
            PVT = I
            AJJ = A( PVT, PVT )
         END IF
      END DO
      IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
         RANK = 0
         INFO = 1
         GO TO 170
      END IF

Since the input matrix is filled with a negative diagonal (and no ZERO and no NAN), xPSTF2/xPSTRF finds a maximum element (which is the largest negative elements on the diagonal), and then this (negative) element is used to compute DSTOP, the stopping value. But then we get a negative DSTOP,

      IF( TOL.LT.ZERO ) THEN
         DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
      ELSE
         DSTOP = TOL
      END IF

and then I stop even thinking on what would happen next. It seems that the code continues and does produce NAN in the first column. Why not. But I think we should catch this case and try to handle it in a different way. Maybe.

** Bug fix **

Replace
      IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
With
      IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
So that if there is any nonpositive element or any NAN on the diagonal of A, we return right away with RANK=0 and INFO=1.

** Email from Nick Higham **

[ Julien wrote ] So that if there is any nonpositive element or any NAN on the diagonal of A, we return right away with RANK=0 and INFO=1?

[ Nick wrote ]

I think so.
This is justified by the documentation, which says

   INFO is INTEGER
   If INFO = -K, the K-th argument had an illegal value,
   = 0: algorithm completed successfully, and
   > 0: the matrix A is either rank deficient with computed rank
   as returned in RANK, or is indefinite.

We are in the indefinite case if any diagonal element is negative.
Actually, that is not quite true: we are in the "not positive definite"
case.  We can only say "indefinite" if there are both positive *and* negative
diagonal entries!  However, if that is going to be said then it would be
necessary to check that all other LAPACK routines are similarly precise in
their wording, so it's simplest, and not too misleading,
to stay with "indefinite".

diff --git a/SRC/cpstf2.f b/SRC/cpstf2.f
index 246c996..623eb62 100644 (file)
--- a/SRC/cpstf2.f
+++ b/SRC/cpstf2.f
@@ -1,4 +1,4 @@
-*> \brief \b CPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric or complex Hermitian positive semi-definite matrix.
+*> \brief \b CPSTF2 computes the Cholesky factorization with complete pivoting of complex Hermitian positive semidefinite matrix.
 *
 *  =========== DOCUMENTATION ===========
 *
@@ -122,8 +122,9 @@
 *>          < 0: If INFO = -K, the K-th argument had an illegal value,
 *>          = 0: algorithm completed successfully, and
 *>          > 0: the matrix A is either rank deficient with computed rank
-*>               as returned in RANK, or is indefinite.  See Section 7 of
-*>               LAPACK Working Note #161 for further information.
+*>               as returned in RANK, or is not positive semidefinite. See
+*>               Section 7 of LAPACK Working Note #161 for further
+*>               information.
 *> \endverbatim
 *
 *  Authors:
@@ -218,7 +219,7 @@
   110 CONTINUE
       PVT = MAXLOC( WORK( 1:N ), 1 )
       AJJ = REAL ( A( PVT, PVT ) )
-      IF( AJJ.EQ.ZERO.OR.SISNAN( AJJ ) ) THEN
+      IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
          RANK = 0
          INFO = 1
          GO TO 200

diff --git a/SRC/cpstrf.f b/SRC/cpstrf.f
index 2afbc57..e0c651a 100644 (file)
--- a/SRC/cpstrf.f
+++ b/SRC/cpstrf.f
@@ -1,4 +1,4 @@
-*> \brief \b CPSTRF
+*> \brief \b CPSTRF computes the Cholesky factorization with complete pivoting of complex Hermitian positive semidefinite matrix.
 *
 *  =========== DOCUMENTATION ===========
 *
@@ -122,8 +122,9 @@
 *>          < 0: If INFO = -K, the K-th argument had an illegal value,
 *>          = 0: algorithm completed successfully, and
 *>          > 0: the matrix A is either rank deficient with computed rank
-*>               as returned in RANK, or is indefinite.  See Section 7 of
-*>               LAPACK Working Note #161 for further information.
+*>               as returned in RANK, or is not positive semidefinite. See
+*>               Section 7 of LAPACK Working Note #161 for further
+*>               information.
 *> \endverbatim
 *
 *  Authors:
@@ -233,7 +234,7 @@
   110    CONTINUE
          PVT = MAXLOC( WORK( 1:N ), 1 )
          AJJ = REAL( A( PVT, PVT ) )
-         IF( AJJ.EQ.ZERO.OR.SISNAN( AJJ ) ) THEN
+         IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
             RANK = 0
             INFO = 1
             GO TO 230

diff --git a/SRC/dpstf2.f b/SRC/dpstf2.f
index a055853..d2cedcd 100644 (file)
--- a/SRC/dpstf2.f
+++ b/SRC/dpstf2.f
@@ -1,4 +1,4 @@
-*> \brief \b DPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric or complex Hermitian positive semi-definite matrix.
+*> \brief \b DPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric positive semidefinite matrix.
 *
 *  =========== DOCUMENTATION ===========
 *
@@ -121,8 +121,9 @@
 *>          < 0: If INFO = -K, the K-th argument had an illegal value,
 *>          = 0: algorithm completed successfully, and
 *>          > 0: the matrix A is either rank deficient with computed rank
-*>               as returned in RANK, or is indefinite.  See Section 7 of
-*>               LAPACK Working Note #161 for further information.
+*>               as returned in RANK, or is not positive semidefinite. See
+*>               Section 7 of LAPACK Working Note #161 for further
+*>               information.
 *> \endverbatim
 *
 *  Authors:
@@ -216,7 +217,7 @@
             AJJ = A( PVT, PVT )
          END IF
       END DO
-      IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
+      IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
          RANK = 0
          INFO = 1
          GO TO 170

diff --git a/SRC/dpstrf.f b/SRC/dpstrf.f
index 2f64aad..066ea85 100644 (file)
--- a/SRC/dpstrf.f
+++ b/SRC/dpstrf.f
@@ -1,4 +1,5 @@
-*> \brief \b DPSTRF
+*> \brief \b DPSTRF computes the Cholesky factorization with complete pivoting of a real symmetric positive semidefinite matrix.
+*
 *
 *  =========== DOCUMENTATION ===========
 *
@@ -121,8 +122,9 @@
 *>          < 0: If INFO = -K, the K-th argument had an illegal value,
 *>          = 0: algorithm completed successfully, and
 *>          > 0: the matrix A is either rank deficient with computed rank
-*>               as returned in RANK, or is indefinite.  See Section 7 of
-*>               LAPACK Working Note #161 for further information.
+*>               as returned in RANK, or is not positive semidefinite. See
+*>               Section 7 of LAPACK Working Note #161 for further
+*>               information.
 *> \endverbatim
 *
 *  Authors:
@@ -230,7 +232,7 @@
                AJJ = A( PVT, PVT )
             END IF
          END DO
-         IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
+         IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
             RANK = 0
             INFO = 1
             GO TO 200

diff --git a/SRC/spstf2.f b/SRC/spstf2.f
index 67dbf7d..3214979 100644 (file)
--- a/SRC/spstf2.f
+++ b/SRC/spstf2.f
@@ -1,4 +1,4 @@
-*> \brief \b SPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric or complex Hermitian positive semi-definite matrix.
+*> \brief \b SPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric positive semidefinite matrix.
 *
 *  =========== DOCUMENTATION ===========
 *
@@ -121,8 +121,9 @@
 *>          < 0: If INFO = -K, the K-th argument had an illegal value,
 *>          = 0: algorithm completed successfully, and
 *>          > 0: the matrix A is either rank deficient with computed rank
-*>               as returned in RANK, or is indefinite.  See Section 7 of
-*>               LAPACK Working Note #161 for further information.
+*>               as returned in RANK, or is not positive semidefinite. See
+*>               Section 7 of LAPACK Working Note #161 for further
+*>               information.
 *> \endverbatim
 *
 *  Authors:
@@ -216,7 +217,7 @@
             AJJ = A( PVT, PVT )
          END IF
       END DO
-      IF( AJJ.EQ.ZERO.OR.SISNAN( AJJ ) ) THEN
+      IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
          RANK = 0
          INFO = 1
          GO TO 170

diff --git a/SRC/spstrf.f b/SRC/spstrf.f
index 76fb870..8049a76 100644 (file)
--- a/SRC/spstrf.f
+++ b/SRC/spstrf.f
@@ -1,4 +1,4 @@
-*> \brief \b SPSTRF
+*> \brief \b SPSTRF computes the Cholesky factorization with complete pivoting of a real symmetric positive semidefinite matrix.
 *
 *  =========== DOCUMENTATION ===========
 *
@@ -121,8 +121,9 @@
 *>          < 0: If INFO = -K, the K-th argument had an illegal value,
 *>          = 0: algorithm completed successfully, and
 *>          > 0: the matrix A is either rank deficient with computed rank
-*>               as returned in RANK, or is indefinite.  See Section 7 of
-*>               LAPACK Working Note #161 for further information.
+*>               as returned in RANK, or is not positive semidefinite. See
+*>               Section 7 of LAPACK Working Note #161 for further
+*>               information.
 *> \endverbatim
 *
 *  Authors:
@@ -230,7 +231,7 @@
                AJJ = A( PVT, PVT )
             END IF
          END DO
-         IF( AJJ.EQ.ZERO.OR.SISNAN( AJJ ) ) THEN
+         IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
             RANK = 0
             INFO = 1
             GO TO 200

diff --git a/SRC/zpstf2.f b/SRC/zpstf2.f
index 4e8a8d0..bb40c45 100644 (file)
--- a/SRC/zpstf2.f
+++ b/SRC/zpstf2.f
@@ -1,4 +1,4 @@
-*> \brief \b ZPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric or complex Hermitian positive semi-definite matrix.
+*> \brief \b ZPSTF2 computes the Cholesky factorization with complete pivoting of a complex Hermitian positive semidefinite matrix.
 *
 *  =========== DOCUMENTATION ===========
 *
@@ -122,8 +122,9 @@
 *>          < 0: If INFO = -K, the K-th argument had an illegal value,
 *>          = 0: algorithm completed successfully, and
 *>          > 0: the matrix A is either rank deficient with computed rank
-*>               as returned in RANK, or is indefinite.  See Section 7 of
-*>               LAPACK Working Note #161 for further information.
+*>               as returned in RANK, or is not positive semidefinite. See
+*>               Section 7 of LAPACK Working Note #161 for further
+*>               information.
 *> \endverbatim
 *
 *  Authors:
@@ -218,7 +219,7 @@
   110 CONTINUE
       PVT = MAXLOC( WORK( 1:N ), 1 )
       AJJ = DBLE( A( PVT, PVT ) )
-      IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
+      IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
          RANK = 0
          INFO = 1
          GO TO 200

diff --git a/SRC/zpstrf.f b/SRC/zpstrf.f
index dd00700..433c826 100644 (file)
--- a/SRC/zpstrf.f
+++ b/SRC/zpstrf.f
@@ -1,4 +1,4 @@
-*> \brief \b ZPSTRF
+*> \brief \b ZPSTRF computes the Cholesky factorization with complete pivoting of a complex Hermitian positive semidefinite matrix.
 *
 *  =========== DOCUMENTATION ===========
 *
@@ -122,8 +122,9 @@
 *>          < 0: If INFO = -K, the K-th argument had an illegal value,
 *>          = 0: algorithm completed successfully, and
 *>          > 0: the matrix A is either rank deficient with computed rank
-*>               as returned in RANK, or is indefinite.  See Section 7 of
-*>               LAPACK Working Note #161 for further information.
+*>               as returned in RANK, or is not positive semidefinite. See
+*>               Section 7 of LAPACK Working Note #161 for further
+*>               information.
 *> \endverbatim
 *
 *  Authors:
@@ -233,7 +234,7 @@
   110    CONTINUE
          PVT = MAXLOC( WORK( 1:N ), 1 )
          AJJ = DBLE( A( PVT, PVT ) )
-         IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
+         IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
             RANK = 0
             INFO = 1
             GO TO 230