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