version 1.20, 2018/05/29 07:18:10
|
version 1.21, 2023/08/07 08:39:12
|
Line 181
|
Line 181
|
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date December 2016 |
|
* |
|
*> \ingroup doubleGEauxiliary |
*> \ingroup doubleGEauxiliary |
* |
* |
*> \par Further Details: |
*> \par Further Details: |
Line 221
|
Line 219
|
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.7.0) -- |
* -- 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..-- |
* December 2016 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
LOGICAL WANTQ, WANTZ |
LOGICAL WANTQ, WANTZ |
Line 250
|
Line 247
|
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 293
|
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 310
|
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 ) |
* |
* |
* THRES has been changed from |
* THRES has been changed from |
* THRESH = MAX( TEN*EPS*SA, SMLNUM ) |
* THRESH = MAX( TEN*EPS*SA, SMLNUM ) |
Line 322
|
Line 323
|
* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by |
* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by |
* Jim Demmel and Guillaume Revy. See forum post 1783. |
* Jim Demmel and Guillaume Revy. See forum post 1783. |
* |
* |
THRESH = MAX( TWENTY*EPS*DNORM, SMLNUM ) |
THRESHA = MAX( TWENTY*EPS*DNORMA, SMLNUM ) |
|
THRESHB = MAX( TWENTY*EPS*DNORMB, SMLNUM ) |
* |
* |
IF( M.EQ.2 ) THEN |
IF( M.EQ.2 ) THEN |
* |
* |
Line 333
|
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 356
|
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**T*S*QR, B-QL**T*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 378
|
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 385
|
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 439
|
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: |
* |
* |
Line 538
|
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 556
|
Line 565
|
IF( WANDS ) THEN |
IF( WANDS ) THEN |
* |
* |
* Strong stability test: |
* Strong stability test: |
* F-norm((A-QL*S*QR**T, B-QL*T*QR**T)) <= 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 567
|
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 574
|
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 |