Fix to Ming Gu's bug in dqds. See post of LAPACK DEV message board for more details

8 files changed, 152 insertions, 19 deletions

diff --git a/SRC/cbdsqr.f b/SRC/cbdsqr.f
index 6142a818..12cbb748 100644
--- a/SRC/cbdsqr.f
+++ b/SRC/cbdsqr.f
@@ -106,7 +106,8 @@
 *          The leading dimension of the array C.
 *          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
 *
-*  RWORK   (workspace) REAL array, dimension (4*N) 
+*  RWORK   (workspace) REAL array, dimension (2*N) 
+*          if NCVT = NRU = NCC = 0, (max(1, 4*N-4)) otherwise
 *
 *  INFO    (output) INTEGER
 *          = 0:  successful exit
@@ -223,7 +224,11 @@
 *
       IF( .NOT.ROTATE ) THEN
          CALL SLASQ1( N, D, E, RWORK, INFO )
-         RETURN
+*
+*     If INFO equals 2, dqds didn't finish, try to finish
+*         
+         IF( INFO .NE. 2 ) RETURN
+         INFO = 0
       END IF
 *
       NM1 = N - 1
diff --git a/SRC/dbdsqr.f b/SRC/dbdsqr.f
index 4881a69f..0c33294a 100644
--- a/SRC/dbdsqr.f
+++ b/SRC/dbdsqr.f
@@ -231,7 +231,11 @@
 *
       IF( .NOT.ROTATE ) THEN
          CALL DLASQ1( N, D, E, WORK, INFO )
-         RETURN
+*
+*     If INFO equals 2, dqds didn't finish, try to finish
+*         
+         IF( INFO .NE. 2 ) RETURN
+         INFO = 0
       END IF
 *
       NM1 = N - 1
diff --git a/SRC/dlasq1.f b/SRC/dlasq1.f
index 53fd3ddd..9ecb07b6 100644
--- a/SRC/dlasq1.f
+++ b/SRC/dlasq1.f
@@ -56,8 +56,11 @@
 *        < 0: if INFO = -i, the i-th argument had an illegal value
 *        > 0: the algorithm failed
 *             = 1, a split was marked by a positive value in E
-*             = 2, current block of Z not diagonalized after 30*N
-*                  iterations (in inner while loop)
+*             = 2, current block of Z not diagonalized after 100*N
+*                  iterations (in inner while loop)  On exit D and E
+*                  represent a matrix with the same singular values
+*                  which the calling subroutine could use to finish the
+*                  computation, or even feed back into DLASQ1
 *             = 3, termination criterion of outer while loop not met 
 *                  (program created more than N unreduced blocks)
 *
@@ -145,6 +148,17 @@
             D( I ) = SQRT( WORK( I ) )
    40    CONTINUE
          CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
+      ELSE IF( INFO.EQ.2 ) THEN
+*
+*     Maximum number of iterations exceeded.  Move data from WORK
+*     into D and E so the calling subroutine can try to finish
+*
+         DO I = 1, N
+            D( I ) = SQRT( WORK( 2*I-1 ) )
+            E( I ) = SQRT( WORK( 2*I ) )
+         END DO
+         CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
+         CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, E, N, IINFO )
       END IF
 *
       RETURN
diff --git a/SRC/dlasq2.f b/SRC/dlasq2.f
index 00bd2b14..82639956 100644
--- a/SRC/dlasq2.f
+++ b/SRC/dlasq2.f
@@ -58,8 +58,9 @@
 *             INFO = -(i*100+j)
 *        > 0: the algorithm failed
 *              = 1, a split was marked by a positive value in E
-*              = 2, current block of Z not diagonalized after 30*N
-*                   iterations (in inner while loop)
+*              = 2, current block of Z not diagonalized after 100*N
+*                   iterations (in inner while loop).  On exit Z holds
+*                   a qd array with the same eigenvalues as the given Z.
 *              = 3, termination criterion of outer while loop not met 
 *                   (program created more than N unreduced blocks)
 *
@@ -86,7 +87,7 @@
       DOUBLE PRECISION   D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN,
      $                   DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX,
      $                   QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL,
-     $                   TOL2, TRACE, ZMAX
+     $                   TOL2, TRACE, ZMAX, TEMPE, TEMPQ
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           DLASQ3, DLASRT, XERBLA
@@ -397,7 +398,7 @@
 *               and that the tests for deflation upon entry in DLASQ3 
 *               should not be performed.
 *
-         NBIG = 30*( N0-I0+1 )
+         NBIG = 100*( N0-I0+1 )
          DO 140 IWHILB = 1, NBIG
             IF( I0.GT.N0 ) 
      $         GO TO 150
@@ -442,6 +443,47 @@
   140    CONTINUE
 *
          INFO = 2
+*       
+*        Maximum number of iterations exceeded, restore the shift 
+*        SIGMA and place the new d's and e's in a qd array.
+*        This might need to be done for several blocks
+*
+         I1 = I0
+         N1 = N0
+ 145     CONTINUE
+         TEMPQ = Z( 4*I0-3 )
+         Z( 4*I0-3 ) = Z( 4*I0-3 ) + SIGMA
+         DO K = I0+1, N0
+            TEMPE = Z( 4*K-5 )
+            Z( 4*K-5 ) = Z( 4*K-5 ) * (TEMPQ / Z( 4*K-7 ))
+            TEMPQ = Z( 4*K-3 )
+            Z( 4*K-3 ) = Z( 4*K-3 ) + SIGMA + TEMPE - Z( 4*K-5 )
+         END DO
+*
+*        Prepare to do this on the previous block if there is one
+*
+         IF( I1.GT.1 ) THEN
+            N1 = I1-1
+            DO WHILE( ( I1.GE.2 ) .AND. ( Z(4*I1-5).GE.ZERO ) )
+               I1 = I1 - 1
+            END DO
+            SIGMA = -Z(4*N1-1)
+            GO TO 145
+         END IF
+
+         DO K = 1, N
+            Z( 2*K-1 ) = Z( 4*K-3 )
+*
+*        Only the block 1..N0 is unfinished.  The rest of the e's
+*        must be essentially zero, although sometimes other data
+*        has been stored in them.
+*
+            IF( K.LT.N0 ) THEN
+               Z( 2*K ) = Z( 4*K-1 )
+            ELSE
+               Z( 2*K ) = 0
+            END IF
+         END DO
          RETURN
 *
 *        end IWHILB
diff --git a/SRC/sbdsqr.f b/SRC/sbdsqr.f
index e3148433..eb3b6a89 100644
--- a/SRC/sbdsqr.f
+++ b/SRC/sbdsqr.f
@@ -231,7 +231,11 @@
 *
       IF( .NOT.ROTATE ) THEN
          CALL SLASQ1( N, D, E, WORK, INFO )
-         RETURN
+*
+*     If INFO equals 2, dqds didn't finish, try to finish
+*         
+         IF( INFO .NE. 2 ) RETURN
+         INFO = 0
       END IF
 *
       NM1 = N - 1
diff --git a/SRC/slasq1.f b/SRC/slasq1.f
index 5a063a7c..56d85d02 100644
--- a/SRC/slasq1.f
+++ b/SRC/slasq1.f
@@ -56,8 +56,11 @@
 *        < 0: if INFO = -i, the i-th argument had an illegal value
 *        > 0: the algorithm failed
 *             = 1, a split was marked by a positive value in E
-*             = 2, current block of Z not diagonalized after 30*N
-*                  iterations (in inner while loop)
+*             = 2, current block of Z not diagonalized after 100*N
+*                  iterations (in inner while loop)  On exit D and E
+*                  represent a matrix with the same singular values
+*                  which the calling subroutine could use to finish the
+*                  computation, or even feed back into SLASQ1
 *             = 3, termination criterion of outer while loop not met 
 *                  (program created more than N unreduced blocks)
 *
@@ -145,6 +148,17 @@
             D( I ) = SQRT( WORK( I ) )
    40    CONTINUE
          CALL SLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
+      ELSE IF( INFO.EQ.2 ) THEN
+*
+*     Maximum number of iterations exceeded.  Move data from WORK
+*     into D and E so the calling subroutine can try to finish
+*
+         DO I = 1, N
+            D( I ) = SQRT( WORK( 2*I-1 ) )
+            E( I ) = SQRT( WORK( 2*I ) )
+         END DO
+         CALL SLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
+         CALL SLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, E, N, IINFO )
       END IF
 *
       RETURN
diff --git a/SRC/slasq2.f b/SRC/slasq2.f
index 05c3e7cb..47af6149 100644
--- a/SRC/slasq2.f
+++ b/SRC/slasq2.f
@@ -58,8 +58,9 @@
 *             INFO = -(i*100+j)
 *        > 0: the algorithm failed
 *              = 1, a split was marked by a positive value in E
-*              = 2, current block of Z not diagonalized after 30*N
-*                   iterations (in inner while loop)
+*              = 2, current block of Z not diagonalized after 100*N
+*                   iterations (in inner while loop).  On exit Z holds
+*                   a qd array with the same eigenvalues as the given Z.
 *              = 3, termination criterion of outer while loop not met 
 *                   (program created more than N unreduced blocks)
 *
@@ -81,11 +82,12 @@
 *     .. Local Scalars ..
       LOGICAL            IEEE
       INTEGER            I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K,
-     $                   KMIN, N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE
+     $                   KMIN, N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE,
+     $                   I1, N1
       REAL               D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN,
      $                   DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX,
      $                   QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL,
-     $                   TOL2, TRACE, ZMAX
+     $                   TOL2, TRACE, ZMAX, TEMPE, TEMPQ
 *     ..
 *     .. External Subroutines ..
       EXTERNAL           SLASQ3, SLASRT, XERBLA
@@ -401,7 +403,7 @@
 *               and that the tests for deflation upon entry in SLASQ3 
 *               should not be performed.
 *
-         NBIG = 30*( N0-I0+1 )
+         NBIG = 100*( N0-I0+1 )
          DO 140 IWHILB = 1, NBIG
             IF( I0.GT.N0 ) 
      $         GO TO 150
@@ -446,6 +448,49 @@
   140    CONTINUE
 *
          INFO = 2
+*       
+*        Maximum number of iterations exceeded, restore the shift 
+*        SIGMA and place the new d's and e's in a qd array.
+*        This might need to be done for several blocks
+*
+         I1 = I0
+         N1 = N0
+ 145     CONTINUE
+         TEMPQ = Z( 4*I0-3 )
+         Z( 4*I0-3 ) = Z( 4*I0-3 ) + SIGMA
+         DO K = I0+1, N0
+            TEMPE = Z( 4*K-5 )
+            Z( 4*K-5 ) = Z( 4*K-5 ) * (TEMPQ / Z( 4*K-7 ))
+            TEMPQ = Z( 4*K-3 )
+            Z( 4*K-3 ) = Z( 4*K-3 ) + SIGMA + TEMPE - Z( 4*K-5 )
+         END DO
+*
+*        Prepare to do this on the previous block if there is one
+*
+         IF( I1.GT.1 ) THEN
+            N1 = I1-1
+            DO WHILE( ( I1.GE.2 ) .AND. ( Z(4*I-5).GE.ZERO ) )
+               I1 = I1 - 1
+            END DO
+            IF( I1.GE.1 ) THEN
+               SIGMA = -Z(4*N1-1)
+               GO TO 145
+            END IF
+         END IF
+
+         DO K = 1, N
+            Z( 2*K-1 ) = Z( 4*K-3 )
+*
+*        Only the block 1..N0 is unfinished.  The rest of the e's
+*        must be essentially zero, although sometimes other data
+*        has been stored in them.
+*
+            IF( K.LT.N0 ) THEN
+               Z( 2*K ) = Z( 4*K-1 )
+            ELSE
+               Z( 2*K ) = 0
+            END IF
+         END DO
          RETURN
 *
 *        end IWHILB
diff --git a/SRC/zbdsqr.f b/SRC/zbdsqr.f
index 4bbc5b7b..6741bd25 100644
--- a/SRC/zbdsqr.f
+++ b/SRC/zbdsqr.f
@@ -106,7 +106,8 @@
 *          The leading dimension of the array C.
 *          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
 *
-*  RWORK   (workspace) DOUBLE PRECISION array, dimension (4*N)
+*  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
+*          if NCVT = NRU = NCC = 0, (max(1, 4*N-4)) otherwise
 *
 *  INFO    (output) INTEGER
 *          = 0:  successful exit
@@ -223,7 +224,11 @@
 *
       IF( .NOT.ROTATE ) THEN
          CALL DLASQ1( N, D, E, RWORK, INFO )
-         RETURN
+*
+*     If INFO equals 2, dqds didn't finish, try to finish
+*         
+         IF( INFO .NE. 2 ) RETURN
+         INFO = 0
       END IF
 *
       NM1 = N - 1