SRC/strexc.f

   1 *> \brief \b STREXC
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download STREXC + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strexc.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strexc.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strexc.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE STREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
  22 *                          INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          COMPQ
  26 *       INTEGER            IFST, ILST, INFO, LDQ, LDT, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       REAL               Q( LDQ, * ), T( LDT, * ), WORK( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> STREXC reorders the real Schur factorization of a real matrix
  39 *> A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
  40 *> moved to row ILST.
  41 *>
  42 *> The real Schur form T is reordered by an orthogonal similarity
  43 *> transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
  44 *> is updated by postmultiplying it with Z.
  45 *>
  46 *> T must be in Schur canonical form (as returned by SHSEQR), that is,
  47 *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
  48 *> 2-by-2 diagonal block has its diagonal elements equal and its
  49 *> off-diagonal elements of opposite sign.
  50 *> \endverbatim
  51 *
  52 *  Arguments:
  53 *  ==========
  54 *
  55 *> \param[in] COMPQ
  56 *> \verbatim
  57 *>          COMPQ is CHARACTER*1
  58 *>          = 'V':  update the matrix Q of Schur vectors;
  59 *>          = 'N':  do not update Q.
  60 *> \endverbatim
  61 *>
  62 *> \param[in] N
  63 *> \verbatim
  64 *>          N is INTEGER
  65 *>          The order of the matrix T. N >= 0.
  66 *>          If N == 0 arguments ILST and IFST may be any value.
  67 *> \endverbatim
  68 *>
  69 *> \param[in,out] T
  70 *> \verbatim
  71 *>          T is REAL array, dimension (LDT,N)
  72 *>          On entry, the upper quasi-triangular matrix T, in Schur
  73 *>          Schur canonical form.
  74 *>          On exit, the reordered upper quasi-triangular matrix, again
  75 *>          in Schur canonical form.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] LDT
  79 *> \verbatim
  80 *>          LDT is INTEGER
  81 *>          The leading dimension of the array T. LDT >= max(1,N).
  82 *> \endverbatim
  83 *>
  84 *> \param[in,out] Q
  85 *> \verbatim
  86 *>          Q is REAL array, dimension (LDQ,N)
  87 *>          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
  88 *>          On exit, if COMPQ = 'V', Q has been postmultiplied by the
  89 *>          orthogonal transformation matrix Z which reorders T.
  90 *>          If COMPQ = 'N', Q is not referenced.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] LDQ
  94 *> \verbatim
  95 *>          LDQ is INTEGER
  96 *>          The leading dimension of the array Q.  LDQ >= 1, and if
  97 *>          COMPQ = 'V', LDQ >= max(1,N).
  98 *> \endverbatim
  99 *>
 100 *> \param[in,out] IFST
 101 *> \verbatim
 102 *>          IFST is INTEGER
 103 *> \endverbatim
 104 *>
 105 *> \param[in,out] ILST
 106 *> \verbatim
 107 *>          ILST is INTEGER
 108 *>
 109 *>          Specify the reordering of the diagonal blocks of T.
 110 *>          The block with row index IFST is moved to row ILST, by a
 111 *>          sequence of transpositions between adjacent blocks.
 112 *>          On exit, if IFST pointed on entry to the second row of a
 113 *>          2-by-2 block, it is changed to point to the first row; ILST
 114 *>          always points to the first row of the block in its final
 115 *>          position (which may differ from its input value by +1 or -1).
 116 *>          1 <= IFST <= N; 1 <= ILST <= N.
 117 *> \endverbatim
 118 *>
 119 *> \param[out] WORK
 120 *> \verbatim
 121 *>          WORK is REAL array, dimension (N)
 122 *> \endverbatim
 123 *>
 124 *> \param[out] INFO
 125 *> \verbatim
 126 *>          INFO is INTEGER
 127 *>          = 0:  successful exit
 128 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 129 *>          = 1:  two adjacent blocks were too close to swap (the problem
 130 *>                is very ill-conditioned); T may have been partially
 131 *>                reordered, and ILST points to the first row of the
 132 *>                current position of the block being moved.
 133 *> \endverbatim
 134 *
 135 *  Authors:
 136 *  ========
 137 *
 138 *> \author Univ. of Tennessee
 139 *> \author Univ. of California Berkeley
 140 *> \author Univ. of Colorado Denver
 141 *> \author NAG Ltd.
 142 *
 143 *> \date November 2011
 144 *
 145 *> \ingroup realOTHERcomputational
 146 *
 147 *  =====================================================================
 148       SUBROUTINE STREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
 149      $                   INFO )
 150 *
 151 *  -- LAPACK computational routine (version 3.4.0) --
 152 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 153 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 154 *     November 2011
 155 *
 156 *     .. Scalar Arguments ..
 157       CHARACTER          COMPQ
 158       INTEGER            IFST, ILST, INFO, LDQ, LDT, N
 159 *     ..
 160 *     .. Array Arguments ..
 161       REAL               Q( LDQ, * ), T( LDT, * ), WORK( * )
 162 *     ..
 163 *
 164 *  =====================================================================
 165 *
 166 *     .. Parameters ..
 167       REAL               ZERO
 168       PARAMETER          ( ZERO = 0.0E+0 )
 169 *     ..
 170 *     .. Local Scalars ..
 171       LOGICAL            WANTQ
 172       INTEGER            HERE, NBF, NBL, NBNEXT
 173 *     ..
 174 *     .. External Functions ..
 175       LOGICAL            LSAME
 176       EXTERNAL           LSAME
 177 *     ..
 178 *     .. External Subroutines ..
 179       EXTERNAL           SLAEXC, XERBLA
 180 *     ..
 181 *     .. Intrinsic Functions ..
 182       INTRINSIC          MAX
 183 *     ..
 184 *     .. Executable Statements ..
 185 *
 186 *     Decode and test the input arguments.
 187 *
 188       INFO = 0
 189       WANTQ = LSAME( COMPQ, 'V' )
 190       IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
 191          INFO = -1
 192       ELSE IF( N.LT.0 ) THEN
 193          INFO = -2
 194       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
 195          INFO = -4
 196       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
 197          INFO = -6
 198       ELSE IF(( IFST.LT.1 .OR. IFST.GT.N ).AND.( N.GT.0 )) THEN
 199          INFO = -7
 200       ELSE IF(( ILST.LT.1 .OR. ILST.GT.N ).AND.( N.GT.0 )) THEN
 201          INFO = -8
 202       END IF
 203       IF( INFO.NE.0 ) THEN
 204          CALL XERBLA( 'STREXC', -INFO )
 205          RETURN
 206       END IF
 207 *
 208 *     Quick return if possible
 209 *
 210       IF( N.LE.1 )
 211      $   RETURN
 212 *
 213 *     Determine the first row of specified block
 214 *     and find out it is 1 by 1 or 2 by 2.
 215 *
 216       IF( IFST.GT.1 ) THEN
 217          IF( T( IFST, IFST-1 ).NE.ZERO )
 218      $      IFST = IFST - 1
 219       END IF
 220       NBF = 1
 221       IF( IFST.LT.N ) THEN
 222          IF( T( IFST+1, IFST ).NE.ZERO )
 223      $      NBF = 2
 224       END IF
 225 *
 226 *     Determine the first row of the final block
 227 *     and find out it is 1 by 1 or 2 by 2.
 228 *
 229       IF( ILST.GT.1 ) THEN
 230          IF( T( ILST, ILST-1 ).NE.ZERO )
 231      $      ILST = ILST - 1
 232       END IF
 233       NBL = 1
 234       IF( ILST.LT.N ) THEN
 235          IF( T( ILST+1, ILST ).NE.ZERO )
 236      $      NBL = 2
 237       END IF
 238 *
 239       IF( IFST.EQ.ILST )
 240      $   RETURN
 241 *
 242       IF( IFST.LT.ILST ) THEN
 243 *
 244 *        Update ILST
 245 *
 246          IF( NBF.EQ.2 .AND. NBL.EQ.1 )
 247      $      ILST = ILST - 1
 248          IF( NBF.EQ.1 .AND. NBL.EQ.2 )
 249      $      ILST = ILST + 1
 250 *
 251          HERE = IFST
 252 *
 253    10    CONTINUE
 254 *
 255 *        Swap block with next one below
 256 *
 257          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
 258 *
 259 *           Current block either 1 by 1 or 2 by 2
 260 *
 261             NBNEXT = 1
 262             IF( HERE+NBF+1.LE.N ) THEN
 263                IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
 264      $            NBNEXT = 2
 265             END IF
 266             CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
 267      $                   WORK, INFO )
 268             IF( INFO.NE.0 ) THEN
 269                ILST = HERE
 270                RETURN
 271             END IF
 272             HERE = HERE + NBNEXT
 273 *
 274 *           Test if 2 by 2 block breaks into two 1 by 1 blocks
 275 *
 276             IF( NBF.EQ.2 ) THEN
 277                IF( T( HERE+1, HERE ).EQ.ZERO )
 278      $            NBF = 3
 279             END IF
 280 *
 281          ELSE
 282 *
 283 *           Current block consists of two 1 by 1 blocks each of which
 284 *           must be swapped individually
 285 *
 286             NBNEXT = 1
 287             IF( HERE+3.LE.N ) THEN
 288                IF( T( HERE+3, HERE+2 ).NE.ZERO )
 289      $            NBNEXT = 2
 290             END IF
 291             CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
 292      $                   WORK, INFO )
 293             IF( INFO.NE.0 ) THEN
 294                ILST = HERE
 295                RETURN
 296             END IF
 297             IF( NBNEXT.EQ.1 ) THEN
 298 *
 299 *              Swap two 1 by 1 blocks, no problems possible
 300 *
 301                CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
 302      $                      WORK, INFO )
 303                HERE = HERE + 1
 304             ELSE
 305 *
 306 *              Recompute NBNEXT in case 2 by 2 split
 307 *
 308                IF( T( HERE+2, HERE+1 ).EQ.ZERO )
 309      $            NBNEXT = 1
 310                IF( NBNEXT.EQ.2 ) THEN
 311 *
 312 *                 2 by 2 Block did not split
 313 *
 314                   CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
 315      $                         NBNEXT, WORK, INFO )
 316                   IF( INFO.NE.0 ) THEN
 317                      ILST = HERE
 318                      RETURN
 319                   END IF
 320                   HERE = HERE + 2
 321                ELSE
 322 *
 323 *                 2 by 2 Block did split
 324 *
 325                   CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
 326      $                         WORK, INFO )
 327                   CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
 328      $                         WORK, INFO )
 329                   HERE = HERE + 2
 330                END IF
 331             END IF
 332          END IF
 333          IF( HERE.LT.ILST )
 334      $      GO TO 10
 335 *
 336       ELSE
 337 *
 338          HERE = IFST
 339    20    CONTINUE
 340 *
 341 *        Swap block with next one above
 342 *
 343          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
 344 *
 345 *           Current block either 1 by 1 or 2 by 2
 346 *
 347             NBNEXT = 1
 348             IF( HERE.GE.3 ) THEN
 349                IF( T( HERE-1, HERE-2 ).NE.ZERO )
 350      $            NBNEXT = 2
 351             END IF
 352             CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
 353      $                   NBF, WORK, INFO )
 354             IF( INFO.NE.0 ) THEN
 355                ILST = HERE
 356                RETURN
 357             END IF
 358             HERE = HERE - NBNEXT
 359 *
 360 *           Test if 2 by 2 block breaks into two 1 by 1 blocks
 361 *
 362             IF( NBF.EQ.2 ) THEN
 363                IF( T( HERE+1, HERE ).EQ.ZERO )
 364      $            NBF = 3
 365             END IF
 366 *
 367          ELSE
 368 *
 369 *           Current block consists of two 1 by 1 blocks each of which
 370 *           must be swapped individually
 371 *
 372             NBNEXT = 1
 373             IF( HERE.GE.3 ) THEN
 374                IF( T( HERE-1, HERE-2 ).NE.ZERO )
 375      $            NBNEXT = 2
 376             END IF
 377             CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
 378      $                   1, WORK, INFO )
 379             IF( INFO.NE.0 ) THEN
 380                ILST = HERE
 381                RETURN
 382             END IF
 383             IF( NBNEXT.EQ.1 ) THEN
 384 *
 385 *              Swap two 1 by 1 blocks, no problems possible
 386 *
 387                CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
 388      $                      WORK, INFO )
 389                HERE = HERE - 1
 390             ELSE
 391 *
 392 *              Recompute NBNEXT in case 2 by 2 split
 393 *
 394                IF( T( HERE, HERE-1 ).EQ.ZERO )
 395      $            NBNEXT = 1
 396                IF( NBNEXT.EQ.2 ) THEN
 397 *
 398 *                 2 by 2 Block did not split
 399 *
 400                   CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
 401      $                         WORK, INFO )
 402                   IF( INFO.NE.0 ) THEN
 403                      ILST = HERE
 404                      RETURN
 405                   END IF
 406                   HERE = HERE - 2
 407                ELSE
 408 *
 409 *                 2 by 2 Block did split
 410 *
 411                   CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
 412      $                         WORK, INFO )
 413                   CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
 414      $                         WORK, INFO )
 415                   HERE = HERE - 2
 416                END IF
 417             END IF
 418          END IF
 419          IF( HERE.GT.ILST )
 420      $      GO TO 20
 421       END IF
 422       ILST = HERE
 423 *
 424       RETURN
 425 *
 426 *     End of STREXC
 427 *
 428       END