diff options
Diffstat (limited to 'SRC/dlasq2.f')
| -rw-r--r-- | SRC/dlasq2.f | 50 |
1 files changed, 46 insertions, 4 deletions
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 |