Created
February 1, 2021 20:03
-
-
Save omus/c25ec6e3f871a062799a4275b40bb448 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
commit c4b5abbe43d7c22215ef36ef4f7c1413c975678c | |
Author: Martin Kroeker <martin@ruby.chemie.uni-freiburg.de> | |
Date: Fri Jan 29 10:45:36 2021 +0100 | |
fix data type | |
commit f87842483eee9d158f44d51d4c09662c3cff7526 | |
Author: Martin Kroeker <martin@ruby.chemie.uni-freiburg.de> | |
Date: Fri Jan 29 09:56:12 2021 +0100 | |
fix calculation of non-exceptional shift (from Reference-LAPACK PR 477) | |
--- | |
diff --git a/lapack-netlib/SRC/chgeqz.f b/lapack-netlib/SRC/chgeqz.f | |
index 1616840e..4725e716 100644 | |
--- a/lapack-netlib/SRC/chgeqz.f | |
+++ b/lapack-netlib/SRC/chgeqz.f | |
@@ -320,12 +320,13 @@ | |
$ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP | |
COMPLEX ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2, | |
$ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1, | |
- $ U12, X | |
+ $ U12, X, ABI12, Y | |
* .. | |
* .. External Functions .. | |
+ COMPLEX CLADIV | |
LOGICAL LSAME | |
REAL CLANHS, SLAMCH | |
- EXTERNAL LSAME, CLANHS, SLAMCH | |
+ EXTERNAL CLADIV, LLSAME, CLANHS, SLAMCH | |
* .. | |
* .. External Subroutines .. | |
EXTERNAL CLARTG, CLASET, CROT, CSCAL, XERBLA | |
@@ -729,15 +730,21 @@ | |
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 | |
* | |
diff --git a/lapack-netlib/SRC/zhgeqz.f b/lapack-netlib/SRC/zhgeqz.f | |
index b21199e9..b28ae47a 100644 | |
--- a/lapack-netlib/SRC/zhgeqz.f | |
+++ b/lapack-netlib/SRC/zhgeqz.f | |
@@ -320,12 +320,13 @@ | |
$ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP | |
COMPLEX*16 ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2, | |
$ CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1, | |
- $ U12, X | |
+ $ 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 | |
@@ -730,15 +731,21 @@ | |
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 | |
* |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment