Skip to content

Instantly share code, notes, and snippets.

@pablosanjose

pablosanjose/openblas-exshift.patch

Last active January 30, 2021 16:12
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save pablosanjose/3b673d14865ac7c7bdef194feca1c966 to your computer and use it in GitHub Desktop.
Save pablosanjose/3b673d14865ac7c7bdef194feca1c966 to your computer and use it in GitHub Desktop.
Download ZIP
Patch for OpenBLAS@0.3.13
Raw
openblas-exshift.patch
diff --git a/lapack-netlib/SRC/chgeqz.f b/lapack-netlib/SRC/chgeqz.f
index 73d35621c..0351c79e1 100644
--- a/lapack-netlib/SRC/chgeqz.f
+++ b/lapack-netlib/SRC/chgeqz.f
@@ -320,12 +320,13 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
$ 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
* ..diff --git a/lapack-netlib/SRC/chgeqz.f b/lapack-netlib/SRC/chgeqz.f
index 1616840ec..4725e7169 100644
--- a/lapack-netlib/SRC/chgeqz.f
+++ b/lapack-netlib/SRC/chgeqz.f
@@ -320,12 +320,13 @@ SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
$ 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, LSAME, CLANHS, SLAMCH
* ..
* .. External Subroutines ..
EXTERNAL CLARTG, CLASET, CROT, CSCAL, XERBLA
@@ -729,15 +730,21 @@ 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
*
diff --git a/lapack-netlib/SRC/zhgeqz.f b/lapack-netlib/SRC/zhgeqz.f
index b21199e9e..b28ae47a4 100644
--- a/lapack-netlib/SRC/zhgeqz.f
+++ b/lapack-netlib/SRC/zhgeqz.f
@@ -320,12 +320,13 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
$ 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 @@ 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
*
* .. 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
@@ -729,22 +730,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
*
diff --git a/lapack-netlib/SRC/zhgeqz.f b/lapack-netlib/SRC/zhgeqz.f
index b51cba4f7..b28ae47a4 100644
--- a/lapack-netlib/SRC/zhgeqz.f
+++ b/lapack-netlib/SRC/zhgeqz.f
@@ -320,12 +320,13 @@ SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
$ 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,22 +731,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
*
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment