Annotation of rpl/lapack/lapack/dtgsen.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
! 2: $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL,
! 3: $ PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO )
! 4: *
! 5: * -- LAPACK routine (version 3.2) --
! 6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 8: * January 2007
! 9: *
! 10: * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
! 11: *
! 12: * .. Scalar Arguments ..
! 13: LOGICAL WANTQ, WANTZ
! 14: INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
! 15: $ M, N
! 16: DOUBLE PRECISION PL, PR
! 17: * ..
! 18: * .. Array Arguments ..
! 19: LOGICAL SELECT( * )
! 20: INTEGER IWORK( * )
! 21: DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
! 22: $ B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ),
! 23: $ WORK( * ), Z( LDZ, * )
! 24: * ..
! 25: *
! 26: * Purpose
! 27: * =======
! 28: *
! 29: * DTGSEN reorders the generalized real Schur decomposition of a real
! 30: * matrix pair (A, B) (in terms of an orthonormal equivalence trans-
! 31: * formation Q' * (A, B) * Z), so that a selected cluster of eigenvalues
! 32: * appears in the leading diagonal blocks of the upper quasi-triangular
! 33: * matrix A and the upper triangular B. The leading columns of Q and
! 34: * Z form orthonormal bases of the corresponding left and right eigen-
! 35: * spaces (deflating subspaces). (A, B) must be in generalized real
! 36: * Schur canonical form (as returned by DGGES), i.e. A is block upper
! 37: * triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper
! 38: * triangular.
! 39: *
! 40: * DTGSEN also computes the generalized eigenvalues
! 41: *
! 42: * w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)
! 43: *
! 44: * of the reordered matrix pair (A, B).
! 45: *
! 46: * Optionally, DTGSEN computes the estimates of reciprocal condition
! 47: * numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
! 48: * (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
! 49: * between the matrix pairs (A11, B11) and (A22,B22) that correspond to
! 50: * the selected cluster and the eigenvalues outside the cluster, resp.,
! 51: * and norms of "projections" onto left and right eigenspaces w.r.t.
! 52: * the selected cluster in the (1,1)-block.
! 53: *
! 54: * Arguments
! 55: * =========
! 56: *
! 57: * IJOB (input) INTEGER
! 58: * Specifies whether condition numbers are required for the
! 59: * cluster of eigenvalues (PL and PR) or the deflating subspaces
! 60: * (Difu and Difl):
! 61: * =0: Only reorder w.r.t. SELECT. No extras.
! 62: * =1: Reciprocal of norms of "projections" onto left and right
! 63: * eigenspaces w.r.t. the selected cluster (PL and PR).
! 64: * =2: Upper bounds on Difu and Difl. F-norm-based estimate
! 65: * (DIF(1:2)).
! 66: * =3: Estimate of Difu and Difl. 1-norm-based estimate
! 67: * (DIF(1:2)).
! 68: * About 5 times as expensive as IJOB = 2.
! 69: * =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic
! 70: * version to get it all.
! 71: * =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above)
! 72: *
! 73: * WANTQ (input) LOGICAL
! 74: * .TRUE. : update the left transformation matrix Q;
! 75: * .FALSE.: do not update Q.
! 76: *
! 77: * WANTZ (input) LOGICAL
! 78: * .TRUE. : update the right transformation matrix Z;
! 79: * .FALSE.: do not update Z.
! 80: *
! 81: * SELECT (input) LOGICAL array, dimension (N)
! 82: * SELECT specifies the eigenvalues in the selected cluster.
! 83: * To select a real eigenvalue w(j), SELECT(j) must be set to
! 84: * .TRUE.. To select a complex conjugate pair of eigenvalues
! 85: * w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
! 86: * either SELECT(j) or SELECT(j+1) or both must be set to
! 87: * .TRUE.; a complex conjugate pair of eigenvalues must be
! 88: * either both included in the cluster or both excluded.
! 89: *
! 90: * N (input) INTEGER
! 91: * The order of the matrices A and B. N >= 0.
! 92: *
! 93: * A (input/output) DOUBLE PRECISION array, dimension(LDA,N)
! 94: * On entry, the upper quasi-triangular matrix A, with (A, B) in
! 95: * generalized real Schur canonical form.
! 96: * On exit, A is overwritten by the reordered matrix A.
! 97: *
! 98: * LDA (input) INTEGER
! 99: * The leading dimension of the array A. LDA >= max(1,N).
! 100: *
! 101: * B (input/output) DOUBLE PRECISION array, dimension(LDB,N)
! 102: * On entry, the upper triangular matrix B, with (A, B) in
! 103: * generalized real Schur canonical form.
! 104: * On exit, B is overwritten by the reordered matrix B.
! 105: *
! 106: * LDB (input) INTEGER
! 107: * The leading dimension of the array B. LDB >= max(1,N).
! 108: *
! 109: * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
! 110: * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
! 111: * BETA (output) DOUBLE PRECISION array, dimension (N)
! 112: * On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
! 113: * be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i
! 114: * and BETA(j),j=1,...,N are the diagonals of the complex Schur
! 115: * form (S,T) that would result if the 2-by-2 diagonal blocks of
! 116: * the real generalized Schur form of (A,B) were further reduced
! 117: * to triangular form using complex unitary transformations.
! 118: * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
! 119: * positive, then the j-th and (j+1)-st eigenvalues are a
! 120: * complex conjugate pair, with ALPHAI(j+1) negative.
! 121: *
! 122: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
! 123: * On entry, if WANTQ = .TRUE., Q is an N-by-N matrix.
! 124: * On exit, Q has been postmultiplied by the left orthogonal
! 125: * transformation matrix which reorder (A, B); The leading M
! 126: * columns of Q form orthonormal bases for the specified pair of
! 127: * left eigenspaces (deflating subspaces).
! 128: * If WANTQ = .FALSE., Q is not referenced.
! 129: *
! 130: * LDQ (input) INTEGER
! 131: * The leading dimension of the array Q. LDQ >= 1;
! 132: * and if WANTQ = .TRUE., LDQ >= N.
! 133: *
! 134: * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
! 135: * On entry, if WANTZ = .TRUE., Z is an N-by-N matrix.
! 136: * On exit, Z has been postmultiplied by the left orthogonal
! 137: * transformation matrix which reorder (A, B); The leading M
! 138: * columns of Z form orthonormal bases for the specified pair of
! 139: * left eigenspaces (deflating subspaces).
! 140: * If WANTZ = .FALSE., Z is not referenced.
! 141: *
! 142: * LDZ (input) INTEGER
! 143: * The leading dimension of the array Z. LDZ >= 1;
! 144: * If WANTZ = .TRUE., LDZ >= N.
! 145: *
! 146: * M (output) INTEGER
! 147: * The dimension of the specified pair of left and right eigen-
! 148: * spaces (deflating subspaces). 0 <= M <= N.
! 149: *
! 150: * PL (output) DOUBLE PRECISION
! 151: * PR (output) DOUBLE PRECISION
! 152: * If IJOB = 1, 4 or 5, PL, PR are lower bounds on the
! 153: * reciprocal of the norm of "projections" onto left and right
! 154: * eigenspaces with respect to the selected cluster.
! 155: * 0 < PL, PR <= 1.
! 156: * If M = 0 or M = N, PL = PR = 1.
! 157: * If IJOB = 0, 2 or 3, PL and PR are not referenced.
! 158: *
! 159: * DIF (output) DOUBLE PRECISION array, dimension (2).
! 160: * If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
! 161: * If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
! 162: * Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based
! 163: * estimates of Difu and Difl.
! 164: * If M = 0 or N, DIF(1:2) = F-norm([A, B]).
! 165: * If IJOB = 0 or 1, DIF is not referenced.
! 166: *
! 167: * WORK (workspace/output) DOUBLE PRECISION array,
! 168: * dimension (MAX(1,LWORK))
! 169: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 170: *
! 171: * LWORK (input) INTEGER
! 172: * The dimension of the array WORK. LWORK >= 4*N+16.
! 173: * If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)).
! 174: * If IJOB = 3 or 5, LWORK >= MAX(4*N+16, 4*M*(N-M)).
! 175: *
! 176: * If LWORK = -1, then a workspace query is assumed; the routine
! 177: * only calculates the optimal size of the WORK array, returns
! 178: * this value as the first entry of the WORK array, and no error
! 179: * message related to LWORK is issued by XERBLA.
! 180: *
! 181: * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
! 182: * IF IJOB = 0, IWORK is not referenced. Otherwise,
! 183: * on exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
! 184: *
! 185: * LIWORK (input) INTEGER
! 186: * The dimension of the array IWORK. LIWORK >= 1.
! 187: * If IJOB = 1, 2 or 4, LIWORK >= N+6.
! 188: * If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-M), N+6).
! 189: *
! 190: * If LIWORK = -1, then a workspace query is assumed; the
! 191: * routine only calculates the optimal size of the IWORK array,
! 192: * returns this value as the first entry of the IWORK array, and
! 193: * no error message related to LIWORK is issued by XERBLA.
! 194: *
! 195: * INFO (output) INTEGER
! 196: * =0: Successful exit.
! 197: * <0: If INFO = -i, the i-th argument had an illegal value.
! 198: * =1: Reordering of (A, B) failed because the transformed
! 199: * matrix pair (A, B) would be too far from generalized
! 200: * Schur form; the problem is very ill-conditioned.
! 201: * (A, B) may have been partially reordered.
! 202: * If requested, 0 is returned in DIF(*), PL and PR.
! 203: *
! 204: * Further Details
! 205: * ===============
! 206: *
! 207: * DTGSEN first collects the selected eigenvalues by computing
! 208: * orthogonal U and W that move them to the top left corner of (A, B).
! 209: * In other words, the selected eigenvalues are the eigenvalues of
! 210: * (A11, B11) in:
! 211: *
! 212: * U'*(A, B)*W = (A11 A12) (B11 B12) n1
! 213: * ( 0 A22),( 0 B22) n2
! 214: * n1 n2 n1 n2
! 215: *
! 216: * where N = n1+n2 and U' means the transpose of U. The first n1 columns
! 217: * of U and W span the specified pair of left and right eigenspaces
! 218: * (deflating subspaces) of (A, B).
! 219: *
! 220: * If (A, B) has been obtained from the generalized real Schur
! 221: * decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then the
! 222: * reordered generalized real Schur form of (C, D) is given by
! 223: *
! 224: * (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',
! 225: *
! 226: * and the first n1 columns of Q*U and Z*W span the corresponding
! 227: * deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
! 228: *
! 229: * Note that if the selected eigenvalue is sufficiently ill-conditioned,
! 230: * then its value may differ significantly from its value before
! 231: * reordering.
! 232: *
! 233: * The reciprocal condition numbers of the left and right eigenspaces
! 234: * spanned by the first n1 columns of U and W (or Q*U and Z*W) may
! 235: * be returned in DIF(1:2), corresponding to Difu and Difl, resp.
! 236: *
! 237: * The Difu and Difl are defined as:
! 238: *
! 239: * Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
! 240: * and
! 241: * Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
! 242: *
! 243: * where sigma-min(Zu) is the smallest singular value of the
! 244: * (2*n1*n2)-by-(2*n1*n2) matrix
! 245: *
! 246: * Zu = [ kron(In2, A11) -kron(A22', In1) ]
! 247: * [ kron(In2, B11) -kron(B22', In1) ].
! 248: *
! 249: * Here, Inx is the identity matrix of size nx and A22' is the
! 250: * transpose of A22. kron(X, Y) is the Kronecker product between
! 251: * the matrices X and Y.
! 252: *
! 253: * When DIF(2) is small, small changes in (A, B) can cause large changes
! 254: * in the deflating subspace. An approximate (asymptotic) bound on the
! 255: * maximum angular error in the computed deflating subspaces is
! 256: *
! 257: * EPS * norm((A, B)) / DIF(2),
! 258: *
! 259: * where EPS is the machine precision.
! 260: *
! 261: * The reciprocal norm of the projectors on the left and right
! 262: * eigenspaces associated with (A11, B11) may be returned in PL and PR.
! 263: * They are computed as follows. First we compute L and R so that
! 264: * P*(A, B)*Q is block diagonal, where
! 265: *
! 266: * P = ( I -L ) n1 Q = ( I R ) n1
! 267: * ( 0 I ) n2 and ( 0 I ) n2
! 268: * n1 n2 n1 n2
! 269: *
! 270: * and (L, R) is the solution to the generalized Sylvester equation
! 271: *
! 272: * A11*R - L*A22 = -A12
! 273: * B11*R - L*B22 = -B12
! 274: *
! 275: * Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
! 276: * An approximate (asymptotic) bound on the average absolute error of
! 277: * the selected eigenvalues is
! 278: *
! 279: * EPS * norm((A, B)) / PL.
! 280: *
! 281: * There are also global error bounds which valid for perturbations up
! 282: * to a certain restriction: A lower bound (x) on the smallest
! 283: * F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
! 284: * coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
! 285: * (i.e. (A + E, B + F), is
! 286: *
! 287: * x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
! 288: *
! 289: * An approximate bound on x can be computed from DIF(1:2), PL and PR.
! 290: *
! 291: * If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
! 292: * (L', R') and unperturbed (L, R) left and right deflating subspaces
! 293: * associated with the selected cluster in the (1,1)-blocks can be
! 294: * bounded as
! 295: *
! 296: * max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
! 297: * max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
! 298: *
! 299: * See LAPACK User's Guide section 4.11 or the following references
! 300: * for more information.
! 301: *
! 302: * Note that if the default method for computing the Frobenius-norm-
! 303: * based estimate DIF is not wanted (see DLATDF), then the parameter
! 304: * IDIFJB (see below) should be changed from 3 to 4 (routine DLATDF
! 305: * (IJOB = 2 will be used)). See DTGSYL for more details.
! 306: *
! 307: * Based on contributions by
! 308: * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
! 309: * Umea University, S-901 87 Umea, Sweden.
! 310: *
! 311: * References
! 312: * ==========
! 313: *
! 314: * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
! 315: * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
! 316: * M.S. Moonen et al (eds), Linear Algebra for Large Scale and
! 317: * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
! 318: *
! 319: * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
! 320: * Eigenvalues of a Regular Matrix Pair (A, B) and Condition
! 321: * Estimation: Theory, Algorithms and Software,
! 322: * Report UMINF - 94.04, Department of Computing Science, Umea
! 323: * University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
! 324: * Note 87. To appear in Numerical Algorithms, 1996.
! 325: *
! 326: * [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
! 327: * for Solving the Generalized Sylvester Equation and Estimating the
! 328: * Separation between Regular Matrix Pairs, Report UMINF - 93.23,
! 329: * Department of Computing Science, Umea University, S-901 87 Umea,
! 330: * Sweden, December 1993, Revised April 1994, Also as LAPACK Working
! 331: * Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
! 332: * 1996.
! 333: *
! 334: * =====================================================================
! 335: *
! 336: * .. Parameters ..
! 337: INTEGER IDIFJB
! 338: PARAMETER ( IDIFJB = 3 )
! 339: DOUBLE PRECISION ZERO, ONE
! 340: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 341: * ..
! 342: * .. Local Scalars ..
! 343: LOGICAL LQUERY, PAIR, SWAP, WANTD, WANTD1, WANTD2,
! 344: $ WANTP
! 345: INTEGER I, IERR, IJB, K, KASE, KK, KS, LIWMIN, LWMIN,
! 346: $ MN2, N1, N2
! 347: DOUBLE PRECISION DSCALE, DSUM, EPS, RDSCAL, SMLNUM
! 348: * ..
! 349: * .. Local Arrays ..
! 350: INTEGER ISAVE( 3 )
! 351: * ..
! 352: * .. External Subroutines ..
! 353: EXTERNAL DLACN2, DLACPY, DLAG2, DLASSQ, DTGEXC, DTGSYL,
! 354: $ XERBLA
! 355: * ..
! 356: * .. External Functions ..
! 357: DOUBLE PRECISION DLAMCH
! 358: EXTERNAL DLAMCH
! 359: * ..
! 360: * .. Intrinsic Functions ..
! 361: INTRINSIC MAX, SIGN, SQRT
! 362: * ..
! 363: * .. Executable Statements ..
! 364: *
! 365: * Decode and test the input parameters
! 366: *
! 367: INFO = 0
! 368: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
! 369: *
! 370: IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN
! 371: INFO = -1
! 372: ELSE IF( N.LT.0 ) THEN
! 373: INFO = -5
! 374: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 375: INFO = -7
! 376: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 377: INFO = -9
! 378: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
! 379: INFO = -14
! 380: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
! 381: INFO = -16
! 382: END IF
! 383: *
! 384: IF( INFO.NE.0 ) THEN
! 385: CALL XERBLA( 'DTGSEN', -INFO )
! 386: RETURN
! 387: END IF
! 388: *
! 389: * Get machine constants
! 390: *
! 391: EPS = DLAMCH( 'P' )
! 392: SMLNUM = DLAMCH( 'S' ) / EPS
! 393: IERR = 0
! 394: *
! 395: WANTP = IJOB.EQ.1 .OR. IJOB.GE.4
! 396: WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4
! 397: WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5
! 398: WANTD = WANTD1 .OR. WANTD2
! 399: *
! 400: * Set M to the dimension of the specified pair of deflating
! 401: * subspaces.
! 402: *
! 403: M = 0
! 404: PAIR = .FALSE.
! 405: DO 10 K = 1, N
! 406: IF( PAIR ) THEN
! 407: PAIR = .FALSE.
! 408: ELSE
! 409: IF( K.LT.N ) THEN
! 410: IF( A( K+1, K ).EQ.ZERO ) THEN
! 411: IF( SELECT( K ) )
! 412: $ M = M + 1
! 413: ELSE
! 414: PAIR = .TRUE.
! 415: IF( SELECT( K ) .OR. SELECT( K+1 ) )
! 416: $ M = M + 2
! 417: END IF
! 418: ELSE
! 419: IF( SELECT( N ) )
! 420: $ M = M + 1
! 421: END IF
! 422: END IF
! 423: 10 CONTINUE
! 424: *
! 425: IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
! 426: LWMIN = MAX( 1, 4*N+16, 2*M*( N-M ) )
! 427: LIWMIN = MAX( 1, N+6 )
! 428: ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
! 429: LWMIN = MAX( 1, 4*N+16, 4*M*( N-M ) )
! 430: LIWMIN = MAX( 1, 2*M*( N-M ), N+6 )
! 431: ELSE
! 432: LWMIN = MAX( 1, 4*N+16 )
! 433: LIWMIN = 1
! 434: END IF
! 435: *
! 436: WORK( 1 ) = LWMIN
! 437: IWORK( 1 ) = LIWMIN
! 438: *
! 439: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
! 440: INFO = -22
! 441: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
! 442: INFO = -24
! 443: END IF
! 444: *
! 445: IF( INFO.NE.0 ) THEN
! 446: CALL XERBLA( 'DTGSEN', -INFO )
! 447: RETURN
! 448: ELSE IF( LQUERY ) THEN
! 449: RETURN
! 450: END IF
! 451: *
! 452: * Quick return if possible.
! 453: *
! 454: IF( M.EQ.N .OR. M.EQ.0 ) THEN
! 455: IF( WANTP ) THEN
! 456: PL = ONE
! 457: PR = ONE
! 458: END IF
! 459: IF( WANTD ) THEN
! 460: DSCALE = ZERO
! 461: DSUM = ONE
! 462: DO 20 I = 1, N
! 463: CALL DLASSQ( N, A( 1, I ), 1, DSCALE, DSUM )
! 464: CALL DLASSQ( N, B( 1, I ), 1, DSCALE, DSUM )
! 465: 20 CONTINUE
! 466: DIF( 1 ) = DSCALE*SQRT( DSUM )
! 467: DIF( 2 ) = DIF( 1 )
! 468: END IF
! 469: GO TO 60
! 470: END IF
! 471: *
! 472: * Collect the selected blocks at the top-left corner of (A, B).
! 473: *
! 474: KS = 0
! 475: PAIR = .FALSE.
! 476: DO 30 K = 1, N
! 477: IF( PAIR ) THEN
! 478: PAIR = .FALSE.
! 479: ELSE
! 480: *
! 481: SWAP = SELECT( K )
! 482: IF( K.LT.N ) THEN
! 483: IF( A( K+1, K ).NE.ZERO ) THEN
! 484: PAIR = .TRUE.
! 485: SWAP = SWAP .OR. SELECT( K+1 )
! 486: END IF
! 487: END IF
! 488: *
! 489: IF( SWAP ) THEN
! 490: KS = KS + 1
! 491: *
! 492: * Swap the K-th block to position KS.
! 493: * Perform the reordering of diagonal blocks in (A, B)
! 494: * by orthogonal transformation matrices and update
! 495: * Q and Z accordingly (if requested):
! 496: *
! 497: KK = K
! 498: IF( K.NE.KS )
! 499: $ CALL DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
! 500: $ Z, LDZ, KK, KS, WORK, LWORK, IERR )
! 501: *
! 502: IF( IERR.GT.0 ) THEN
! 503: *
! 504: * Swap is rejected: exit.
! 505: *
! 506: INFO = 1
! 507: IF( WANTP ) THEN
! 508: PL = ZERO
! 509: PR = ZERO
! 510: END IF
! 511: IF( WANTD ) THEN
! 512: DIF( 1 ) = ZERO
! 513: DIF( 2 ) = ZERO
! 514: END IF
! 515: GO TO 60
! 516: END IF
! 517: *
! 518: IF( PAIR )
! 519: $ KS = KS + 1
! 520: END IF
! 521: END IF
! 522: 30 CONTINUE
! 523: IF( WANTP ) THEN
! 524: *
! 525: * Solve generalized Sylvester equation for R and L
! 526: * and compute PL and PR.
! 527: *
! 528: N1 = M
! 529: N2 = N - M
! 530: I = N1 + 1
! 531: IJB = 0
! 532: CALL DLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
! 533: CALL DLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
! 534: $ N1 )
! 535: CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
! 536: $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1,
! 537: $ DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
! 538: $ LWORK-2*N1*N2, IWORK, IERR )
! 539: *
! 540: * Estimate the reciprocal of norms of "projections" onto left
! 541: * and right eigenspaces.
! 542: *
! 543: RDSCAL = ZERO
! 544: DSUM = ONE
! 545: CALL DLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM )
! 546: PL = RDSCAL*SQRT( DSUM )
! 547: IF( PL.EQ.ZERO ) THEN
! 548: PL = ONE
! 549: ELSE
! 550: PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) )
! 551: END IF
! 552: RDSCAL = ZERO
! 553: DSUM = ONE
! 554: CALL DLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM )
! 555: PR = RDSCAL*SQRT( DSUM )
! 556: IF( PR.EQ.ZERO ) THEN
! 557: PR = ONE
! 558: ELSE
! 559: PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) )
! 560: END IF
! 561: END IF
! 562: *
! 563: IF( WANTD ) THEN
! 564: *
! 565: * Compute estimates of Difu and Difl.
! 566: *
! 567: IF( WANTD1 ) THEN
! 568: N1 = M
! 569: N2 = N - M
! 570: I = N1 + 1
! 571: IJB = IDIFJB
! 572: *
! 573: * Frobenius norm-based Difu-estimate.
! 574: *
! 575: CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
! 576: $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
! 577: $ N1, DSCALE, DIF( 1 ), WORK( 2*N1*N2+1 ),
! 578: $ LWORK-2*N1*N2, IWORK, IERR )
! 579: *
! 580: * Frobenius norm-based Difl-estimate.
! 581: *
! 582: CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
! 583: $ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
! 584: $ N2, DSCALE, DIF( 2 ), WORK( 2*N1*N2+1 ),
! 585: $ LWORK-2*N1*N2, IWORK, IERR )
! 586: ELSE
! 587: *
! 588: *
! 589: * Compute 1-norm-based estimates of Difu and Difl using
! 590: * reversed communication with DLACN2. In each step a
! 591: * generalized Sylvester equation or a transposed variant
! 592: * is solved.
! 593: *
! 594: KASE = 0
! 595: N1 = M
! 596: N2 = N - M
! 597: I = N1 + 1
! 598: IJB = 0
! 599: MN2 = 2*N1*N2
! 600: *
! 601: * 1-norm-based estimate of Difu.
! 602: *
! 603: 40 CONTINUE
! 604: CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 1 ),
! 605: $ KASE, ISAVE )
! 606: IF( KASE.NE.0 ) THEN
! 607: IF( KASE.EQ.1 ) THEN
! 608: *
! 609: * Solve generalized Sylvester equation.
! 610: *
! 611: CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
! 612: $ WORK, N1, B, LDB, B( I, I ), LDB,
! 613: $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
! 614: $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
! 615: $ IERR )
! 616: ELSE
! 617: *
! 618: * Solve the transposed variant.
! 619: *
! 620: CALL DTGSYL( 'T', IJB, N1, N2, A, LDA, A( I, I ), LDA,
! 621: $ WORK, N1, B, LDB, B( I, I ), LDB,
! 622: $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
! 623: $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
! 624: $ IERR )
! 625: END IF
! 626: GO TO 40
! 627: END IF
! 628: DIF( 1 ) = DSCALE / DIF( 1 )
! 629: *
! 630: * 1-norm-based estimate of Difl.
! 631: *
! 632: 50 CONTINUE
! 633: CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 2 ),
! 634: $ KASE, ISAVE )
! 635: IF( KASE.NE.0 ) THEN
! 636: IF( KASE.EQ.1 ) THEN
! 637: *
! 638: * Solve generalized Sylvester equation.
! 639: *
! 640: CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
! 641: $ WORK, N2, B( I, I ), LDB, B, LDB,
! 642: $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
! 643: $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
! 644: $ IERR )
! 645: ELSE
! 646: *
! 647: * Solve the transposed variant.
! 648: *
! 649: CALL DTGSYL( 'T', IJB, N2, N1, A( I, I ), LDA, A, LDA,
! 650: $ WORK, N2, B( I, I ), LDB, B, LDB,
! 651: $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
! 652: $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
! 653: $ IERR )
! 654: END IF
! 655: GO TO 50
! 656: END IF
! 657: DIF( 2 ) = DSCALE / DIF( 2 )
! 658: *
! 659: END IF
! 660: END IF
! 661: *
! 662: 60 CONTINUE
! 663: *
! 664: * Compute generalized eigenvalues of reordered pair (A, B) and
! 665: * normalize the generalized Schur form.
! 666: *
! 667: PAIR = .FALSE.
! 668: DO 80 K = 1, N
! 669: IF( PAIR ) THEN
! 670: PAIR = .FALSE.
! 671: ELSE
! 672: *
! 673: IF( K.LT.N ) THEN
! 674: IF( A( K+1, K ).NE.ZERO ) THEN
! 675: PAIR = .TRUE.
! 676: END IF
! 677: END IF
! 678: *
! 679: IF( PAIR ) THEN
! 680: *
! 681: * Compute the eigenvalue(s) at position K.
! 682: *
! 683: WORK( 1 ) = A( K, K )
! 684: WORK( 2 ) = A( K+1, K )
! 685: WORK( 3 ) = A( K, K+1 )
! 686: WORK( 4 ) = A( K+1, K+1 )
! 687: WORK( 5 ) = B( K, K )
! 688: WORK( 6 ) = B( K+1, K )
! 689: WORK( 7 ) = B( K, K+1 )
! 690: WORK( 8 ) = B( K+1, K+1 )
! 691: CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA( K ),
! 692: $ BETA( K+1 ), ALPHAR( K ), ALPHAR( K+1 ),
! 693: $ ALPHAI( K ) )
! 694: ALPHAI( K+1 ) = -ALPHAI( K )
! 695: *
! 696: ELSE
! 697: *
! 698: IF( SIGN( ONE, B( K, K ) ).LT.ZERO ) THEN
! 699: *
! 700: * If B(K,K) is negative, make it positive
! 701: *
! 702: DO 70 I = 1, N
! 703: A( K, I ) = -A( K, I )
! 704: B( K, I ) = -B( K, I )
! 705: IF( WANTQ ) Q( I, K ) = -Q( I, K )
! 706: 70 CONTINUE
! 707: END IF
! 708: *
! 709: ALPHAR( K ) = A( K, K )
! 710: ALPHAI( K ) = ZERO
! 711: BETA( K ) = B( K, K )
! 712: *
! 713: END IF
! 714: END IF
! 715: 80 CONTINUE
! 716: *
! 717: WORK( 1 ) = LWMIN
! 718: IWORK( 1 ) = LIWMIN
! 719: *
! 720: RETURN
! 721: *
! 722: * End of DTGSEN
! 723: *
! 724: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtgsen.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack