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