**** Fix bug 124 ****

**** 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".