441 lines
14 KiB
FortranFixed
441 lines
14 KiB
FortranFixed
SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
|
|
$ LDZ, IFST, ILST, WORK, LWORK, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.1) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
|
* November 2006
|
|
*
|
|
* .. Scalar Arguments ..
|
|
LOGICAL WANTQ, WANTZ
|
|
INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
|
|
$ WORK( * ), Z( LDZ, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DTGEXC reorders the generalized real Schur decomposition of a real
|
|
* matrix pair (A,B) using an orthogonal equivalence transformation
|
|
*
|
|
* (A, B) = Q * (A, B) * Z',
|
|
*
|
|
* so that the diagonal block of (A, B) with row index IFST is moved
|
|
* to row ILST.
|
|
*
|
|
* (A, B) must be in generalized real Schur canonical form (as returned
|
|
* by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
|
|
* diagonal blocks. B is upper triangular.
|
|
*
|
|
* Optionally, the matrices Q and Z of generalized Schur vectors are
|
|
* updated.
|
|
*
|
|
* Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
|
|
* Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
|
|
*
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* WANTQ (input) LOGICAL
|
|
* .TRUE. : update the left transformation matrix Q;
|
|
* .FALSE.: do not update Q.
|
|
*
|
|
* WANTZ (input) LOGICAL
|
|
* .TRUE. : update the right transformation matrix Z;
|
|
* .FALSE.: do not update Z.
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrices A and B. N >= 0.
|
|
*
|
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* On entry, the matrix A in generalized real Schur canonical
|
|
* form.
|
|
* On exit, the updated matrix A, again in generalized
|
|
* real Schur canonical form.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
|
|
* On entry, the matrix B in generalized real Schur canonical
|
|
* form (A,B).
|
|
* On exit, the updated matrix B, again in generalized
|
|
* real Schur canonical form (A,B).
|
|
*
|
|
* LDB (input) INTEGER
|
|
* The leading dimension of the array B. LDB >= max(1,N).
|
|
*
|
|
* Q (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
|
|
* On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
|
|
* On exit, the updated matrix Q.
|
|
* If WANTQ = .FALSE., Q is not referenced.
|
|
*
|
|
* LDQ (input) INTEGER
|
|
* The leading dimension of the array Q. LDQ >= 1.
|
|
* If WANTQ = .TRUE., LDQ >= N.
|
|
*
|
|
* Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
|
|
* On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
|
|
* On exit, the updated matrix Z.
|
|
* If WANTZ = .FALSE., Z is not referenced.
|
|
*
|
|
* LDZ (input) INTEGER
|
|
* The leading dimension of the array Z. LDZ >= 1.
|
|
* If WANTZ = .TRUE., LDZ >= N.
|
|
*
|
|
* IFST (input/output) INTEGER
|
|
* ILST (input/output) INTEGER
|
|
* Specify the reordering of the diagonal blocks of (A, B).
|
|
* The block with row index IFST is moved to row ILST, by a
|
|
* sequence of swapping 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, ILST <= N.
|
|
*
|
|
* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
|
|
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
|
|
*
|
|
* LWORK (input) INTEGER
|
|
* The dimension of the array WORK.
|
|
* LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
|
|
*
|
|
* If LWORK = -1, then a workspace query is assumed; the routine
|
|
* only calculates the optimal size of the WORK array, returns
|
|
* this value as the first entry of the WORK array, and no error
|
|
* message related to LWORK is issued by XERBLA.
|
|
*
|
|
* INFO (output) INTEGER
|
|
* =0: successful exit.
|
|
* <0: if INFO = -i, the i-th argument had an illegal value.
|
|
* =1: The transformed matrix pair (A, B) would be too far
|
|
* from generalized Schur form; the problem is ill-
|
|
* conditioned. (A, B) may have been partially reordered,
|
|
* and ILST points to the first row of the current
|
|
* position of the block being moved.
|
|
*
|
|
* Further Details
|
|
* ===============
|
|
*
|
|
* Based on contributions by
|
|
* Bo Kagstrom and Peter Poromaa, Department of Computing Science,
|
|
* Umea University, S-901 87 Umea, Sweden.
|
|
*
|
|
* [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
|
|
* Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
|
|
* M.S. Moonen et al (eds), Linear Algebra for Large Scale and
|
|
* Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ZERO
|
|
PARAMETER ( ZERO = 0.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL LQUERY
|
|
INTEGER HERE, LWMIN, NBF, NBL, NBNEXT
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DTGEX2, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Decode and test input arguments.
|
|
*
|
|
INFO = 0
|
|
LQUERY = ( LWORK.EQ.-1 )
|
|
IF( N.LT.0 ) THEN
|
|
INFO = -3
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -5
|
|
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
|
|
INFO = -7
|
|
ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
|
|
INFO = -9
|
|
ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
|
|
INFO = -11
|
|
ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
|
|
INFO = -12
|
|
ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
|
|
INFO = -13
|
|
END IF
|
|
*
|
|
IF( INFO.EQ.0 ) THEN
|
|
IF( N.LE.1 ) THEN
|
|
LWMIN = 1
|
|
ELSE
|
|
LWMIN = 4*N + 16
|
|
END IF
|
|
WORK(1) = LWMIN
|
|
*
|
|
IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
|
|
INFO = -15
|
|
END IF
|
|
END IF
|
|
*
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DTGEXC', -INFO )
|
|
RETURN
|
|
ELSE IF( LQUERY ) THEN
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.LE.1 )
|
|
$ RETURN
|
|
*
|
|
* Determine the first row of the specified block and find out
|
|
* if it is 1-by-1 or 2-by-2.
|
|
*
|
|
IF( IFST.GT.1 ) THEN
|
|
IF( A( IFST, IFST-1 ).NE.ZERO )
|
|
$ IFST = IFST - 1
|
|
END IF
|
|
NBF = 1
|
|
IF( IFST.LT.N ) THEN
|
|
IF( A( IFST+1, IFST ).NE.ZERO )
|
|
$ NBF = 2
|
|
END IF
|
|
*
|
|
* Determine the first row of the final block
|
|
* and find out if it is 1-by-1 or 2-by-2.
|
|
*
|
|
IF( ILST.GT.1 ) THEN
|
|
IF( A( ILST, ILST-1 ).NE.ZERO )
|
|
$ ILST = ILST - 1
|
|
END IF
|
|
NBL = 1
|
|
IF( ILST.LT.N ) THEN
|
|
IF( A( 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 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( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
|
|
$ NBNEXT = 2
|
|
END IF
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
|
|
$ LDZ, HERE, NBF, NBNEXT, WORK, LWORK, 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( A( 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( A( HERE+3, HERE+2 ).NE.ZERO )
|
|
$ NBNEXT = 2
|
|
END IF
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
|
|
$ LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
ILST = HERE
|
|
RETURN
|
|
END IF
|
|
IF( NBNEXT.EQ.1 ) THEN
|
|
*
|
|
* Swap two 1-by-1 blocks.
|
|
*
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
|
|
$ LDZ, HERE, 1, 1, WORK, LWORK, INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
ILST = HERE
|
|
RETURN
|
|
END IF
|
|
HERE = HERE + 1
|
|
*
|
|
ELSE
|
|
*
|
|
* Recompute NBNEXT in case of 2-by-2 split.
|
|
*
|
|
IF( A( HERE+2, HERE+1 ).EQ.ZERO )
|
|
$ NBNEXT = 1
|
|
IF( NBNEXT.EQ.2 ) THEN
|
|
*
|
|
* 2-by-2 block did not split.
|
|
*
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
|
|
$ Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
|
|
$ INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
ILST = HERE
|
|
RETURN
|
|
END IF
|
|
HERE = HERE + 2
|
|
ELSE
|
|
*
|
|
* 2-by-2 block did split.
|
|
*
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
|
|
$ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
ILST = HERE
|
|
RETURN
|
|
END IF
|
|
HERE = HERE + 1
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
|
|
$ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
ILST = HERE
|
|
RETURN
|
|
END IF
|
|
HERE = HERE + 1
|
|
END IF
|
|
*
|
|
END IF
|
|
END IF
|
|
IF( HERE.LT.ILST )
|
|
$ GO TO 10
|
|
ELSE
|
|
HERE = IFST
|
|
*
|
|
20 CONTINUE
|
|
*
|
|
* Swap 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.GE.3 ) THEN
|
|
IF( A( HERE-1, HERE-2 ).NE.ZERO )
|
|
$ NBNEXT = 2
|
|
END IF
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
|
|
$ LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
|
|
$ 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( A( 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( A( HERE-1, HERE-2 ).NE.ZERO )
|
|
$ NBNEXT = 2
|
|
END IF
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
|
|
$ LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
|
|
$ INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
ILST = HERE
|
|
RETURN
|
|
END IF
|
|
IF( NBNEXT.EQ.1 ) THEN
|
|
*
|
|
* Swap two 1-by-1 blocks.
|
|
*
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
|
|
$ LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
ILST = HERE
|
|
RETURN
|
|
END IF
|
|
HERE = HERE - 1
|
|
ELSE
|
|
*
|
|
* Recompute NBNEXT in case of 2-by-2 split.
|
|
*
|
|
IF( A( HERE, HERE-1 ).EQ.ZERO )
|
|
$ NBNEXT = 1
|
|
IF( NBNEXT.EQ.2 ) THEN
|
|
*
|
|
* 2-by-2 block did not split.
|
|
*
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
|
|
$ Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
ILST = HERE
|
|
RETURN
|
|
END IF
|
|
HERE = HERE - 2
|
|
ELSE
|
|
*
|
|
* 2-by-2 block did split.
|
|
*
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
|
|
$ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
ILST = HERE
|
|
RETURN
|
|
END IF
|
|
HERE = HERE - 1
|
|
CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
|
|
$ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
|
|
IF( INFO.NE.0 ) THEN
|
|
ILST = HERE
|
|
RETURN
|
|
END IF
|
|
HERE = HERE - 1
|
|
END IF
|
|
END IF
|
|
END IF
|
|
IF( HERE.GT.ILST )
|
|
$ GO TO 20
|
|
END IF
|
|
ILST = HERE
|
|
WORK( 1 ) = LWMIN
|
|
RETURN
|
|
*
|
|
* End of DTGEXC
|
|
*
|
|
END
|