Annotation of rpl/lapack/lapack/dtgexc.f, revision 1.1
1.1 ! bertrand 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) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 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 (LDZ,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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtgexc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack