Skip to content

Commit

Permalink
Fix bug of *kteqr about convergence when subdiagonal elements are equ…
Browse files Browse the repository at this point in the history 
…al or elements are small
  • Loading branch information
@sh-zheng
sh-zheng committed Mar 15, 2024
1 parent 6b076a4 commit d7c7829
Show file tree
Hide file tree
Showing 2 changed files with 118 additions and 18 deletions.
68 changes: 59 additions & 9 deletions SRC/dkteqr.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,8 @@ SUBROUTINE DKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
$ LENDM1, LENDP1, LENDSV, LM3, LSV, M, MM, MM1,
$ NM1, NMAXIT
DOUBLE PRECISION ANORM, B, EPS, EPS2, P, R, VA, VB, VC, VD, E3,
$ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST, TEMP
$ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST, TEMP,
$ SQRT2
* ..
* .. External Functions ..
LOGICAL LSAME
Expand Down Expand Up @@ -237,6 +238,7 @@ SUBROUTINE DKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
SAFMAX = ONE / SAFMIN
SSFMAX = SQRT( SAFMAX ) / THREE
SSFMIN = SQRT( SAFMIN ) / EPS2
SQRT2 = SQRT( TWO )
*
* Compute the eigenvalues and eigenvectors of the tridiagonal
* matrix.
Expand Down Expand Up @@ -371,12 +373,36 @@ SUBROUTINE DKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
R = E(M-1)*E(M-2)
S = SIGN(DLAPY2( P, R ), P)
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
VD = ZERO
E(M-1) = ZERO
ELSEIF(E(M-2).EQ.E(M-1)) THEN
E(M-2) = SQRT2*E(M-2)
E(M-1) = ZERO
IF( ICOMPZ.GT.0 ) THEN
DO J = 1, N
TEMP = Z( J, M-2 )
Z( J, M-2 ) = -Z( J, M-1 )
Z( J, M-1 ) = (SQRT2*TEMP - SQRT2*Z(J, M))/TWO
Z( J, M ) = (SQRT2*TEMP + SQRT2*Z(J, M))/TWO
END DO
END IF
GO TO 40
ELSEIF(E(M-2).EQ.-E(M-1)) THEN
E(M-2) = SQRT2*E(M-2)
E(M-1) = ZERO
IF( ICOMPZ.GT.0 ) THEN
DO J = 1, N
TEMP = Z( J, M-2 )
Z( J, M-2 ) = -Z( J, M-1 )
Z( J, M-1 ) = (SQRT2*TEMP + SQRT2*Z(J, M))/TWO
Z( J, M ) = (-SQRT2*TEMP + SQRT2*Z(J, M))/TWO
END DO
END IF
GO TO 40
ELSE
VA = -P/S
VB = -VA + (R/S)*(R/(P+S))
Expand Down Expand Up @@ -426,7 +452,7 @@ SUBROUTINE DKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
R = E(M-1)*E(M-2)
S = SIGN(DLAPY2( P, R ), P)
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down Expand Up @@ -470,7 +496,7 @@ SUBROUTINE DKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
S = SIGN(DLAPY2( P, R ), P)
E(I) = -S
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down Expand Up @@ -513,7 +539,7 @@ SUBROUTINE DKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
S = SIGN(DLAPY2( P, R ), P)
E(I) = -S
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down Expand Up @@ -624,12 +650,36 @@ SUBROUTINE DKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
R = E(M)*E(M+1)
S = SIGN(DLAPY2( P, R ), P)
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
VD = ZERO
E(M) = ZERO
ELSEIF(E(M).EQ.E(M+1)) THEN
E(M) = SQRT2*E(M)
E(M+1) = ZERO
IF( ICOMPZ.GT.0 ) THEN
DO J = 1, N
TEMP = Z( J, M )
Z( J, M ) = -Z( J, M+1 )
Z( J, M+1 ) = (SQRT2*TEMP - SQRT2*Z(J, M+2))/TWO
Z( J, M+2 ) = (SQRT2*TEMP + SQRT2*Z(J, M+2))/TWO
END DO
END IF
GO TO 90
ELSEIF(E(M).EQ.-E(M+1)) THEN
E(M) = SQRT2*E(M)
E(M+1) = ZERO
IF( ICOMPZ.GT.0 ) THEN
DO J = 1, N
TEMP = Z( J, M )
Z( J, M ) = -Z( J, M+1 )
Z( J, M+1 ) = (SQRT2*TEMP + SQRT2*Z(J, M+2))/TWO
Z( J, M+2 ) = (-SQRT2*TEMP + SQRT2*Z(J, M+2))/TWO
END DO
END IF
GO TO 90
ELSE
VA = -P/S
VB = -VA + (R/S)*(R/(P+S))
Expand Down Expand Up @@ -679,7 +729,7 @@ SUBROUTINE DKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
R = E(M)*E(M+1)
S = SIGN(DLAPY2( P, R ), P)
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down Expand Up @@ -723,7 +773,7 @@ SUBROUTINE DKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
S = SIGN(DLAPY2( P, R ), P)
E(I-1) = -S
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down Expand Up @@ -766,7 +816,7 @@ SUBROUTINE DKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
S = SIGN(DLAPY2( P, R ), P)
E(I-1) = -S
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down
68 changes: 59 additions & 9 deletions SRC/skteqr.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,8 @@ SUBROUTINE SKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
$ LENDM1, LENDP1, LENDSV, LM3, LSV, M, MM, MM1,
$ NM1, NMAXIT
REAL ANORM, B, EPS, EPS2, P, R, VA, VB, VC, VD, E3,
$ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST, TEMP
$ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST, TEMP,
$ SQRT2
* ..
* .. External Functions ..
LOGICAL LSAME
Expand Down Expand Up @@ -237,6 +238,7 @@ SUBROUTINE SKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
SAFMAX = ONE / SAFMIN
SSFMAX = SQRT( SAFMAX ) / THREE
SSFMIN = SQRT( SAFMIN ) / EPS2
SQRT2 = SQRT( TWO )
*
* Compute the eigenvalues and eigenvectors of the tridiagonal
* matrix.
Expand Down Expand Up @@ -371,12 +373,36 @@ SUBROUTINE SKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
R = E(M-1)*E(M-2)
S = SIGN(SLAPY2( P, R ), P)
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
VD = ZERO
E(M-1) = ZERO
ELSEIF(E(M-2).EQ.E(M-1)) THEN
E(M-2) = SQRT2*E(M-2)
E(M-1) = ZERO
IF( ICOMPZ.GT.0 ) THEN
DO J = 1, N
TEMP = Z( J, M-2 )
Z( J, M-2 ) = -Z( J, M-1 )
Z( J, M-1 ) = (SQRT2*TEMP - SQRT2*Z(J, M))/TWO
Z( J, M ) = (SQRT2*TEMP + SQRT2*Z(J, M))/TWO
END DO
END IF
GO TO 40
ELSEIF(E(M-2).EQ.-E(M-1)) THEN
E(M-2) = SQRT2*E(M-2)
E(M-1) = ZERO
IF( ICOMPZ.GT.0 ) THEN
DO J = 1, N
TEMP = Z( J, M-2 )
Z( J, M-2 ) = -Z( J, M-1 )
Z( J, M-1 ) = (SQRT2*TEMP + SQRT2*Z(J, M))/TWO
Z( J, M ) = (-SQRT2*TEMP + SQRT2*Z(J, M))/TWO
END DO
END IF
GO TO 40
ELSE
VA = -P/S
VB = -VA + (R/S)*(R/(P+S))
Expand Down Expand Up @@ -426,7 +452,7 @@ SUBROUTINE SKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
R = E(M-1)*E(M-2)
S = SIGN(SLAPY2( P, R ), P)
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down Expand Up @@ -470,7 +496,7 @@ SUBROUTINE SKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
S = SIGN(SLAPY2( P, R ), P)
E(I) = -S
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down Expand Up @@ -513,7 +539,7 @@ SUBROUTINE SKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
S = SIGN(SLAPY2( P, R ), P)
E(I) = -S
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down Expand Up @@ -624,12 +650,36 @@ SUBROUTINE SKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
R = E(M)*E(M+1)
S = SIGN(SLAPY2( P, R ), P)
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
VD = ZERO
E(M) = ZERO
ELSEIF(E(M).EQ.E(M+1)) THEN
E(M) = SQRT2*E(M)
E(M+1) = ZERO
IF( ICOMPZ.GT.0 ) THEN
DO J = 1, N
TEMP = Z( J, M )
Z( J, M ) = -Z( J, M+1 )
Z( J, M+1 ) = (SQRT2*TEMP - SQRT2*Z(J, M+2))/TWO
Z( J, M+2 ) = (SQRT2*TEMP + SQRT2*Z(J, M+2))/TWO
END DO
END IF
GO TO 90
ELSEIF(E(M).EQ.-E(M+1)) THEN
E(M) = SQRT2*E(M)
E(M+1) = ZERO
IF( ICOMPZ.GT.0 ) THEN
DO J = 1, N
TEMP = Z( J, M )
Z( J, M ) = -Z( J, M+1 )
Z( J, M+1 ) = (SQRT2*TEMP + SQRT2*Z(J, M+2))/TWO
Z( J, M+2 ) = (-SQRT2*TEMP + SQRT2*Z(J, M+2))/TWO
END DO
END IF
GO TO 90
ELSE
VA = -P/S
VB = -VA + (R/S)*(R/(P+S))
Expand Down Expand Up @@ -679,7 +729,7 @@ SUBROUTINE SKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
R = E(M)*E(M+1)
S = SIGN(SLAPY2( P, R ), P)
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down Expand Up @@ -723,7 +773,7 @@ SUBROUTINE SKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
S = SIGN(SLAPY2( P, R ), P)
E(I-1) = -S
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down Expand Up @@ -766,7 +816,7 @@ SUBROUTINE SKTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
S = SIGN(SLAPY2( P, R ), P)
E(I-1) = -S
*
IF(S*(P+S).EQ.ZERO) THEN
IF(S.EQ.ZERO .OR. R/S.EQ.ZERO) THEN
VA = -ONE
VB = ONE
VC = ONE
Expand Down

0 comments on commit d7c7829

Please sign in to comment.