SRC/dtrexc.f

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