LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dtrexc.f
Go to the documentation of this file.
00001 *> \brief \b DTREXC
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DTREXC + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrexc.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrexc.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrexc.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
00022 *                          INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          COMPQ
00026 *       INTEGER            IFST, ILST, INFO, LDQ, LDT, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> DTREXC reorders the real Schur factorization of a real matrix
00039 *> A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
00040 *> moved to row ILST.
00041 *>
00042 *> The real Schur form T is reordered by an orthogonal similarity
00043 *> transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
00044 *> is updated by postmultiplying it with Z.
00045 *>
00046 *> T must be in Schur canonical form (as returned by DHSEQR), that is,
00047 *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
00048 *> 2-by-2 diagonal block has its diagonal elements equal and its
00049 *> off-diagonal elements of opposite sign.
00050 *> \endverbatim
00051 *
00052 *  Arguments:
00053 *  ==========
00054 *
00055 *> \param[in] COMPQ
00056 *> \verbatim
00057 *>          COMPQ is CHARACTER*1
00058 *>          = 'V':  update the matrix Q of Schur vectors;
00059 *>          = 'N':  do not update Q.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] N
00063 *> \verbatim
00064 *>          N is INTEGER
00065 *>          The order of the matrix T. N >= 0.
00066 *> \endverbatim
00067 *>
00068 *> \param[in,out] T
00069 *> \verbatim
00070 *>          T is DOUBLE PRECISION array, dimension (LDT,N)
00071 *>          On entry, the upper quasi-triangular matrix T, in Schur
00072 *>          Schur canonical form.
00073 *>          On exit, the reordered upper quasi-triangular matrix, again
00074 *>          in Schur canonical form.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] LDT
00078 *> \verbatim
00079 *>          LDT is INTEGER
00080 *>          The leading dimension of the array T. LDT >= max(1,N).
00081 *> \endverbatim
00082 *>
00083 *> \param[in,out] Q
00084 *> \verbatim
00085 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
00086 *>          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
00087 *>          On exit, if COMPQ = 'V', Q has been postmultiplied by the
00088 *>          orthogonal transformation matrix Z which reorders T.
00089 *>          If COMPQ = 'N', Q is not referenced.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] LDQ
00093 *> \verbatim
00094 *>          LDQ is INTEGER
00095 *>          The leading dimension of the array Q.  LDQ >= max(1,N).
00096 *> \endverbatim
00097 *>
00098 *> \param[in,out] IFST
00099 *> \verbatim
00100 *>          IFST is INTEGER
00101 *> \endverbatim
00102 *>
00103 *> \param[in,out] ILST
00104 *> \verbatim
00105 *>          ILST is INTEGER
00106 *>
00107 *>          Specify the reordering of the diagonal blocks of T.
00108 *>          The block with row index IFST is moved to row ILST, by a
00109 *>          sequence of transpositions between adjacent blocks.
00110 *>          On exit, if IFST pointed on entry to the second row of a
00111 *>          2-by-2 block, it is changed to point to the first row; ILST
00112 *>          always points to the first row of the block in its final
00113 *>          position (which may differ from its input value by +1 or -1).
00114 *>          1 <= IFST <= N; 1 <= ILST <= N.
00115 *> \endverbatim
00116 *>
00117 *> \param[out] WORK
00118 *> \verbatim
00119 *>          WORK is DOUBLE PRECISION array, dimension (N)
00120 *> \endverbatim
00121 *>
00122 *> \param[out] INFO
00123 *> \verbatim
00124 *>          INFO is INTEGER
00125 *>          = 0:  successful exit
00126 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00127 *>          = 1:  two adjacent blocks were too close to swap (the problem
00128 *>                is very ill-conditioned); T may have been partially
00129 *>                reordered, and ILST points to the first row of the
00130 *>                current position of the block being moved.
00131 *> \endverbatim
00132 *
00133 *  Authors:
00134 *  ========
00135 *
00136 *> \author Univ. of Tennessee 
00137 *> \author Univ. of California Berkeley 
00138 *> \author Univ. of Colorado Denver 
00139 *> \author NAG Ltd. 
00140 *
00141 *> \date November 2011
00142 *
00143 *> \ingroup doubleOTHERcomputational
00144 *
00145 *  =====================================================================
00146       SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
00147      $                   INFO )
00148 *
00149 *  -- LAPACK computational routine (version 3.4.0) --
00150 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00151 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00152 *     November 2011
00153 *
00154 *     .. Scalar Arguments ..
00155       CHARACTER          COMPQ
00156       INTEGER            IFST, ILST, INFO, LDQ, LDT, N
00157 *     ..
00158 *     .. Array Arguments ..
00159       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
00160 *     ..
00161 *
00162 *  =====================================================================
00163 *
00164 *     .. Parameters ..
00165       DOUBLE PRECISION   ZERO
00166       PARAMETER          ( ZERO = 0.0D+0 )
00167 *     ..
00168 *     .. Local Scalars ..
00169       LOGICAL            WANTQ
00170       INTEGER            HERE, NBF, NBL, NBNEXT
00171 *     ..
00172 *     .. External Functions ..
00173       LOGICAL            LSAME
00174       EXTERNAL           LSAME
00175 *     ..
00176 *     .. External Subroutines ..
00177       EXTERNAL           DLAEXC, XERBLA
00178 *     ..
00179 *     .. Intrinsic Functions ..
00180       INTRINSIC          MAX
00181 *     ..
00182 *     .. Executable Statements ..
00183 *
00184 *     Decode and test the input arguments.
00185 *
00186       INFO = 0
00187       WANTQ = LSAME( COMPQ, 'V' )
00188       IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
00189          INFO = -1
00190       ELSE IF( N.LT.0 ) THEN
00191          INFO = -2
00192       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00193          INFO = -4
00194       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
00195          INFO = -6
00196       ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
00197          INFO = -7
00198       ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
00199          INFO = -8
00200       END IF
00201       IF( INFO.NE.0 ) THEN
00202          CALL XERBLA( 'DTREXC', -INFO )
00203          RETURN
00204       END IF
00205 *
00206 *     Quick return if possible
00207 *
00208       IF( N.LE.1 )
00209      $   RETURN
00210 *
00211 *     Determine the first row of specified block
00212 *     and find out it is 1 by 1 or 2 by 2.
00213 *
00214       IF( IFST.GT.1 ) THEN
00215          IF( T( IFST, IFST-1 ).NE.ZERO )
00216      $      IFST = IFST - 1
00217       END IF
00218       NBF = 1
00219       IF( IFST.LT.N ) THEN
00220          IF( T( IFST+1, IFST ).NE.ZERO )
00221      $      NBF = 2
00222       END IF
00223 *
00224 *     Determine the first row of the final block
00225 *     and find out it is 1 by 1 or 2 by 2.
00226 *
00227       IF( ILST.GT.1 ) THEN
00228          IF( T( ILST, ILST-1 ).NE.ZERO )
00229      $      ILST = ILST - 1
00230       END IF
00231       NBL = 1
00232       IF( ILST.LT.N ) THEN
00233          IF( T( ILST+1, ILST ).NE.ZERO )
00234      $      NBL = 2
00235       END IF
00236 *
00237       IF( IFST.EQ.ILST )
00238      $   RETURN
00239 *
00240       IF( IFST.LT.ILST ) THEN
00241 *
00242 *        Update ILST
00243 *
00244          IF( NBF.EQ.2 .AND. NBL.EQ.1 )
00245      $      ILST = ILST - 1
00246          IF( NBF.EQ.1 .AND. NBL.EQ.2 )
00247      $      ILST = ILST + 1
00248 *
00249          HERE = IFST
00250 *
00251    10    CONTINUE
00252 *
00253 *        Swap block with next one below
00254 *
00255          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
00256 *
00257 *           Current block either 1 by 1 or 2 by 2
00258 *
00259             NBNEXT = 1
00260             IF( HERE+NBF+1.LE.N ) THEN
00261                IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
00262      $            NBNEXT = 2
00263             END IF
00264             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
00265      $                   WORK, INFO )
00266             IF( INFO.NE.0 ) THEN
00267                ILST = HERE
00268                RETURN
00269             END IF
00270             HERE = HERE + NBNEXT
00271 *
00272 *           Test if 2 by 2 block breaks into two 1 by 1 blocks
00273 *
00274             IF( NBF.EQ.2 ) THEN
00275                IF( T( HERE+1, HERE ).EQ.ZERO )
00276      $            NBF = 3
00277             END IF
00278 *
00279          ELSE
00280 *
00281 *           Current block consists of two 1 by 1 blocks each of which
00282 *           must be swapped individually
00283 *
00284             NBNEXT = 1
00285             IF( HERE+3.LE.N ) THEN
00286                IF( T( HERE+3, HERE+2 ).NE.ZERO )
00287      $            NBNEXT = 2
00288             END IF
00289             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
00290      $                   WORK, INFO )
00291             IF( INFO.NE.0 ) THEN
00292                ILST = HERE
00293                RETURN
00294             END IF
00295             IF( NBNEXT.EQ.1 ) THEN
00296 *
00297 *              Swap two 1 by 1 blocks, no problems possible
00298 *
00299                CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
00300      $                      WORK, INFO )
00301                HERE = HERE + 1
00302             ELSE
00303 *
00304 *              Recompute NBNEXT in case 2 by 2 split
00305 *
00306                IF( T( HERE+2, HERE+1 ).EQ.ZERO )
00307      $            NBNEXT = 1
00308                IF( NBNEXT.EQ.2 ) THEN
00309 *
00310 *                 2 by 2 Block did not split
00311 *
00312                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
00313      $                         NBNEXT, WORK, INFO )
00314                   IF( INFO.NE.0 ) THEN
00315                      ILST = HERE
00316                      RETURN
00317                   END IF
00318                   HERE = HERE + 2
00319                ELSE
00320 *
00321 *                 2 by 2 Block did split
00322 *
00323                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
00324      $                         WORK, INFO )
00325                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
00326      $                         WORK, INFO )
00327                   HERE = HERE + 2
00328                END IF
00329             END IF
00330          END IF
00331          IF( HERE.LT.ILST )
00332      $      GO TO 10
00333 *
00334       ELSE
00335 *
00336          HERE = IFST
00337    20    CONTINUE
00338 *
00339 *        Swap block with next one above
00340 *
00341          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
00342 *
00343 *           Current block either 1 by 1 or 2 by 2
00344 *
00345             NBNEXT = 1
00346             IF( HERE.GE.3 ) THEN
00347                IF( T( HERE-1, HERE-2 ).NE.ZERO )
00348      $            NBNEXT = 2
00349             END IF
00350             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
00351      $                   NBF, WORK, INFO )
00352             IF( INFO.NE.0 ) THEN
00353                ILST = HERE
00354                RETURN
00355             END IF
00356             HERE = HERE - NBNEXT
00357 *
00358 *           Test if 2 by 2 block breaks into two 1 by 1 blocks
00359 *
00360             IF( NBF.EQ.2 ) THEN
00361                IF( T( HERE+1, HERE ).EQ.ZERO )
00362      $            NBF = 3
00363             END IF
00364 *
00365          ELSE
00366 *
00367 *           Current block consists of two 1 by 1 blocks each of which
00368 *           must be swapped individually
00369 *
00370             NBNEXT = 1
00371             IF( HERE.GE.3 ) THEN
00372                IF( T( HERE-1, HERE-2 ).NE.ZERO )
00373      $            NBNEXT = 2
00374             END IF
00375             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
00376      $                   1, WORK, INFO )
00377             IF( INFO.NE.0 ) THEN
00378                ILST = HERE
00379                RETURN
00380             END IF
00381             IF( NBNEXT.EQ.1 ) THEN
00382 *
00383 *              Swap two 1 by 1 blocks, no problems possible
00384 *
00385                CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
00386      $                      WORK, INFO )
00387                HERE = HERE - 1
00388             ELSE
00389 *
00390 *              Recompute NBNEXT in case 2 by 2 split
00391 *
00392                IF( T( HERE, HERE-1 ).EQ.ZERO )
00393      $            NBNEXT = 1
00394                IF( NBNEXT.EQ.2 ) THEN
00395 *
00396 *                 2 by 2 Block did not split
00397 *
00398                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
00399      $                         WORK, INFO )
00400                   IF( INFO.NE.0 ) THEN
00401                      ILST = HERE
00402                      RETURN
00403                   END IF
00404                   HERE = HERE - 2
00405                ELSE
00406 *
00407 *                 2 by 2 Block did split
00408 *
00409                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
00410      $                         WORK, INFO )
00411                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
00412      $                         WORK, INFO )
00413                   HERE = HERE - 2
00414                END IF
00415             END IF
00416          END IF
00417          IF( HERE.GT.ILST )
00418      $      GO TO 20
00419       END IF
00420       ILST = HERE
00421 *
00422       RETURN
00423 *
00424 *     End of DTREXC
00425 *
00426       END
 All Files Functions