dtrexc.f

Go to the documentation of this file.

*> \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.
*
*> \ingroup trexc
*
*  =====================================================================
      SUBROUTINE dtrexc( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
     $                   INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. 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