Skip to content

Commit e56b31e

Browse files
langoulangou
authored
Merge pull request #477 from thijssteel/fix-475-eigenvalue-convergence-failure
add extra exceptional shift to solve rare convergence issues
2 parents 2b9f0f5 + 6bb0869 commit e56b31e

File tree

2 files changed

+58
-26
lines changed

2 files changed

+58
-26
lines changed

‎SRC/chgeqz.f

Lines changed: 29 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -319,13 +319,14 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
319319
REAL ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
320320
$ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
321321
COMPLEX ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
322-
$ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,
323-
$ U12, X
322+
$ CTEMP3, ESHIFT, S, SHIFT, SIGNBC,
323+
$ U12, X, ABI12, Y
324324
* ..
325325
* .. External Functions ..
326+
COMPLEX CLADIV
326327
LOGICAL LSAME
327328
REAL CLANHS, SLAMCH
328-
EXTERNAL LSAME, CLANHS, SLAMCH
329+
EXTERNAL CLADIV, LSAME, CLANHS, SLAMCH
329330
* ..
330331
* .. External Subroutines ..
331332
EXTERNAL CLARTG, CLASET, CROT, CSCAL, XERBLA
@@ -350,6 +351,7 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
350351
ILSCHR = .TRUE.
351352
ISCHUR = 2
352353
ELSE
354+
ILSCHR = .TRUE.
353355
ISCHUR = 0
354356
END IF
355357
*
@@ -363,6 +365,7 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
363365
ILQ = .TRUE.
364366
ICOMPQ = 3
365367
ELSE
368+
ILQ = .TRUE.
366369
ICOMPQ = 0
367370
END IF
368371
*
@@ -376,6 +379,7 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
376379
ILZ = .TRUE.
377380
ICOMPZ = 3
378381
ELSE
382+
ILZ = .TRUE.
379383
ICOMPZ = 0
380384
END IF
381385
*
@@ -729,22 +733,34 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
729733
AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
730734
$ ( BSCALE*T( ILAST, ILAST ) )
731735
ABI22 = AD22 - U12*AD21
736+
ABI12 = AD12 - U12*AD11
732737
*
733-
T1 = HALF*( AD11+ABI22 )
734-
RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
735-
TEMP = REAL( T1-ABI22 )*REAL( RTDISC ) +
736-
$ AIMAG( T1-ABI22 )*AIMAG( RTDISC )
737-
IF( TEMP.LE.ZERO ) THEN
738-
SHIFT = T1 + RTDISC
739-
ELSE
740-
SHIFT = T1 - RTDISC
738+
SHIFT = ABI22
739+
CTEMP = SQRT( ABI12 )*SQRT( AD21 )
740+
TEMP = ABS1( CTEMP )
741+
IF( CTEMP.NE.ZERO ) THEN
742+
X = HALF*( AD11-SHIFT )
743+
TEMP2 = ABS1( X )
744+
TEMP = MAX( TEMP, ABS1( X ) )
745+
Y = TEMP*SQRT( ( X / TEMP )**2+( CTEMP / TEMP )**2 )
746+
IF( TEMP2.GT.ZERO ) THEN
747+
IF( REAL( X / TEMP2 )*REAL( Y )+
748+
$ AIMAG( X / TEMP2 )*AIMAG( Y ).LT.ZERO )Y = -Y
749+
END IF
750+
SHIFT = SHIFT - CTEMP*CLADIV( CTEMP, ( X+Y ) )
741751
END IF
742752
ELSE
743753
*
744754
* Exceptional shift. Chosen for no particularly good reason.
745755
*
746-
ESHIFT = ESHIFT + (ASCALE*H(ILAST,ILAST-1))/
747-
$ (BSCALE*T(ILAST-1,ILAST-1))
756+
IF( ( IITER / 20 )*20.EQ.IITER .AND.
757+
$ BSCALE*ABS1(T( ILAST, ILAST )).GT.SAFMIN ) THEN
758+
ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
759+
$ ILAST ) )/( BSCALE*T( ILAST, ILAST ) )
760+
ELSE
761+
ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
762+
$ ILAST-1 ) )/( BSCALE*T( ILAST-1, ILAST-1 ) )
763+
END IF
748764
SHIFT = ESHIFT
749765
END IF
750766
*

‎SRC/zhgeqz.f

Lines changed: 29 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -319,13 +319,14 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
319319
DOUBLE PRECISION ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
320320
$ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
321321
COMPLEX*16 ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
322-
$ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,
323-
$ U12, X
322+
$ CTEMP3, ESHIFT, S, SHIFT, SIGNBC,
323+
$ U12, X, ABI12, Y
324324
* ..
325325
* .. External Functions ..
326+
COMPLEX*16 ZLADIV
326327
LOGICAL LSAME
327328
DOUBLE PRECISION DLAMCH, ZLANHS
328-
EXTERNAL LSAME, DLAMCH, ZLANHS
329+
EXTERNAL ZLADIV, LSAME, DLAMCH, ZLANHS
329330
* ..
330331
* .. External Subroutines ..
331332
EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL
@@ -351,6 +352,7 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
351352
ILSCHR = .TRUE.
352353
ISCHUR = 2
353354
ELSE
355+
ILSCHR = .TRUE.
354356
ISCHUR = 0
355357
END IF
356358
*
@@ -364,6 +366,7 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
364366
ILQ = .TRUE.
365367
ICOMPQ = 3
366368
ELSE
369+
ILQ = .TRUE.
367370
ICOMPQ = 0
368371
END IF
369372
*
@@ -377,6 +380,7 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
377380
ILZ = .TRUE.
378381
ICOMPZ = 3
379382
ELSE
383+
ILZ = .TRUE.
380384
ICOMPZ = 0
381385
END IF
382386
*
@@ -730,22 +734,34 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
730734
AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
731735
$ ( BSCALE*T( ILAST, ILAST ) )
732736
ABI22 = AD22 - U12*AD21
737+
ABI12 = AD12 - U12*AD11
733738
*
734-
T1 = HALF*( AD11+ABI22 )
735-
RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
736-
TEMP = DBLE( T1-ABI22 )*DBLE( RTDISC ) +
737-
$ DIMAG( T1-ABI22 )*DIMAG( RTDISC )
738-
IF( TEMP.LE.ZERO ) THEN
739-
SHIFT = T1 + RTDISC
740-
ELSE
741-
SHIFT = T1 - RTDISC
739+
SHIFT = ABI22
740+
CTEMP = SQRT( ABI12 )*SQRT( AD21 )
741+
TEMP = ABS1( CTEMP )
742+
IF( CTEMP.NE.ZERO ) THEN
743+
X = HALF*( AD11-SHIFT )
744+
TEMP2 = ABS1( X )
745+
TEMP = MAX( TEMP, ABS1( X ) )
746+
Y = TEMP*SQRT( ( X / TEMP )**2+( CTEMP / TEMP )**2 )
747+
IF( TEMP2.GT.ZERO ) THEN
748+
IF( DBLE( X / TEMP2 )*DBLE( Y )+
749+
$ DIMAG( X / TEMP2 )*DIMAG( Y ).LT.ZERO )Y = -Y
750+
END IF
751+
SHIFT = SHIFT - CTEMP*ZLADIV( CTEMP, ( X+Y ) )
742752
END IF
743753
ELSE
744754
*
745755
* Exceptional shift. Chosen for no particularly good reason.
746756
*
747-
ESHIFT = ESHIFT + (ASCALE*H(ILAST,ILAST-1))/
748-
$ (BSCALE*T(ILAST-1,ILAST-1))
757+
IF( ( IITER / 20 )*20.EQ.IITER .AND.
758+
$ BSCALE*ABS1(T( ILAST, ILAST )).GT.SAFMIN ) THEN
759+
ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
760+
$ ILAST ) )/( BSCALE*T( ILAST, ILAST ) )
761+
ELSE
762+
ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
763+
$ ILAST-1 ) )/( BSCALE*T( ILAST-1, ILAST-1 ) )
764+
END IF
749765
SHIFT = ESHIFT
750766
END IF
751767
*

0 commit comments

Comments
 (0)