--- rpl/lapack/lapack/dhgeqz.f 2012/07/31 11:06:35 1.11
+++ rpl/lapack/lapack/dhgeqz.f 2017/06/17 10:53:51 1.18
@@ -2,18 +2,18 @@
*
* =========== DOCUMENTATION ===========
*
-* Online html documentation available at
-* http://www.netlib.org/lapack/explore-html/
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
-*> Download DHGEQZ + dependencies
-*>
-*> [TGZ]
-*>
-*> [ZIP]
-*>
+*> Download DHGEQZ + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
*> [TXT]
-*> \endhtmlonly
+*> \endhtmlonly
*
* Definition:
* ===========
@@ -21,7 +21,7 @@
* SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
* ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
* LWORK, INFO )
-*
+*
* .. Scalar Arguments ..
* CHARACTER COMPQ, COMPZ, JOB
* INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
@@ -31,7 +31,7 @@
* $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
* $ WORK( * ), Z( LDZ, * )
* ..
-*
+*
*
*> \par Purpose:
* =============
@@ -50,9 +50,9 @@
*>
*> If JOB='S', then the Hessenberg-triangular pair (H,T) is
*> also reduced to generalized Schur form,
-*>
+*>
*> H = Q*S*Z**T, T = Q*P*Z**T,
-*>
+*>
*> where Q and Z are orthogonal matrices, P is an upper triangular
*> matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
*> diagonal blocks.
@@ -75,7 +75,7 @@
*> generalized Schur factorization of (A,B):
*>
*> A = (Q1*Q)*S*(Z1*Z)**T, B = (Q1*Q)*P*(Z1*Z)**T.
-*>
+*>
*> To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
*> of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
*> complex and beta real.
@@ -86,7 +86,7 @@
*> alternate form of the GNEP
*> mu*A*y = B*y.
*> Real eigenvalues can be read directly from the generalized Schur
-*> form:
+*> form:
*> alpha = S(i,i), beta = P(i,i).
*>
*> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
@@ -101,7 +101,7 @@
*> \verbatim
*> JOB is CHARACTER*1
*> = 'E': Compute eigenvalues only;
-*> = 'S': Compute eigenvalues and the Schur form.
+*> = 'S': Compute eigenvalues and the Schur form.
*> \endverbatim
*>
*> \param[in] COMPQ
@@ -211,12 +211,12 @@
*> \param[in,out] Q
*> \verbatim
*> Q is DOUBLE PRECISION array, dimension (LDQ, N)
-*> On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in
+*> On entry, if COMPQ = 'V', the orthogonal matrix Q1 used in
*> the reduction of (A,B) to generalized Hessenberg form.
-*> On exit, if COMPZ = 'I', the orthogonal matrix of left Schur
-*> vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix
+*> On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
+*> vectors of (H,T), and if COMPQ = 'V', the orthogonal matrix
*> of left Schur vectors of (A,B).
-*> Not referenced if COMPZ = 'N'.
+*> Not referenced if COMPQ = 'N'.
*> \endverbatim
*>
*> \param[in] LDQ
@@ -277,12 +277,12 @@
* Authors:
* ========
*
-*> \author Univ. of Tennessee
-*> \author Univ. of California Berkeley
-*> \author Univ. of Colorado Denver
-*> \author NAG Ltd.
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
*
-*> \date April 2012
+*> \date June 2016
*
*> \ingroup doubleGEcomputational
*
@@ -304,10 +304,10 @@
$ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
$ LWORK, INFO )
*
-* -- LAPACK computational routine (version 3.4.1) --
+* -- LAPACK computational routine (version 3.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* April 2012
+* June 2016
*
* .. Scalar Arguments ..
CHARACTER COMPQ, COMPZ, JOB
@@ -739,9 +739,9 @@
* Exceptional shift. Chosen for no particularly good reason.
* (Single shift only.)
*
- IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT.
+ IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST, ILAST-1 ) ).LT.
$ ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
- ESHIFT = ESHIFT + H( ILAST, ILAST-1 ) /
+ ESHIFT = H( ILAST, ILAST-1 ) /
$ T( ILAST-1, ILAST-1 )
ELSE
ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
@@ -759,6 +759,16 @@
$ T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
$ S2, WR, WR2, WI )
*
+ IF ( ABS( (WR/S1)*T( ILAST, ILAST ) - H( ILAST, ILAST ) )
+ $ .GT. ABS( (WR2/S2)*T( ILAST, ILAST )
+ $ - H( ILAST, ILAST ) ) ) THEN
+ TEMP = WR
+ WR = WR2
+ WR2 = TEMP
+ TEMP = S1
+ S1 = S2
+ S2 = TEMP
+ END IF
TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )
IF( WI.NE.ZERO )
$ GO TO 200