#include "ESMF_LapackBlas.inc" !> \brief \b DTREXC ! ! =========== DOCUMENTATION =========== ! ! Online html documentation available at ! http://www.netlib.org/lapack/explore-html/ ! !> \htmlonly !> Download DTREXC + dependencies !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrexc.f"> !> [TGZ]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrexc.f"> !> [ZIP]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrexc.f"> !> [TXT]</a> !> \endhtmlonly ! ! Definition: ! =========== ! ! SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, ! INFO ) ! ! .. Scalar Arguments .. ! CHARACTER COMPQ ! INTEGER IFST, ILST, INFO, LDQ, LDT, N ! .. ! .. Array Arguments .. ! DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * ) ! .. ! ! !> \par Purpose: ! ============= !> !> \verbatim !> !> 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. !> \endverbatim ! ! Arguments: ! ========== ! !> \param[in] COMPQ !> \verbatim !> COMPQ is CHARACTER*1 !> = 'V': update the matrix Q of Schur vectors; !> = 'N': do not update Q. !> \endverbatim !> !> \param[in] N !> \verbatim !> N is INTEGER !> The order of the matrix T. N >= 0. !> If N == 0 arguments ILST and IFST may be any value. !> \endverbatim !> !> \param[in,out] T !> \verbatim !> T is 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. !> \endverbatim !> !> \param[in] LDT !> \verbatim !> LDT is INTEGER !> The leading dimension of the array T. LDT >= max(1,N). !> \endverbatim !> !> \param[in,out] Q !> \verbatim !> Q is 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. !> \endverbatim !> !> \param[in] LDQ !> \verbatim !> LDQ is INTEGER !> The leading dimension of the array Q. LDQ >= 1, and if !> COMPQ = 'V', LDQ >= max(1,N). !> \endverbatim !> !> \param[in,out] IFST !> \verbatim !> IFST is INTEGER !> \endverbatim !> !> \param[in,out] ILST !> \verbatim !> ILST is 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. !> \endverbatim !> !> \param[out] WORK !> \verbatim !> WORK is DOUBLE PRECISION array, dimension (N) !> \endverbatim !> !> \param[out] INFO !> \verbatim !> INFO is 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. !> \endverbatim ! ! Authors: ! ======== ! !> \author Univ. of Tennessee !> \author Univ. of California Berkeley !> \author Univ. of Colorado Denver !> \author NAG Ltd. ! !> \date December 2016 ! !> \ingroup doubleOTHERcomputational ! ! ===================================================================== SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, & INFO ) ! ! -- LAPACK computational routine (version 3.7.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! December 2016 ! ! .. Scalar Arguments .. CHARACTER COMPQ INTEGER IFST, ILST, INFO, LDQ, LDT, N ! .. ! .. Array Arguments .. DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * ) ! .. ! ! ===================================================================== ! ! .. 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 ).AND.( N.GT.0 )) THEN INFO = -7 ELSE IF(( ILST.LT.1 .OR. ILST.GT.N ).AND.( N.GT.0 )) 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