Skip to content

add extra exceptional shift to solve rare convergence issues #477

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 7 commits into from
Feb 9, 2021
Merged

add extra exceptional shift to solve rare convergence issues #477

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
d789f2e
add extra exceptional shift to solve rare convergence issues
thijssteel Jan 26, 2021
b8ab53f
some extra overflow protection (better safe than sorry)
thijssteel Jan 26, 2021
441b21e
remove unnecessary overflow check
thijssteel Jan 26, 2021
4e89f98
fix error in calculation of shifts in complex single shift qz
thijssteel Jan 29, 2021
c0964cf
fix type error in chgeqz
thijssteel Jan 29, 2021
5a893d0
fix typo in chgeqz
thijssteel Jan 30, 2021
6bb0869
fix compiler warnings in qz code
thijssteel Jan 30, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 29 additions & 13 deletions SRC/chgeqz.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -319,13 +319,14 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
REAL ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
$ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
COMPLEX ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
$ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,
$ U12, X
$ CTEMP3, ESHIFT, S, SHIFT, SIGNBC,
$ U12, X, ABI12, Y
* ..
* .. External Functions ..
COMPLEX CLADIV
LOGICAL LSAME
REAL CLANHS, SLAMCH
EXTERNAL LSAME, CLANHS, SLAMCH
EXTERNAL CLADIV, LSAME, CLANHS, SLAMCH
* ..
* .. External Subroutines ..
EXTERNAL CLARTG, CLASET, CROT, CSCAL, XERBLA
Expand All @@ -350,6 +351,7 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
ILSCHR = .TRUE.
ISCHUR = 2
ELSE
ILSCHR = .TRUE.
ISCHUR = 0
END IF
*
Expand All @@ -363,6 +365,7 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
ILQ = .TRUE.
ICOMPQ = 3
ELSE
ILQ = .TRUE.
ICOMPQ = 0
END IF
*
Expand All @@ -376,6 +379,7 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
ILZ = .TRUE.
ICOMPZ = 3
ELSE
ILZ = .TRUE.
ICOMPZ = 0
END IF
*
Expand Down Expand Up @@ -729,22 +733,34 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
$ ( BSCALE*T( ILAST, ILAST ) )
ABI22 = AD22 - U12*AD21
ABI12 = AD12 - U12*AD11
*
T1 = HALF*( AD11+ABI22 )
RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
TEMP = REAL( T1-ABI22 )*REAL( RTDISC ) +
$ AIMAG( T1-ABI22 )*AIMAG( RTDISC )
IF( TEMP.LE.ZERO ) THEN
SHIFT = T1 + RTDISC
ELSE
SHIFT = T1 - RTDISC
SHIFT = ABI22
CTEMP = SQRT( ABI12 )*SQRT( AD21 )
TEMP = ABS1( CTEMP )
IF( CTEMP.NE.ZERO ) THEN
X = HALF*( AD11-SHIFT )
TEMP2 = ABS1( X )
TEMP = MAX( TEMP, ABS1( X ) )
Y = TEMP*SQRT( ( X / TEMP )**2+( CTEMP / TEMP )**2 )
IF( TEMP2.GT.ZERO ) THEN
IF( REAL( X / TEMP2 )*REAL( Y )+
$ AIMAG( X / TEMP2 )*AIMAG( Y ).LT.ZERO )Y = -Y
END IF
SHIFT = SHIFT - CTEMP*CLADIV( CTEMP, ( X+Y ) )
END IF
ELSE
*
* Exceptional shift. Chosen for no particularly good reason.
*
ESHIFT = ESHIFT + (ASCALE*H(ILAST,ILAST-1))/
$ (BSCALE*T(ILAST-1,ILAST-1))
IF( ( IITER / 20 )*20.EQ.IITER .AND.
$ BSCALE*ABS1(T( ILAST, ILAST )).GT.SAFMIN ) THEN
ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
$ ILAST ) )/( BSCALE*T( ILAST, ILAST ) )
ELSE
ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
$ ILAST-1 ) )/( BSCALE*T( ILAST-1, ILAST-1 ) )
END IF
SHIFT = ESHIFT
END IF
*
Expand Down
42 changes: 29 additions & 13 deletions SRC/zhgeqz.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -319,13 +319,14 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
DOUBLE PRECISION ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
$ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
COMPLEX*16 ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
$ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,
$ U12, X
$ CTEMP3, ESHIFT, S, SHIFT, SIGNBC,
$ U12, X, ABI12, Y
* ..
* .. External Functions ..
COMPLEX*16 ZLADIV
LOGICAL LSAME
DOUBLE PRECISION DLAMCH, ZLANHS
EXTERNAL LSAME, DLAMCH, ZLANHS
EXTERNAL ZLADIV, LSAME, DLAMCH, ZLANHS
* ..
* .. External Subroutines ..
EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL
Expand All @@ -351,6 +352,7 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
ILSCHR = .TRUE.
ISCHUR = 2
ELSE
ILSCHR = .TRUE.
ISCHUR = 0
END IF
*
Expand All @@ -364,6 +366,7 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
ILQ = .TRUE.
ICOMPQ = 3
ELSE
ILQ = .TRUE.
ICOMPQ = 0
END IF
*
Expand All @@ -377,6 +380,7 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
ILZ = .TRUE.
ICOMPZ = 3
ELSE
ILZ = .TRUE.
ICOMPZ = 0
END IF
*
Expand Down Expand Up @@ -730,22 +734,34 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
$ ( BSCALE*T( ILAST, ILAST ) )
ABI22 = AD22 - U12*AD21
ABI12 = AD12 - U12*AD11
*
T1 = HALF*( AD11+ABI22 )
RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
TEMP = DBLE( T1-ABI22 )*DBLE( RTDISC ) +
$ DIMAG( T1-ABI22 )*DIMAG( RTDISC )
IF( TEMP.LE.ZERO ) THEN
SHIFT = T1 + RTDISC
ELSE
SHIFT = T1 - RTDISC
SHIFT = ABI22
CTEMP = SQRT( ABI12 )*SQRT( AD21 )
TEMP = ABS1( CTEMP )
IF( CTEMP.NE.ZERO ) THEN
X = HALF*( AD11-SHIFT )
TEMP2 = ABS1( X )
TEMP = MAX( TEMP, ABS1( X ) )
Y = TEMP*SQRT( ( X / TEMP )**2+( CTEMP / TEMP )**2 )
IF( TEMP2.GT.ZERO ) THEN
IF( DBLE( X / TEMP2 )*DBLE( Y )+
$ DIMAG( X / TEMP2 )*DIMAG( Y ).LT.ZERO )Y = -Y
END IF
SHIFT = SHIFT - CTEMP*ZLADIV( CTEMP, ( X+Y ) )
END IF
ELSE
*
* Exceptional shift. Chosen for no particularly good reason.
*
ESHIFT = ESHIFT + (ASCALE*H(ILAST,ILAST-1))/
$ (BSCALE*T(ILAST-1,ILAST-1))
IF( ( IITER / 20 )*20.EQ.IITER .AND.
$ BSCALE*ABS1(T( ILAST, ILAST )).GT.SAFMIN ) THEN
ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
$ ILAST ) )/( BSCALE*T( ILAST, ILAST ) )
ELSE
ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
$ ILAST-1 ) )/( BSCALE*T( ILAST-1, ILAST-1 ) )
END IF
SHIFT = ESHIFT
END IF
*
Expand Down