Annotation of rpl/lapack/lapack/dtrexc.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
! 2: $ 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: CHARACTER COMPQ
! 11: INTEGER IFST, ILST, INFO, LDQ, LDT, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
! 15: * ..
! 16: *
! 17: * Purpose
! 18: * =======
! 19: *
! 20: * DTREXC reorders the real Schur factorization of a real matrix
! 21: * A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
! 22: * moved to row ILST.
! 23: *
! 24: * The real Schur form T is reordered by an orthogonal similarity
! 25: * transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
! 26: * is updated by postmultiplying it with Z.
! 27: *
! 28: * T must be in Schur canonical form (as returned by DHSEQR), that is,
! 29: * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
! 30: * 2-by-2 diagonal block has its diagonal elements equal and its
! 31: * off-diagonal elements of opposite sign.
! 32: *
! 33: * Arguments
! 34: * =========
! 35: *
! 36: * COMPQ (input) CHARACTER*1
! 37: * = 'V': update the matrix Q of Schur vectors;
! 38: * = 'N': do not update Q.
! 39: *
! 40: * N (input) INTEGER
! 41: * The order of the matrix T. N >= 0.
! 42: *
! 43: * T (input/output) DOUBLE PRECISION array, dimension (LDT,N)
! 44: * On entry, the upper quasi-triangular matrix T, in Schur
! 45: * Schur canonical form.
! 46: * On exit, the reordered upper quasi-triangular matrix, again
! 47: * in Schur canonical form.
! 48: *
! 49: * LDT (input) INTEGER
! 50: * The leading dimension of the array T. LDT >= max(1,N).
! 51: *
! 52: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
! 53: * On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
! 54: * On exit, if COMPQ = 'V', Q has been postmultiplied by the
! 55: * orthogonal transformation matrix Z which reorders T.
! 56: * If COMPQ = 'N', Q is not referenced.
! 57: *
! 58: * LDQ (input) INTEGER
! 59: * The leading dimension of the array Q. LDQ >= max(1,N).
! 60: *
! 61: * IFST (input/output) INTEGER
! 62: * ILST (input/output) INTEGER
! 63: * Specify the reordering of the diagonal blocks of T.
! 64: * The block with row index IFST is moved to row ILST, by a
! 65: * sequence of transpositions between adjacent blocks.
! 66: * On exit, if IFST pointed on entry to the second row of a
! 67: * 2-by-2 block, it is changed to point to the first row; ILST
! 68: * always points to the first row of the block in its final
! 69: * position (which may differ from its input value by +1 or -1).
! 70: * 1 <= IFST <= N; 1 <= ILST <= N.
! 71: *
! 72: * WORK (workspace) DOUBLE PRECISION array, dimension (N)
! 73: *
! 74: * INFO (output) INTEGER
! 75: * = 0: successful exit
! 76: * < 0: if INFO = -i, the i-th argument had an illegal value
! 77: * = 1: two adjacent blocks were too close to swap (the problem
! 78: * is very ill-conditioned); T may have been partially
! 79: * reordered, and ILST points to the first row of the
! 80: * current position of the block being moved.
! 81: *
! 82: * =====================================================================
! 83: *
! 84: * .. Parameters ..
! 85: DOUBLE PRECISION ZERO
! 86: PARAMETER ( ZERO = 0.0D+0 )
! 87: * ..
! 88: * .. Local Scalars ..
! 89: LOGICAL WANTQ
! 90: INTEGER HERE, NBF, NBL, NBNEXT
! 91: * ..
! 92: * .. External Functions ..
! 93: LOGICAL LSAME
! 94: EXTERNAL LSAME
! 95: * ..
! 96: * .. External Subroutines ..
! 97: EXTERNAL DLAEXC, XERBLA
! 98: * ..
! 99: * .. Intrinsic Functions ..
! 100: INTRINSIC MAX
! 101: * ..
! 102: * .. Executable Statements ..
! 103: *
! 104: * Decode and test the input arguments.
! 105: *
! 106: INFO = 0
! 107: WANTQ = LSAME( COMPQ, 'V' )
! 108: IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
! 109: INFO = -1
! 110: ELSE IF( N.LT.0 ) THEN
! 111: INFO = -2
! 112: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
! 113: INFO = -4
! 114: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
! 115: INFO = -6
! 116: ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
! 117: INFO = -7
! 118: ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
! 119: INFO = -8
! 120: END IF
! 121: IF( INFO.NE.0 ) THEN
! 122: CALL XERBLA( 'DTREXC', -INFO )
! 123: RETURN
! 124: END IF
! 125: *
! 126: * Quick return if possible
! 127: *
! 128: IF( N.LE.1 )
! 129: $ RETURN
! 130: *
! 131: * Determine the first row of specified block
! 132: * and find out it is 1 by 1 or 2 by 2.
! 133: *
! 134: IF( IFST.GT.1 ) THEN
! 135: IF( T( IFST, IFST-1 ).NE.ZERO )
! 136: $ IFST = IFST - 1
! 137: END IF
! 138: NBF = 1
! 139: IF( IFST.LT.N ) THEN
! 140: IF( T( IFST+1, IFST ).NE.ZERO )
! 141: $ NBF = 2
! 142: END IF
! 143: *
! 144: * Determine the first row of the final block
! 145: * and find out it is 1 by 1 or 2 by 2.
! 146: *
! 147: IF( ILST.GT.1 ) THEN
! 148: IF( T( ILST, ILST-1 ).NE.ZERO )
! 149: $ ILST = ILST - 1
! 150: END IF
! 151: NBL = 1
! 152: IF( ILST.LT.N ) THEN
! 153: IF( T( ILST+1, ILST ).NE.ZERO )
! 154: $ NBL = 2
! 155: END IF
! 156: *
! 157: IF( IFST.EQ.ILST )
! 158: $ RETURN
! 159: *
! 160: IF( IFST.LT.ILST ) THEN
! 161: *
! 162: * Update ILST
! 163: *
! 164: IF( NBF.EQ.2 .AND. NBL.EQ.1 )
! 165: $ ILST = ILST - 1
! 166: IF( NBF.EQ.1 .AND. NBL.EQ.2 )
! 167: $ ILST = ILST + 1
! 168: *
! 169: HERE = IFST
! 170: *
! 171: 10 CONTINUE
! 172: *
! 173: * Swap block with next one below
! 174: *
! 175: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
! 176: *
! 177: * Current block either 1 by 1 or 2 by 2
! 178: *
! 179: NBNEXT = 1
! 180: IF( HERE+NBF+1.LE.N ) THEN
! 181: IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
! 182: $ NBNEXT = 2
! 183: END IF
! 184: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
! 185: $ WORK, INFO )
! 186: IF( INFO.NE.0 ) THEN
! 187: ILST = HERE
! 188: RETURN
! 189: END IF
! 190: HERE = HERE + NBNEXT
! 191: *
! 192: * Test if 2 by 2 block breaks into two 1 by 1 blocks
! 193: *
! 194: IF( NBF.EQ.2 ) THEN
! 195: IF( T( HERE+1, HERE ).EQ.ZERO )
! 196: $ NBF = 3
! 197: END IF
! 198: *
! 199: ELSE
! 200: *
! 201: * Current block consists of two 1 by 1 blocks each of which
! 202: * must be swapped individually
! 203: *
! 204: NBNEXT = 1
! 205: IF( HERE+3.LE.N ) THEN
! 206: IF( T( HERE+3, HERE+2 ).NE.ZERO )
! 207: $ NBNEXT = 2
! 208: END IF
! 209: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
! 210: $ WORK, INFO )
! 211: IF( INFO.NE.0 ) THEN
! 212: ILST = HERE
! 213: RETURN
! 214: END IF
! 215: IF( NBNEXT.EQ.1 ) THEN
! 216: *
! 217: * Swap two 1 by 1 blocks, no problems possible
! 218: *
! 219: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
! 220: $ WORK, INFO )
! 221: HERE = HERE + 1
! 222: ELSE
! 223: *
! 224: * Recompute NBNEXT in case 2 by 2 split
! 225: *
! 226: IF( T( HERE+2, HERE+1 ).EQ.ZERO )
! 227: $ NBNEXT = 1
! 228: IF( NBNEXT.EQ.2 ) THEN
! 229: *
! 230: * 2 by 2 Block did not split
! 231: *
! 232: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
! 233: $ NBNEXT, WORK, INFO )
! 234: IF( INFO.NE.0 ) THEN
! 235: ILST = HERE
! 236: RETURN
! 237: END IF
! 238: HERE = HERE + 2
! 239: ELSE
! 240: *
! 241: * 2 by 2 Block did split
! 242: *
! 243: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
! 244: $ WORK, INFO )
! 245: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
! 246: $ WORK, INFO )
! 247: HERE = HERE + 2
! 248: END IF
! 249: END IF
! 250: END IF
! 251: IF( HERE.LT.ILST )
! 252: $ GO TO 10
! 253: *
! 254: ELSE
! 255: *
! 256: HERE = IFST
! 257: 20 CONTINUE
! 258: *
! 259: * Swap block with next one above
! 260: *
! 261: IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
! 262: *
! 263: * Current block either 1 by 1 or 2 by 2
! 264: *
! 265: NBNEXT = 1
! 266: IF( HERE.GE.3 ) THEN
! 267: IF( T( HERE-1, HERE-2 ).NE.ZERO )
! 268: $ NBNEXT = 2
! 269: END IF
! 270: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
! 271: $ NBF, WORK, INFO )
! 272: IF( INFO.NE.0 ) THEN
! 273: ILST = HERE
! 274: RETURN
! 275: END IF
! 276: HERE = HERE - NBNEXT
! 277: *
! 278: * Test if 2 by 2 block breaks into two 1 by 1 blocks
! 279: *
! 280: IF( NBF.EQ.2 ) THEN
! 281: IF( T( HERE+1, HERE ).EQ.ZERO )
! 282: $ NBF = 3
! 283: END IF
! 284: *
! 285: ELSE
! 286: *
! 287: * Current block consists of two 1 by 1 blocks each of which
! 288: * must be swapped individually
! 289: *
! 290: NBNEXT = 1
! 291: IF( HERE.GE.3 ) THEN
! 292: IF( T( HERE-1, HERE-2 ).NE.ZERO )
! 293: $ NBNEXT = 2
! 294: END IF
! 295: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
! 296: $ 1, WORK, INFO )
! 297: IF( INFO.NE.0 ) THEN
! 298: ILST = HERE
! 299: RETURN
! 300: END IF
! 301: IF( NBNEXT.EQ.1 ) THEN
! 302: *
! 303: * Swap two 1 by 1 blocks, no problems possible
! 304: *
! 305: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
! 306: $ WORK, INFO )
! 307: HERE = HERE - 1
! 308: ELSE
! 309: *
! 310: * Recompute NBNEXT in case 2 by 2 split
! 311: *
! 312: IF( T( HERE, HERE-1 ).EQ.ZERO )
! 313: $ NBNEXT = 1
! 314: IF( NBNEXT.EQ.2 ) THEN
! 315: *
! 316: * 2 by 2 Block did not split
! 317: *
! 318: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
! 319: $ WORK, INFO )
! 320: IF( INFO.NE.0 ) THEN
! 321: ILST = HERE
! 322: RETURN
! 323: END IF
! 324: HERE = HERE - 2
! 325: ELSE
! 326: *
! 327: * 2 by 2 Block did split
! 328: *
! 329: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
! 330: $ WORK, INFO )
! 331: CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
! 332: $ WORK, INFO )
! 333: HERE = HERE - 2
! 334: END IF
! 335: END IF
! 336: END IF
! 337: IF( HERE.GT.ILST )
! 338: $ GO TO 20
! 339: END IF
! 340: ILST = HERE
! 341: *
! 342: RETURN
! 343: *
! 344: * End of DTREXC
! 345: *
! 346: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrexc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack