--- rpl/lapack/lapack/ztgex2.f 2010/04/21 13:45:39 1.2
+++ rpl/lapack/lapack/ztgex2.f 2023/08/07 08:39:40 1.21
@@ -1,10 +1,196 @@
+*> \brief \b ZTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equivalence transformation.
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZTGEX2 + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
+* LDZ, J1, INFO )
+*
+* .. Scalar Arguments ..
+* LOGICAL WANTQ, WANTZ
+* INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
+* $ Z( LDZ, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
+*> in an upper triangular matrix pair (A, B) by an unitary equivalence
+*> transformation.
+*>
+*> (A, B) must be in generalized Schur canonical form, that is, A and
+*> B are both upper triangular.
+*>
+*> Optionally, the matrices Q and Z of generalized Schur vectors are
+*> updated.
+*>
+*> Q(in) * A(in) * Z(in)**H = Q(out) * A(out) * Z(out)**H
+*> Q(in) * B(in) * Z(in)**H = Q(out) * B(out) * Z(out)**H
+*>
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] WANTQ
+*> \verbatim
+*> WANTQ is LOGICAL
+*> .TRUE. : update the left transformation matrix Q;
+*> .FALSE.: do not update Q.
+*> \endverbatim
+*>
+*> \param[in] WANTZ
+*> \verbatim
+*> WANTZ is LOGICAL
+*> .TRUE. : update the right transformation matrix Z;
+*> .FALSE.: do not update Z.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrices A and B. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimensions (LDA,N)
+*> On entry, the matrix A in the pair (A, B).
+*> On exit, the updated matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*> B is COMPLEX*16 array, dimensions (LDB,N)
+*> On entry, the matrix B in the pair (A, B).
+*> On exit, the updated matrix B.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of the array B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] Q
+*> \verbatim
+*> Q is COMPLEX*16 array, dimension (LDQ,N)
+*> If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
+*> the updated matrix Q.
+*> Not referenced if WANTQ = .FALSE..
+*> \endverbatim
+*>
+*> \param[in] LDQ
+*> \verbatim
+*> LDQ is INTEGER
+*> The leading dimension of the array Q. LDQ >= 1;
+*> If WANTQ = .TRUE., LDQ >= N.
+*> \endverbatim
+*>
+*> \param[in,out] Z
+*> \verbatim
+*> Z is COMPLEX*16 array, dimension (LDZ,N)
+*> If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
+*> the updated matrix Z.
+*> Not referenced if WANTZ = .FALSE..
+*> \endverbatim
+*>
+*> \param[in] LDZ
+*> \verbatim
+*> LDZ is INTEGER
+*> The leading dimension of the array Z. LDZ >= 1;
+*> If WANTZ = .TRUE., LDZ >= N.
+*> \endverbatim
+*>
+*> \param[in] J1
+*> \verbatim
+*> J1 is INTEGER
+*> The index to the first block (A11, B11).
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> =0: Successful exit.
+*> =1: The transformed matrix pair (A, B) would be too far
+*> from generalized Schur form; the problem is ill-
+*> conditioned.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup complex16GEauxiliary
+*
+*> \par Further Details:
+* =====================
+*>
+*> In the current code both weak and strong stability tests are
+*> performed. The user can omit the strong stability test by changing
+*> the internal logical parameter WANDS to .FALSE.. See ref. [2] for
+*> details.
+*
+*> \par Contributors:
+* ==================
+*>
+*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
+*> Umea University, S-901 87 Umea, Sweden.
+*
+*> \par References:
+* ================
+*>
+*> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
+*> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
+*> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
+*> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
+*> \n
+*> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
+*> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
+*> Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
+*> Department of Computing Science, Umea University, S-901 87 Umea,
+*> Sweden, 1994. Also as LAPACK Working Note 87. To appear in
+*> Numerical Algorithms, 1996.
+*>
+* =====================================================================
SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
$ LDZ, J1, INFO )
*
-* -- LAPACK auxiliary routine (version 3.2) --
+* -- LAPACK auxiliary routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2006
*
* .. Scalar Arguments ..
LOGICAL WANTQ, WANTZ
@@ -15,121 +201,24 @@
$ Z( LDZ, * )
* ..
*
-* Purpose
-* =======
-*
-* ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
-* in an upper triangular matrix pair (A, B) by an unitary equivalence
-* transformation.
-*
-* (A, B) must be in generalized Schur canonical form, that is, A and
-* B are both upper triangular.
-*
-* Optionally, the matrices Q and Z of generalized Schur vectors are
-* updated.
-*
-* Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
-* Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
-*
-*
-* Arguments
-* =========
-*
-* WANTQ (input) LOGICAL
-* .TRUE. : update the left transformation matrix Q;
-* .FALSE.: do not update Q.
-*
-* WANTZ (input) LOGICAL
-* .TRUE. : update the right transformation matrix Z;
-* .FALSE.: do not update Z.
-*
-* N (input) INTEGER
-* The order of the matrices A and B. N >= 0.
-*
-* A (input/output) COMPLEX*16 arrays, dimensions (LDA,N)
-* On entry, the matrix A in the pair (A, B).
-* On exit, the updated matrix A.
-*
-* LDA (input) INTEGER
-* The leading dimension of the array A. LDA >= max(1,N).
-*
-* B (input/output) COMPLEX*16 arrays, dimensions (LDB,N)
-* On entry, the matrix B in the pair (A, B).
-* On exit, the updated matrix B.
-*
-* LDB (input) INTEGER
-* The leading dimension of the array B. LDB >= max(1,N).
-*
-* Q (input/output) COMPLEX*16 array, dimension (LDZ,N)
-* If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit,
-* the updated matrix Q.
-* Not referenced if WANTQ = .FALSE..
-*
-* LDQ (input) INTEGER
-* The leading dimension of the array Q. LDQ >= 1;
-* If WANTQ = .TRUE., LDQ >= N.
-*
-* Z (input/output) COMPLEX*16 array, dimension (LDZ,N)
-* If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit,
-* the updated matrix Z.
-* Not referenced if WANTZ = .FALSE..
-*
-* LDZ (input) INTEGER
-* The leading dimension of the array Z. LDZ >= 1;
-* If WANTZ = .TRUE., LDZ >= N.
-*
-* J1 (input) INTEGER
-* The index to the first block (A11, B11).
-*
-* INFO (output) INTEGER
-* =0: Successful exit.
-* =1: The transformed matrix pair (A, B) would be too far
-* from generalized Schur form; the problem is ill-
-* conditioned.
-*
-*
-* Further Details
-* ===============
-*
-* Based on contributions by
-* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
-* Umea University, S-901 87 Umea, Sweden.
-*
-* In the current code both weak and strong stability tests are
-* performed. The user can omit the strong stability test by changing
-* the internal logical parameter WANDS to .FALSE.. See ref. [2] for
-* details.
-*
-* [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
-* Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
-* M.S. Moonen et al (eds), Linear Algebra for Large Scale and
-* Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
-*
-* [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
-* Eigenvalues of a Regular Matrix Pair (A, B) and Condition
-* Estimation: Theory, Algorithms and Software, Report UMINF-94.04,
-* Department of Computing Science, Umea University, S-901 87 Umea,
-* Sweden, 1994. Also as LAPACK Working Note 87. To appear in
-* Numerical Algorithms, 1996.
-*
* =====================================================================
*
* .. Parameters ..
COMPLEX*16 CZERO, CONE
PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
$ CONE = ( 1.0D+0, 0.0D+0 ) )
- DOUBLE PRECISION TEN
- PARAMETER ( TEN = 10.0D+0 )
+ DOUBLE PRECISION TWENTY
+ PARAMETER ( TWENTY = 2.0D+1 )
INTEGER LDST
PARAMETER ( LDST = 2 )
LOGICAL WANDS
PARAMETER ( WANDS = .TRUE. )
* ..
* .. Local Scalars ..
- LOGICAL DTRONG, WEAK
+ LOGICAL STRONG, WEAK
INTEGER I, M
- DOUBLE PRECISION CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SS, SUM,
- $ THRESH, WS
+ DOUBLE PRECISION CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SUM,
+ $ THRESHA, THRESHB
COMPLEX*16 CDUM, F, G, SQ, SZ
* ..
* .. Local Arrays ..
@@ -156,7 +245,7 @@
*
M = LDST
WEAK = .FALSE.
- DTRONG = .FALSE.
+ STRONG = .FALSE.
*
* Make a local copy of selected block in (A, B)
*
@@ -171,17 +260,31 @@
SUM = DBLE( CONE )
CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
- CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
+ CALL ZLASSQ( M*M, WORK, 1, SCALE, SUM )
SA = SCALE*SQRT( SUM )
- THRESH = MAX( TEN*EPS*SA, SMLNUM )
+ SCALE = DBLE( CZERO )
+ SUM = DBLE( CONE )
+ CALL ZLASSQ( M*M, WORK(M*M+1), 1, SCALE, SUM )
+ SB = SCALE*SQRT( SUM )
+*
+* THRES has been changed from
+* THRESH = MAX( TEN*EPS*SA, SMLNUM )
+* to
+* THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
+* on 04/01/10.
+* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
+* Jim Demmel and Guillaume Revy. See forum post 1783.
+*
+ THRESHA = MAX( TWENTY*EPS*SA, SMLNUM )
+ THRESHB = MAX( TWENTY*EPS*SB, SMLNUM )
*
* Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
* using Givens rotations and perform the swap tentatively.
*
F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 )
G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 )
- SA = ABS( S( 2, 2 ) )
- SB = ABS( T( 2, 2 ) )
+ SA = ABS( S( 2, 2 ) ) * ABS( T( 1, 1 ) )
+ SB = ABS( S( 1, 1 ) ) * ABS( T( 2, 2 ) )
CALL ZLARTG( G, F, CZ, SZ, CDUM )
SZ = -SZ
CALL ZROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, CZ, DCONJG( SZ ) )
@@ -194,17 +297,20 @@
CALL ZROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, CQ, SQ )
CALL ZROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, CQ, SQ )
*
-* Weak stability test: |S21| + |T21| <= O(EPS F-norm((S, T)))
+* Weak stability test: |S21| <= O(EPS F-norm((A)))
+* and |T21| <= O(EPS F-norm((B)))
*
- WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) )
- WEAK = WS.LE.THRESH
+ WEAK = ABS( S( 2, 1 ) ).LE.THRESHA .AND.
+ $ ABS( T( 2, 1 ) ).LE.THRESHB
IF( .NOT.WEAK )
$ GO TO 20
*
IF( WANDS ) THEN
*
* Strong stability test:
-* F-norm((A-QL'*S*QR, B-QL'*T*QR)) <= O(EPS*F-norm((A, B)))
+* F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
+* and
+* F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
*
CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M )
CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M )
@@ -220,10 +326,14 @@
10 CONTINUE
SCALE = DBLE( CZERO )
SUM = DBLE( CONE )
- CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM )
- SS = SCALE*SQRT( SUM )
- DTRONG = SS.LE.THRESH
- IF( .NOT.DTRONG )
+ CALL ZLASSQ( M*M, WORK, 1, SCALE, SUM )
+ SA = SCALE*SQRT( SUM )
+ SCALE = DBLE( CZERO )
+ SUM = DBLE( CONE )
+ CALL ZLASSQ( M*M, WORK(M*M+1), 1, SCALE, SUM )
+ SB = SCALE*SQRT( SUM )
+ STRONG = SA.LE.THRESHA .AND. SB.LE.THRESHB
+ IF( .NOT.STRONG )
$ GO TO 20
END IF
*