version 1.1.1.1, 2010/01/26 15:22:45
|
version 1.9, 2011/07/22 07:38:12
|
Line 1
|
Line 1
|
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 (version 3.3.1) -- |
* -- 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 |
* -- April 2011 -- |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
LOGICAL WANTQ, WANTZ |
LOGICAL WANTQ, WANTZ |
Line 29
|
Line 29
|
* Optionally, the matrices Q and Z of generalized Schur vectors are |
* Optionally, the matrices Q and Z of generalized Schur vectors are |
* updated. |
* updated. |
* |
* |
* Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)' |
* Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T |
* Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)' |
* Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T |
* |
* |
* |
* |
* Arguments |
* Arguments |
Line 47
|
Line 47
|
* N (input) INTEGER |
* N (input) INTEGER |
* The order of the matrices A and B. N >= 0. |
* The order of the matrices A and B. N >= 0. |
* |
* |
* A (input/output) DOUBLE PRECISION arrays, dimensions (LDA,N) |
* A (input/output) DOUBLE PRECISION array, dimensions (LDA,N) |
* On entry, the matrix A in the pair (A, B). |
* On entry, the matrix A in the pair (A, B). |
* On exit, the updated matrix A. |
* On exit, the updated matrix A. |
* |
* |
* LDA (input) INTEGER |
* LDA (input) INTEGER |
* The leading dimension of the array A. LDA >= max(1,N). |
* The leading dimension of the array A. LDA >= max(1,N). |
* |
* |
* B (input/output) DOUBLE PRECISION arrays, dimensions (LDB,N) |
* B (input/output) DOUBLE PRECISION array, dimensions (LDB,N) |
* On entry, the matrix B in the pair (A, B). |
* On entry, the matrix B in the pair (A, B). |
* On exit, the updated matrix B. |
* On exit, the updated matrix B. |
* |
* |
* LDB (input) INTEGER |
* LDB (input) INTEGER |
* The leading dimension of the array B. LDB >= max(1,N). |
* The leading dimension of the array B. LDB >= max(1,N). |
* |
* |
* Q (input/output) DOUBLE PRECISION array, dimension (LDZ,N) |
* Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) |
* On entry, if WANTQ = .TRUE., the orthogonal matrix Q. |
* On entry, if WANTQ = .TRUE., the orthogonal matrix Q. |
* On exit, the updated matrix Q. |
* On exit, the updated matrix Q. |
* Not referenced if WANTQ = .FALSE.. |
* Not referenced if WANTQ = .FALSE.. |
Line 134
|
Line 134
|
* .. 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 |
Line 205
|
Line 205
|
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 ) |
DNORM = 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. |
|
* |
|
THRESH = MAX( TWENTY*EPS*DNORM, SMLNUM ) |
* |
* |
IF( M.EQ.2 ) THEN |
IF( M.EQ.2 ) THEN |
* |
* |
Line 250
|
Line 259
|
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**T*S*QR, B-QL**T*T*QR)) <= O(EPS*F-norm((A,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 325
|
Line 334
|
* |
* |
* 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 353
|
* |
* |
* 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 439
|
Line 448
|
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*S*QR**T, B-QL*T*QR**T)) <= O(EPS*F-norm((A,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 ) |