version 1.4, 2010/08/06 15:32:36
|
version 1.21, 2023/08/07 08:39:12
|
Line 1
|
Line 1
|
|
*> \brief \b DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equivalence transformation. |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download DTGEX2 + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgex2.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgex2.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgex2.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, |
|
* LDZ, J1, N1, N2, WORK, LWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* LOGICAL WANTQ, WANTZ |
|
* INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, LWORK, N, N1, N2 |
|
* .. |
|
* .. Array Arguments .. |
|
* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ), |
|
* $ WORK( * ), Z( LDZ, * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DTGEX2 swaps adjacent diagonal blocks (A11, B11) and (A22, B22) |
|
*> of size 1-by-1 or 2-by-2 in an upper (quasi) triangular matrix pair |
|
*> (A, B) by an orthogonal equivalence transformation. |
|
*> |
|
*> (A, B) must be in generalized real Schur canonical form (as returned |
|
*> by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2 |
|
*> diagonal blocks. B is upper triangular. |
|
*> |
|
*> Optionally, the matrices Q and Z of generalized Schur vectors are |
|
*> updated. |
|
*> |
|
*> Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T |
|
*> Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T |
|
*> |
|
*> \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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDQ,N) |
|
*> On entry, if WANTQ = .TRUE., the orthogonal 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 DOUBLE PRECISION array, dimension (LDZ,N) |
|
*> On entry, if WANTZ =.TRUE., the orthogonal 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). 1 <= J1 <= N. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] N1 |
|
*> \verbatim |
|
*> N1 is INTEGER |
|
*> The order of the first block (A11, B11). N1 = 0, 1 or 2. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] N2 |
|
*> \verbatim |
|
*> N2 is INTEGER |
|
*> The order of the second block (A22, B22). N2 = 0, 1 or 2. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)). |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LWORK |
|
*> \verbatim |
|
*> LWORK is INTEGER |
|
*> The dimension of the array WORK. |
|
*> LWORK >= MAX( 1, N*(N2+N1), (N2+N1)*(N2+N1)*2 ) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> =0: Successful exit |
|
*> >0: If INFO = 1, the transformed matrix (A, B) would be |
|
*> too far from generalized Schur form; the blocks are |
|
*> not swapped and (A, B) and (Q, Z) are unchanged. |
|
*> The problem of swapping is too ill-conditioned. |
|
*> <0: If INFO = -16: LWORK is too small. Appropriate value |
|
*> for LWORK is returned in WORK(1). |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \ingroup doubleGEauxiliary |
|
* |
|
*> \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: |
|
* ================ |
|
*> |
|
*> \verbatim |
|
*> |
|
*> [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. |
|
*> \endverbatim |
|
*> |
|
* ===================================================================== |
SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, |
SUBROUTINE DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, |
$ LDZ, J1, N1, N2, WORK, LWORK, INFO ) |
$ LDZ, J1, N1, N2, WORK, LWORK, INFO ) |
* |
* |
* -- LAPACK auxiliary routine (version 3.2) -- |
* -- LAPACK auxiliary routine -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2006 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
LOGICAL WANTQ, WANTZ |
LOGICAL WANTQ, WANTZ |
Line 15
|
Line 232
|
$ WORK( * ), Z( LDZ, * ) |
$ WORK( * ), Z( LDZ, * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* DTGEX2 swaps adjacent diagonal blocks (A11, B11) and (A22, B22) |
|
* of size 1-by-1 or 2-by-2 in an upper (quasi) triangular matrix pair |
|
* (A, B) by an orthogonal equivalence transformation. |
|
* |
|
* (A, B) must be in generalized real Schur canonical form (as returned |
|
* by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2 |
|
* diagonal blocks. B is 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDZ,N) |
|
* On entry, if WANTQ = .TRUE., the orthogonal 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) DOUBLE PRECISION array, dimension (LDZ,N) |
|
* On entry, if WANTZ =.TRUE., the orthogonal 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). 1 <= J1 <= N. |
|
* |
|
* N1 (input) INTEGER |
|
* The order of the first block (A11, B11). N1 = 0, 1 or 2. |
|
* |
|
* N2 (input) INTEGER |
|
* The order of the second block (A22, B22). N2 = 0, 1 or 2. |
|
* |
|
* WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)). |
|
* |
|
* LWORK (input) INTEGER |
|
* The dimension of the array WORK. |
|
* LWORK >= MAX( 1, N*(N2+N1), (N2+N1)*(N2+N1)*2 ) |
|
* |
|
* INFO (output) INTEGER |
|
* =0: Successful exit |
|
* >0: If INFO = 1, the transformed matrix (A, B) would be |
|
* too far from generalized Schur form; the blocks are |
|
* not swapped and (A, B) and (Q, Z) are unchanged. |
|
* The problem of swapping is too ill-conditioned. |
|
* <0: If INFO = -16: LWORK is too small. Appropriate value |
|
* for LWORK is returned in WORK(1). |
|
* |
|
* 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. |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* Replaced various illegal calls to DCOPY by calls to DLASET, or by DO |
* Replaced various illegal calls to DCOPY by calls to DLASET, or by DO |
* loops. Sven Hammarling, 1/5/02. |
* loops. Sven Hammarling, 1/5/02. |
Line 134
|
Line 239
|
* .. Parameters .. |
* .. Parameters .. |
DOUBLE PRECISION ZERO, ONE |
DOUBLE PRECISION ZERO, ONE |
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
DOUBLE PRECISION TEN |
DOUBLE PRECISION TWENTY |
PARAMETER ( TEN = 1.0D+01 ) |
PARAMETER ( TWENTY = 2.0D+01 ) |
INTEGER LDST |
INTEGER LDST |
PARAMETER ( LDST = 4 ) |
PARAMETER ( LDST = 4 ) |
LOGICAL WANDS |
LOGICAL WANDS |
PARAMETER ( WANDS = .TRUE. ) |
PARAMETER ( WANDS = .TRUE. ) |
* .. |
* .. |
* .. Local Scalars .. |
* .. Local Scalars .. |
LOGICAL DTRONG, WEAK |
LOGICAL STRONG, WEAK |
INTEGER I, IDUM, LINFO, M |
INTEGER I, IDUM, LINFO, M |
DOUBLE PRECISION BQRA21, BRQA21, DDUM, DNORM, DSCALE, DSUM, EPS, |
DOUBLE PRECISION BQRA21, BRQA21, DDUM, DNORMA, DNORMB, DSCALE, |
$ F, G, SA, SB, SCALE, SMLNUM, SS, THRESH, WS |
$ DSUM, EPS, F, G, SA, SB, SCALE, SMLNUM, |
|
$ THRESHA, THRESHB |
* .. |
* .. |
* .. Local Arrays .. |
* .. Local Arrays .. |
INTEGER IWORK( LDST ) |
INTEGER IWORK( LDST ) |
Line 185
|
Line 291
|
END IF |
END IF |
* |
* |
WEAK = .FALSE. |
WEAK = .FALSE. |
DTRONG = .FALSE. |
STRONG = .FALSE. |
* |
* |
* Make a local copy of selected block |
* Make a local copy of selected block |
* |
* |
Line 202
|
Line 308
|
DSUM = ONE |
DSUM = ONE |
CALL DLACPY( 'Full', M, M, S, LDST, WORK, M ) |
CALL DLACPY( 'Full', M, M, S, LDST, WORK, M ) |
CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM ) |
CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM ) |
|
DNORMA = DSCALE*SQRT( DSUM ) |
|
DSCALE = ZERO |
|
DSUM = ONE |
CALL DLACPY( 'Full', M, M, T, LDST, WORK, M ) |
CALL DLACPY( 'Full', M, M, T, LDST, WORK, M ) |
CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM ) |
CALL DLASSQ( M*M, WORK, 1, DSCALE, DSUM ) |
DNORM = DSCALE*SQRT( DSUM ) |
DNORMB = DSCALE*SQRT( DSUM ) |
THRESH = MAX( TEN*EPS*DNORM, SMLNUM ) |
* |
|
* 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*DNORMA, SMLNUM ) |
|
THRESHB = MAX( TWENTY*EPS*DNORMB, SMLNUM ) |
* |
* |
IF( M.EQ.2 ) THEN |
IF( M.EQ.2 ) THEN |
* |
* |
Line 216
|
Line 335
|
* |
* |
F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 ) |
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 ) |
G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 ) |
SB = ABS( T( 2, 2 ) ) |
SA = ABS( S( 2, 2 ) ) * ABS( T( 1, 1 ) ) |
SA = ABS( S( 2, 2 ) ) |
SB = ABS( S( 1, 1 ) ) * ABS( T( 2, 2 ) ) |
CALL DLARTG( F, G, IR( 1, 2 ), IR( 1, 1 ), DDUM ) |
CALL DLARTG( F, G, IR( 1, 2 ), IR( 1, 1 ), DDUM ) |
IR( 2, 1 ) = -IR( 1, 2 ) |
IR( 2, 1 ) = -IR( 1, 2 ) |
IR( 2, 2 ) = IR( 1, 1 ) |
IR( 2, 2 ) = IR( 1, 1 ) |
Line 239
|
Line 358
|
LI( 2, 2 ) = LI( 1, 1 ) |
LI( 2, 2 ) = LI( 1, 1 ) |
LI( 1, 2 ) = -LI( 2, 1 ) |
LI( 1, 2 ) = -LI( 2, 1 ) |
* |
* |
* Weak stability test: |
* Weak stability test: |S21| <= O(EPS F-norm((A))) |
* |S21| + |T21| <= O(EPS * F-norm((S, T))) |
* and |T21| <= O(EPS F-norm((B))) |
* |
* |
WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) ) |
WEAK = ABS( S( 2, 1 ) ) .LE. THRESHA .AND. |
WEAK = WS.LE.THRESH |
$ ABS( T( 2, 1 ) ) .LE. THRESHB |
IF( .NOT.WEAK ) |
IF( .NOT.WEAK ) |
$ GO TO 70 |
$ GO TO 70 |
* |
* |
IF( WANDS ) THEN |
IF( WANDS ) THEN |
* |
* |
* Strong stability test: |
* 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 DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ), |
CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ), |
$ M ) |
$ M ) |
Line 261
|
Line 382
|
DSCALE = ZERO |
DSCALE = ZERO |
DSUM = ONE |
DSUM = ONE |
CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) |
CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) |
|
SA = DSCALE*SQRT( DSUM ) |
* |
* |
CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ), |
CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ), |
$ M ) |
$ M ) |
Line 268
|
Line 390
|
$ WORK, M ) |
$ WORK, M ) |
CALL DGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE, |
CALL DGEMM( 'N', 'T', M, M, M, -ONE, WORK, M, IR, LDST, ONE, |
$ WORK( M*M+1 ), M ) |
$ WORK( M*M+1 ), M ) |
|
DSCALE = ZERO |
|
DSUM = ONE |
CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) |
CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) |
SS = DSCALE*SQRT( DSUM ) |
SB = DSCALE*SQRT( DSUM ) |
DTRONG = SS.LE.THRESH |
STRONG = SA.LE.THRESHA .AND. SB.LE.THRESHB |
IF( .NOT.DTRONG ) |
IF( .NOT.STRONG ) |
$ GO TO 70 |
$ GO TO 70 |
END IF |
END IF |
* |
* |
Line 322
|
Line 446
|
$ IR( N2+1, N1+1 ), LDST, T, LDST, T( N1+1, N1+1 ), |
$ IR( N2+1, N1+1 ), LDST, T, LDST, T( N1+1, N1+1 ), |
$ LDST, LI, LDST, SCALE, DSUM, DSCALE, IWORK, IDUM, |
$ LDST, LI, LDST, SCALE, DSUM, DSCALE, IWORK, IDUM, |
$ LINFO ) |
$ LINFO ) |
|
IF( LINFO.NE.0 ) |
|
$ GO TO 70 |
* |
* |
* Compute orthogonal matrix QL: |
* Compute orthogonal matrix QL: |
* |
* |
* QL' * LI = [ TL ] |
* QL**T * LI = [ TL ] |
* [ 0 ] |
* [ 0 ] |
* where |
* where |
* LI = [ -L ] |
* LI = [ -L ] |
* [ SCALE * identity(N2) ] |
* [ SCALE * identity(N2) ] |
Line 344
|
Line 470
|
* |
* |
* Compute orthogonal matrix RQ: |
* Compute orthogonal matrix RQ: |
* |
* |
* IR * RQ' = [ 0 TR], |
* IR * RQ**T = [ 0 TR], |
* |
* |
* where IR = [ SCALE * identity(N1), R ] |
* where IR = [ SCALE * identity(N1), R ] |
* |
* |
Line 421
|
Line 547
|
* |
* |
* Decide which method to use. |
* Decide which method to use. |
* Weak stability test: |
* Weak stability test: |
* F-norm(S21) <= O(EPS * F-norm((S, T))) |
* F-norm(S21) <= O(EPS * F-norm((S))) |
* |
* |
IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESH ) THEN |
IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESHA ) THEN |
CALL DLACPY( 'F', M, M, SCPY, LDST, S, LDST ) |
CALL DLACPY( 'F', M, M, SCPY, LDST, S, LDST ) |
CALL DLACPY( 'F', M, M, TCPY, LDST, T, LDST ) |
CALL DLACPY( 'F', M, M, TCPY, LDST, T, LDST ) |
CALL DLACPY( 'F', M, M, IRCOP, LDST, IR, LDST ) |
CALL DLACPY( 'F', M, M, IRCOP, LDST, IR, LDST ) |
CALL DLACPY( 'F', M, M, LICOP, LDST, LI, LDST ) |
CALL DLACPY( 'F', M, M, LICOP, LDST, LI, LDST ) |
ELSE IF( BRQA21.GE.THRESH ) THEN |
ELSE IF( BRQA21.GE.THRESHA ) THEN |
GO TO 70 |
GO TO 70 |
END IF |
END IF |
* |
* |
Line 439
|
Line 565
|
IF( WANDS ) THEN |
IF( WANDS ) THEN |
* |
* |
* Strong stability test: |
* 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 DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ), |
CALL DLACPY( 'Full', M, M, A( J1, J1 ), LDA, WORK( M*M+1 ), |
$ M ) |
$ M ) |
Line 450
|
Line 578
|
DSCALE = ZERO |
DSCALE = ZERO |
DSUM = ONE |
DSUM = ONE |
CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) |
CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) |
|
SA = DSCALE*SQRT( DSUM ) |
* |
* |
CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ), |
CALL DLACPY( 'Full', M, M, B( J1, J1 ), LDB, WORK( M*M+1 ), |
$ M ) |
$ M ) |
Line 457
|
Line 586
|
$ WORK, M ) |
$ WORK, M ) |
CALL DGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE, |
CALL DGEMM( 'N', 'N', M, M, M, -ONE, WORK, M, IR, LDST, ONE, |
$ WORK( M*M+1 ), M ) |
$ WORK( M*M+1 ), M ) |
|
DSCALE = ZERO |
|
DSUM = ONE |
CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) |
CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) |
SS = DSCALE*SQRT( DSUM ) |
SB = DSCALE*SQRT( DSUM ) |
DTRONG = ( SS.LE.THRESH ) |
STRONG = SA.LE.THRESHA .AND. SB.LE.THRESHB |
IF( .NOT.DTRONG ) |
IF( .NOT.STRONG ) |
$ GO TO 70 |
$ GO TO 70 |
* |
* |
END IF |
END IF |
Line 478
|
Line 609
|
* |
* |
* Standardize existing 2-by-2 blocks. |
* Standardize existing 2-by-2 blocks. |
* |
* |
DO 50 I = 1, M*M |
CALL DLASET( 'Full', M, M, ZERO, ZERO, WORK, M ) |
WORK(I) = ZERO |
|
50 CONTINUE |
|
WORK( 1 ) = ONE |
WORK( 1 ) = ONE |
T( 1, 1 ) = ONE |
T( 1, 1 ) = ONE |
IDUM = LWORK - M*M - 2 |
IDUM = LWORK - M*M - 2 |
Line 551
|
Line 680
|
$ A( J1, I ), LDA, ZERO, WORK, M ) |
$ A( J1, I ), LDA, ZERO, WORK, M ) |
CALL DLACPY( 'Full', M, N-I+1, WORK, M, A( J1, I ), LDA ) |
CALL DLACPY( 'Full', M, N-I+1, WORK, M, A( J1, I ), LDA ) |
CALL DGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST, |
CALL DGEMM( 'T', 'N', M, N-I+1, M, ONE, LI, LDST, |
$ B( J1, I ), LDA, ZERO, WORK, M ) |
$ B( J1, I ), LDB, ZERO, WORK, M ) |
CALL DLACPY( 'Full', M, N-I+1, WORK, M, B( J1, I ), LDB ) |
CALL DLACPY( 'Full', M, N-I+1, WORK, M, B( J1, I ), LDB ) |
END IF |
END IF |
I = J1 - 1 |
I = J1 - 1 |