LAPACK: SRC/dtrexc.f Source File

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