SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
$ INFO )
*
* -- LAPACK routine (version 3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
* .. Scalar Arguments ..
CHARACTER COMPQ
INTEGER IFST, ILST, INFO, LDQ, LDT, N
* ..
* .. Array Arguments ..
DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
* ..
*
* Purpose
* =======
*
* DTREXC reorders the real Schur factorization of a real matrix
* A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
* moved to row ILST.
*
* The real Schur form T is reordered by an orthogonal similarity
* transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
* is updated by postmultiplying it with Z.
*
* T must be in Schur canonical form (as returned by DHSEQR), that is,
* block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
* 2-by-2 diagonal block has its diagonal elements equal and its
* off-diagonal elements of opposite sign.
*
* Arguments
* =========
*
* COMPQ (input) CHARACTER*1
* = 'V': update the matrix Q of Schur vectors;
* = 'N': do not update Q.
*
* N (input) INTEGER
* The order of the matrix T. N >= 0.
*
* T (input/output) DOUBLE PRECISION array, dimension (LDT,N)
* On entry, the upper quasi-triangular matrix T, in Schur
* Schur canonical form.
* On exit, the reordered upper quasi-triangular matrix, again
* in Schur canonical form.
*
* LDT (input) INTEGER
* The leading dimension of the array T. LDT >= max(1,N).
*
* Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
* On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
* On exit, if COMPQ = 'V', Q has been postmultiplied by the
* orthogonal transformation matrix Z which reorders T.
* If COMPQ = 'N', Q is not referenced.
*
* LDQ (input) INTEGER
* The leading dimension of the array Q. LDQ >= max(1,N).
*
* IFST (input/output) INTEGER
* ILST (input/output) INTEGER
* Specify the reordering of the diagonal blocks of T.
* The block with row index IFST is moved to row ILST, by a
* sequence of transpositions between adjacent blocks.
* On exit, if IFST pointed on entry to the second row of a
* 2-by-2 block, it is changed to point to the first row; ILST
* always points to the first row of the block in its final
* position (which may differ from its input value by +1 or -1).
* 1 <= IFST <= N; 1 <= ILST <= N.
*
* WORK (workspace) DOUBLE PRECISION array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* = 1: two adjacent blocks were too close to swap (the problem
* is very ill-conditioned); T may have been partially
* reordered, and ILST points to the first row of the
* current position of the block being moved.
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL WANTQ
INTEGER HERE, NBF, NBL, NBNEXT
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL DLAEXC, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Decode and test the input arguments.
*
INFO = 0
WANTQ = LSAME( COMPQ, 'V' )
IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
INFO = -4
ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
INFO = -6
ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
INFO = -7
ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
INFO = -8
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DTREXC', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.LE.1 )
$ RETURN
*
* Determine the first row of specified block
* and find out it is 1 by 1 or 2 by 2.
*
IF( IFST.GT.1 ) THEN
IF( T( IFST, IFST-1 ).NE.ZERO )
$ IFST = IFST - 1
END IF
NBF = 1
IF( IFST.LT.N ) THEN
IF( T( IFST+1, IFST ).NE.ZERO )
$ NBF = 2
END IF
*
* Determine the first row of the final block
* and find out it is 1 by 1 or 2 by 2.
*
IF( ILST.GT.1 ) THEN
IF( T( ILST, ILST-1 ).NE.ZERO )
$ ILST = ILST - 1
END IF
NBL = 1
IF( ILST.LT.N ) THEN
IF( T( ILST+1, ILST ).NE.ZERO )
$ NBL = 2
END IF
*
IF( IFST.EQ.ILST )
$ RETURN
*
IF( IFST.LT.ILST ) THEN
*
* Update ILST
*
IF( NBF.EQ.2 .AND. NBL.EQ.1 )
$ ILST = ILST - 1
IF( NBF.EQ.1 .AND. NBL.EQ.2 )
$ ILST = ILST + 1
*
HERE = IFST
*
10 CONTINUE
*
* Swap block with next one below
*
IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
*
* Current block either 1 by 1 or 2 by 2
*
NBNEXT = 1
IF( HERE+NBF+1.LE.N ) THEN
IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
$ NBNEXT = 2
END IF
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
$ WORK, INFO )
IF( INFO.NE.0 ) THEN
ILST = HERE
RETURN
END IF
HERE = HERE + NBNEXT
*
* Test if 2 by 2 block breaks into two 1 by 1 blocks
*
IF( NBF.EQ.2 ) THEN
IF( T( HERE+1, HERE ).EQ.ZERO )
$ NBF = 3
END IF
*
ELSE
*
* Current block consists of two 1 by 1 blocks each of which
* must be swapped individually
*
NBNEXT = 1
IF( HERE+3.LE.N ) THEN
IF( T( HERE+3, HERE+2 ).NE.ZERO )
$ NBNEXT = 2
END IF
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
$ WORK, INFO )
IF( INFO.NE.0 ) THEN
ILST = HERE
RETURN
END IF
IF( NBNEXT.EQ.1 ) THEN
*
* Swap two 1 by 1 blocks, no problems possible
*
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
$ WORK, INFO )
HERE = HERE + 1
ELSE
*
* Recompute NBNEXT in case 2 by 2 split
*
IF( T( HERE+2, HERE+1 ).EQ.ZERO )
$ NBNEXT = 1
IF( NBNEXT.EQ.2 ) THEN
*
* 2 by 2 Block did not split
*
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
$ NBNEXT, WORK, INFO )
IF( INFO.NE.0 ) THEN
ILST = HERE
RETURN
END IF
HERE = HERE + 2
ELSE
*
* 2 by 2 Block did split
*
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
$ WORK, INFO )
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
$ WORK, INFO )
HERE = HERE + 2
END IF
END IF
END IF
IF( HERE.LT.ILST )
$ GO TO 10
*
ELSE
*
HERE = IFST
20 CONTINUE
*
* Swap block with next one above
*
IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
*
* Current block either 1 by 1 or 2 by 2
*
NBNEXT = 1
IF( HERE.GE.3 ) THEN
IF( T( HERE-1, HERE-2 ).NE.ZERO )
$ NBNEXT = 2
END IF
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
$ NBF, WORK, INFO )
IF( INFO.NE.0 ) THEN
ILST = HERE
RETURN
END IF
HERE = HERE - NBNEXT
*
* Test if 2 by 2 block breaks into two 1 by 1 blocks
*
IF( NBF.EQ.2 ) THEN
IF( T( HERE+1, HERE ).EQ.ZERO )
$ NBF = 3
END IF
*
ELSE
*
* Current block consists of two 1 by 1 blocks each of which
* must be swapped individually
*
NBNEXT = 1
IF( HERE.GE.3 ) THEN
IF( T( HERE-1, HERE-2 ).NE.ZERO )
$ NBNEXT = 2
END IF
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
$ 1, WORK, INFO )
IF( INFO.NE.0 ) THEN
ILST = HERE
RETURN
END IF
IF( NBNEXT.EQ.1 ) THEN
*
* Swap two 1 by 1 blocks, no problems possible
*
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
$ WORK, INFO )
HERE = HERE - 1
ELSE
*
* Recompute NBNEXT in case 2 by 2 split
*
IF( T( HERE, HERE-1 ).EQ.ZERO )
$ NBNEXT = 1
IF( NBNEXT.EQ.2 ) THEN
*
* 2 by 2 Block did not split
*
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
$ WORK, INFO )
IF( INFO.NE.0 ) THEN
ILST = HERE
RETURN
END IF
HERE = HERE - 2
ELSE
*
* 2 by 2 Block did split
*
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
$ WORK, INFO )
CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
$ WORK, INFO )
HERE = HERE - 2
END IF
END IF
END IF
IF( HERE.GT.ILST )
$ GO TO 20
END IF
ILST = HERE
*
RETURN
*
* End of DTREXC
*
END