![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 2: $ LDZ, IFST, ILST, WORK, LWORK, INFO ) 3: * 4: * -- LAPACK routine (version 3.2.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * June 2010 8: * 9: * .. Scalar Arguments .. 10: LOGICAL WANTQ, WANTZ 11: INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N 12: * .. 13: * .. Array Arguments .. 14: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 15: $ WORK( * ), Z( LDZ, * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * DTGEXC reorders the generalized real Schur decomposition of a real 22: * matrix pair (A,B) using an orthogonal equivalence transformation 23: * 24: * (A, B) = Q * (A, B) * Z', 25: * 26: * so that the diagonal block of (A, B) with row index IFST is moved 27: * to row ILST. 28: * 29: * (A, B) must be in generalized real Schur canonical form (as returned 30: * by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2 31: * diagonal blocks. B is upper triangular. 32: * 33: * Optionally, the matrices Q and Z of generalized Schur vectors are 34: * updated. 35: * 36: * Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)' 37: * Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)' 38: * 39: * 40: * Arguments 41: * ========= 42: * 43: * WANTQ (input) LOGICAL 44: * .TRUE. : update the left transformation matrix Q; 45: * .FALSE.: do not update Q. 46: * 47: * WANTZ (input) LOGICAL 48: * .TRUE. : update the right transformation matrix Z; 49: * .FALSE.: do not update Z. 50: * 51: * N (input) INTEGER 52: * The order of the matrices A and B. N >= 0. 53: * 54: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 55: * On entry, the matrix A in generalized real Schur canonical 56: * form. 57: * On exit, the updated matrix A, again in generalized 58: * real Schur canonical form. 59: * 60: * LDA (input) INTEGER 61: * The leading dimension of the array A. LDA >= max(1,N). 62: * 63: * B (input/output) DOUBLE PRECISION array, dimension (LDB,N) 64: * On entry, the matrix B in generalized real Schur canonical 65: * form (A,B). 66: * On exit, the updated matrix B, again in generalized 67: * real Schur canonical form (A,B). 68: * 69: * LDB (input) INTEGER 70: * The leading dimension of the array B. LDB >= max(1,N). 71: * 72: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) 73: * On entry, if WANTQ = .TRUE., the orthogonal matrix Q. 74: * On exit, the updated matrix Q. 75: * If WANTQ = .FALSE., Q is not referenced. 76: * 77: * LDQ (input) INTEGER 78: * The leading dimension of the array Q. LDQ >= 1. 79: * If WANTQ = .TRUE., LDQ >= N. 80: * 81: * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) 82: * On entry, if WANTZ = .TRUE., the orthogonal matrix Z. 83: * On exit, the updated matrix Z. 84: * If WANTZ = .FALSE., Z is not referenced. 85: * 86: * LDZ (input) INTEGER 87: * The leading dimension of the array Z. LDZ >= 1. 88: * If WANTZ = .TRUE., LDZ >= N. 89: * 90: * IFST (input/output) INTEGER 91: * ILST (input/output) INTEGER 92: * Specify the reordering of the diagonal blocks of (A, B). 93: * The block with row index IFST is moved to row ILST, by a 94: * sequence of swapping between adjacent blocks. 95: * On exit, if IFST pointed on entry to the second row of 96: * a 2-by-2 block, it is changed to point to the first row; 97: * ILST always points to the first row of the block in its 98: * final position (which may differ from its input value by 99: * +1 or -1). 1 <= IFST, ILST <= N. 100: * 101: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 102: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 103: * 104: * LWORK (input) INTEGER 105: * The dimension of the array WORK. 106: * LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16. 107: * 108: * If LWORK = -1, then a workspace query is assumed; the routine 109: * only calculates the optimal size of the WORK array, returns 110: * this value as the first entry of the WORK array, and no error 111: * message related to LWORK is issued by XERBLA. 112: * 113: * INFO (output) INTEGER 114: * =0: successful exit. 115: * <0: if INFO = -i, the i-th argument had an illegal value. 116: * =1: The transformed matrix pair (A, B) would be too far 117: * from generalized Schur form; the problem is ill- 118: * conditioned. (A, B) may have been partially reordered, 119: * and ILST points to the first row of the current 120: * position of the block being moved. 121: * 122: * Further Details 123: * =============== 124: * 125: * Based on contributions by 126: * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 127: * Umea University, S-901 87 Umea, Sweden. 128: * 129: * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 130: * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 131: * M.S. Moonen et al (eds), Linear Algebra for Large Scale and 132: * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 133: * 134: * ===================================================================== 135: * 136: * .. Parameters .. 137: DOUBLE PRECISION ZERO 138: PARAMETER ( ZERO = 0.0D+0 ) 139: * .. 140: * .. Local Scalars .. 141: LOGICAL LQUERY 142: INTEGER HERE, LWMIN, NBF, NBL, NBNEXT 143: * .. 144: * .. External Subroutines .. 145: EXTERNAL DTGEX2, XERBLA 146: * .. 147: * .. Intrinsic Functions .. 148: INTRINSIC MAX 149: * .. 150: * .. Executable Statements .. 151: * 152: * Decode and test input arguments. 153: * 154: INFO = 0 155: LQUERY = ( LWORK.EQ.-1 ) 156: IF( N.LT.0 ) THEN 157: INFO = -3 158: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 159: INFO = -5 160: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 161: INFO = -7 162: ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN 163: INFO = -9 164: ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN 165: INFO = -11 166: ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN 167: INFO = -12 168: ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN 169: INFO = -13 170: END IF 171: * 172: IF( INFO.EQ.0 ) THEN 173: IF( N.LE.1 ) THEN 174: LWMIN = 1 175: ELSE 176: LWMIN = 4*N + 16 177: END IF 178: WORK(1) = LWMIN 179: * 180: IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN 181: INFO = -15 182: END IF 183: END IF 184: * 185: IF( INFO.NE.0 ) THEN 186: CALL XERBLA( 'DTGEXC', -INFO ) 187: RETURN 188: ELSE IF( LQUERY ) THEN 189: RETURN 190: END IF 191: * 192: * Quick return if possible 193: * 194: IF( N.LE.1 ) 195: $ RETURN 196: * 197: * Determine the first row of the specified block and find out 198: * if it is 1-by-1 or 2-by-2. 199: * 200: IF( IFST.GT.1 ) THEN 201: IF( A( IFST, IFST-1 ).NE.ZERO ) 202: $ IFST = IFST - 1 203: END IF 204: NBF = 1 205: IF( IFST.LT.N ) THEN 206: IF( A( IFST+1, IFST ).NE.ZERO ) 207: $ NBF = 2 208: END IF 209: * 210: * Determine the first row of the final block 211: * and find out if it is 1-by-1 or 2-by-2. 212: * 213: IF( ILST.GT.1 ) THEN 214: IF( A( ILST, ILST-1 ).NE.ZERO ) 215: $ ILST = ILST - 1 216: END IF 217: NBL = 1 218: IF( ILST.LT.N ) THEN 219: IF( A( ILST+1, ILST ).NE.ZERO ) 220: $ NBL = 2 221: END IF 222: IF( IFST.EQ.ILST ) 223: $ RETURN 224: * 225: IF( IFST.LT.ILST ) THEN 226: * 227: * Update ILST. 228: * 229: IF( NBF.EQ.2 .AND. NBL.EQ.1 ) 230: $ ILST = ILST - 1 231: IF( NBF.EQ.1 .AND. NBL.EQ.2 ) 232: $ ILST = ILST + 1 233: * 234: HERE = IFST 235: * 236: 10 CONTINUE 237: * 238: * Swap with next one below. 239: * 240: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN 241: * 242: * Current block either 1-by-1 or 2-by-2. 243: * 244: NBNEXT = 1 245: IF( HERE+NBF+1.LE.N ) THEN 246: IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO ) 247: $ NBNEXT = 2 248: END IF 249: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 250: $ LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO ) 251: IF( INFO.NE.0 ) THEN 252: ILST = HERE 253: RETURN 254: END IF 255: HERE = HERE + NBNEXT 256: * 257: * Test if 2-by-2 block breaks into two 1-by-1 blocks. 258: * 259: IF( NBF.EQ.2 ) THEN 260: IF( A( HERE+1, HERE ).EQ.ZERO ) 261: $ NBF = 3 262: END IF 263: * 264: ELSE 265: * 266: * Current block consists of two 1-by-1 blocks, each of which 267: * must be swapped individually. 268: * 269: NBNEXT = 1 270: IF( HERE+3.LE.N ) THEN 271: IF( A( HERE+3, HERE+2 ).NE.ZERO ) 272: $ NBNEXT = 2 273: END IF 274: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 275: $ LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO ) 276: IF( INFO.NE.0 ) THEN 277: ILST = HERE 278: RETURN 279: END IF 280: IF( NBNEXT.EQ.1 ) THEN 281: * 282: * Swap two 1-by-1 blocks. 283: * 284: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 285: $ LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 286: IF( INFO.NE.0 ) THEN 287: ILST = HERE 288: RETURN 289: END IF 290: HERE = HERE + 1 291: * 292: ELSE 293: * 294: * Recompute NBNEXT in case of 2-by-2 split. 295: * 296: IF( A( HERE+2, HERE+1 ).EQ.ZERO ) 297: $ NBNEXT = 1 298: IF( NBNEXT.EQ.2 ) THEN 299: * 300: * 2-by-2 block did not split. 301: * 302: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 303: $ Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK, 304: $ INFO ) 305: IF( INFO.NE.0 ) THEN 306: ILST = HERE 307: RETURN 308: END IF 309: HERE = HERE + 2 310: ELSE 311: * 312: * 2-by-2 block did split. 313: * 314: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 315: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 316: IF( INFO.NE.0 ) THEN 317: ILST = HERE 318: RETURN 319: END IF 320: HERE = HERE + 1 321: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 322: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 323: IF( INFO.NE.0 ) THEN 324: ILST = HERE 325: RETURN 326: END IF 327: HERE = HERE + 1 328: END IF 329: * 330: END IF 331: END IF 332: IF( HERE.LT.ILST ) 333: $ GO TO 10 334: ELSE 335: HERE = IFST 336: * 337: 20 CONTINUE 338: * 339: * Swap with next one below. 340: * 341: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN 342: * 343: * Current block either 1-by-1 or 2-by-2. 344: * 345: NBNEXT = 1 346: IF( HERE.GE.3 ) THEN 347: IF( A( HERE-1, HERE-2 ).NE.ZERO ) 348: $ NBNEXT = 2 349: END IF 350: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 351: $ LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK, 352: $ INFO ) 353: IF( INFO.NE.0 ) THEN 354: ILST = HERE 355: RETURN 356: END IF 357: HERE = HERE - NBNEXT 358: * 359: * Test if 2-by-2 block breaks into two 1-by-1 blocks. 360: * 361: IF( NBF.EQ.2 ) THEN 362: IF( A( HERE+1, HERE ).EQ.ZERO ) 363: $ NBF = 3 364: END IF 365: * 366: ELSE 367: * 368: * Current block consists of two 1-by-1 blocks, each of which 369: * must be swapped individually. 370: * 371: NBNEXT = 1 372: IF( HERE.GE.3 ) THEN 373: IF( A( HERE-1, HERE-2 ).NE.ZERO ) 374: $ NBNEXT = 2 375: END IF 376: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 377: $ LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK, 378: $ INFO ) 379: IF( INFO.NE.0 ) THEN 380: ILST = HERE 381: RETURN 382: END IF 383: IF( NBNEXT.EQ.1 ) THEN 384: * 385: * Swap two 1-by-1 blocks. 386: * 387: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 388: $ LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO ) 389: IF( INFO.NE.0 ) THEN 390: ILST = HERE 391: RETURN 392: END IF 393: HERE = HERE - 1 394: ELSE 395: * 396: * Recompute NBNEXT in case of 2-by-2 split. 397: * 398: IF( A( HERE, HERE-1 ).EQ.ZERO ) 399: $ NBNEXT = 1 400: IF( NBNEXT.EQ.2 ) THEN 401: * 402: * 2-by-2 block did not split. 403: * 404: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 405: $ Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO ) 406: IF( INFO.NE.0 ) THEN 407: ILST = HERE 408: RETURN 409: END IF 410: HERE = HERE - 2 411: ELSE 412: * 413: * 2-by-2 block did split. 414: * 415: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 416: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 417: IF( INFO.NE.0 ) THEN 418: ILST = HERE 419: RETURN 420: END IF 421: HERE = HERE - 1 422: CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 423: $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 424: IF( INFO.NE.0 ) THEN 425: ILST = HERE 426: RETURN 427: END IF 428: HERE = HERE - 1 429: END IF 430: END IF 431: END IF 432: IF( HERE.GT.ILST ) 433: $ GO TO 20 434: END IF 435: ILST = HERE 436: WORK( 1 ) = LWMIN 437: RETURN 438: * 439: * End of DTGEXC 440: * 441: END