summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorlipshitz <lipshitz@8a072113-8704-0410-8d35-dd094bca7971>2011-10-03 19:24:45 +0000
committerlipshitz <lipshitz@8a072113-8704-0410-8d35-dd094bca7971>2011-10-03 19:24:45 +0000
commit3b27d3196c299551e634bd40e0662d244ff693e1 (patch)
treec88b0f664840b9227fc9f4a7f11f725d8d4bb884
parent0cd8057636e95b333c27913ee8d0dd3143f8af73 (diff)
downloadlapack-3b27d3196c299551e634bd40e0662d244ff693e1.tar.gz
lapack-3b27d3196c299551e634bd40e0662d244ff693e1.tar.bz2
lapack-3b27d3196c299551e634bd40e0662d244ff693e1.zip
Fix to Ming Gu's bug in dqds. See post of LAPACK DEV message board for more details
-rw-r--r--SRC/cbdsqr.f9
-rw-r--r--SRC/dbdsqr.f6
-rw-r--r--SRC/dlasq1.f18
-rw-r--r--SRC/dlasq2.f50
-rw-r--r--SRC/sbdsqr.f6
-rw-r--r--SRC/slasq1.f18
-rw-r--r--SRC/slasq2.f55
-rw-r--r--SRC/zbdsqr.f9
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