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