SRC/dtgexc.f

   1 *> \brief \b DTGEXC
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DTGEXC + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgexc.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgexc.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgexc.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
  22 *                          LDZ, IFST, ILST, WORK, LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       LOGICAL            WANTQ, WANTZ
  26 *       INTEGER            IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
  30 *      $                   WORK( * ), Z( LDZ, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DTGEXC reorders the generalized real Schur decomposition of a real
  40 *> matrix pair (A,B) using an orthogonal equivalence transformation
  41 *>
  42 *>                (A, B) = Q * (A, B) * Z**T,
  43 *>
  44 *> so that the diagonal block of (A, B) with row index IFST is moved
  45 *> to row ILST.
  46 *>
  47 *> (A, B) must be in generalized real Schur canonical form (as returned
  48 *> by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
  49 *> diagonal blocks. B is upper triangular.
  50 *>
  51 *> Optionally, the matrices Q and Z of generalized Schur vectors are
  52 *> updated.
  53 *>
  54 *>        Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
  55 *>        Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
  56 *>
  57 *> \endverbatim
  58 *
  59 *  Arguments:
  60 *  ==========
  61 *
  62 *> \param[in] WANTQ
  63 *> \verbatim
  64 *>          WANTQ is LOGICAL
  65 *>          .TRUE. : update the left transformation matrix Q;
  66 *>          .FALSE.: do not update Q.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] WANTZ
  70 *> \verbatim
  71 *>          WANTZ is LOGICAL
  72 *>          .TRUE. : update the right transformation matrix Z;
  73 *>          .FALSE.: do not update Z.
  74 *> \endverbatim
  75 *>
  76 *> \param[in] N
  77 *> \verbatim
  78 *>          N is INTEGER
  79 *>          The order of the matrices A and B. N >= 0.
  80 *> \endverbatim
  81 *>
  82 *> \param[in,out] A
  83 *> \verbatim
  84 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  85 *>          On entry, the matrix A in generalized real Schur canonical
  86 *>          form.
  87 *>          On exit, the updated matrix A, again in generalized
  88 *>          real Schur canonical form.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] LDA
  92 *> \verbatim
  93 *>          LDA is INTEGER
  94 *>          The leading dimension of the array A. LDA >= max(1,N).
  95 *> \endverbatim
  96 *>
  97 *> \param[in,out] B
  98 *> \verbatim
  99 *>          B is DOUBLE PRECISION array, dimension (LDB,N)
 100 *>          On entry, the matrix B in generalized real Schur canonical
 101 *>          form (A,B).
 102 *>          On exit, the updated matrix B, again in generalized
 103 *>          real Schur canonical form (A,B).
 104 *> \endverbatim
 105 *>
 106 *> \param[in] LDB
 107 *> \verbatim
 108 *>          LDB is INTEGER
 109 *>          The leading dimension of the array B. LDB >= max(1,N).
 110 *> \endverbatim
 111 *>
 112 *> \param[in,out] Q
 113 *> \verbatim
 114 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
 115 *>          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
 116 *>          On exit, the updated matrix Q.
 117 *>          If WANTQ = .FALSE., Q is not referenced.
 118 *> \endverbatim
 119 *>
 120 *> \param[in] LDQ
 121 *> \verbatim
 122 *>          LDQ is INTEGER
 123 *>          The leading dimension of the array Q. LDQ >= 1.
 124 *>          If WANTQ = .TRUE., LDQ >= N.
 125 *> \endverbatim
 126 *>
 127 *> \param[in,out] Z
 128 *> \verbatim
 129 *>          Z is DOUBLE PRECISION array, dimension (LDZ,N)
 130 *>          On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
 131 *>          On exit, the updated matrix Z.
 132 *>          If WANTZ = .FALSE., Z is not referenced.
 133 *> \endverbatim
 134 *>
 135 *> \param[in] LDZ
 136 *> \verbatim
 137 *>          LDZ is INTEGER
 138 *>          The leading dimension of the array Z. LDZ >= 1.
 139 *>          If WANTZ = .TRUE., LDZ >= N.
 140 *> \endverbatim
 141 *>
 142 *> \param[in,out] IFST
 143 *> \verbatim
 144 *>          IFST is INTEGER
 145 *> \endverbatim
 146 *>
 147 *> \param[in,out] ILST
 148 *> \verbatim
 149 *>          ILST is INTEGER
 150 *>          Specify the reordering of the diagonal blocks of (A, B).
 151 *>          The block with row index IFST is moved to row ILST, by a
 152 *>          sequence of swapping between adjacent blocks.
 153 *>          On exit, if IFST pointed on entry to the second row of
 154 *>          a 2-by-2 block, it is changed to point to the first row;
 155 *>          ILST always points to the first row of the block in its
 156 *>          final position (which may differ from its input value by
 157 *>          +1 or -1). 1 <= IFST, ILST <= N.
 158 *> \endverbatim
 159 *>
 160 *> \param[out] WORK
 161 *> \verbatim
 162 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 163 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 164 *> \endverbatim
 165 *>
 166 *> \param[in] LWORK
 167 *> \verbatim
 168 *>          LWORK is INTEGER
 169 *>          The dimension of the array WORK.
 170 *>          LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
 171 *>
 172 *>          If LWORK = -1, then a workspace query is assumed; the routine
 173 *>          only calculates the optimal size of the WORK array, returns
 174 *>          this value as the first entry of the WORK array, and no error
 175 *>          message related to LWORK is issued by XERBLA.
 176 *> \endverbatim
 177 *>
 178 *> \param[out] INFO
 179 *> \verbatim
 180 *>          INFO is INTEGER
 181 *>           =0:  successful exit.
 182 *>           <0:  if INFO = -i, the i-th argument had an illegal value.
 183 *>           =1:  The transformed matrix pair (A, B) would be too far
 184 *>                from generalized Schur form; the problem is ill-
 185 *>                conditioned. (A, B) may have been partially reordered,
 186 *>                and ILST points to the first row of the current
 187 *>                position of the block being moved.
 188 *> \endverbatim
 189 *
 190 *  Authors:
 191 *  ========
 192 *
 193 *> \author Univ. of Tennessee
 194 *> \author Univ. of California Berkeley
 195 *> \author Univ. of Colorado Denver
 196 *> \author NAG Ltd.
 197 *
 198 *> \date November 2011
 199 *
 200 *> \ingroup doubleGEcomputational
 201 *
 202 *> \par Contributors:
 203 *  ==================
 204 *>
 205 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 206 *>     Umea University, S-901 87 Umea, Sweden.
 207 *
 208 *> \par References:
 209 *  ================
 210 *>
 211 *> \verbatim
 212 *>
 213 *>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
 214 *>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
 215 *>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
 216 *>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
 217 *> \endverbatim
 218 *>
 219 *  =====================================================================
 220       SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
 221      $                   LDZ, IFST, ILST, WORK, LWORK, INFO )
 222 *
 223 *  -- LAPACK computational routine (version 3.4.0) --
 224 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 225 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 226 *     November 2011
 227 *
 228 *     .. Scalar Arguments ..
 229       LOGICAL            WANTQ, WANTZ
 230       INTEGER            IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
 231 *     ..
 232 *     .. Array Arguments ..
 233       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
 234      $                   WORK( * ), Z( LDZ, * )
 235 *     ..
 236 *
 237 *  =====================================================================
 238 *
 239 *     .. Parameters ..
 240       DOUBLE PRECISION   ZERO
 241       PARAMETER          ( ZERO = 0.0D+0 )
 242 *     ..
 243 *     .. Local Scalars ..
 244       LOGICAL            LQUERY
 245       INTEGER            HERE, LWMIN, NBF, NBL, NBNEXT
 246 *     ..
 247 *     .. External Subroutines ..
 248       EXTERNAL           DTGEX2, XERBLA
 249 *     ..
 250 *     .. Intrinsic Functions ..
 251       INTRINSIC          MAX
 252 *     ..
 253 *     .. Executable Statements ..
 254 *
 255 *     Decode and test input arguments.
 256 *
 257       INFO = 0
 258       LQUERY = ( LWORK.EQ.-1 )
 259       IF( N.LT.0 ) THEN
 260          INFO = -3
 261       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 262          INFO = -5
 263       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 264          INFO = -7
 265       ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
 266          INFO = -9
 267       ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
 268          INFO = -11
 269       ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
 270          INFO = -12
 271       ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
 272          INFO = -13
 273       END IF
 274 *
 275       IF( INFO.EQ.0 ) THEN
 276          IF( N.LE.1 ) THEN
 277             LWMIN = 1
 278          ELSE
 279             LWMIN = 4*N + 16
 280          END IF
 281          WORK(1) = LWMIN
 282 *
 283          IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
 284             INFO = -15
 285          END IF
 286       END IF
 287 *
 288       IF( INFO.NE.0 ) THEN
 289          CALL XERBLA( 'DTGEXC', -INFO )
 290          RETURN
 291       ELSE IF( LQUERY ) THEN
 292          RETURN
 293       END IF
 294 *
 295 *     Quick return if possible
 296 *
 297       IF( N.LE.1 )
 298      $   RETURN
 299 *
 300 *     Determine the first row of the specified block and find out
 301 *     if it is 1-by-1 or 2-by-2.
 302 *
 303       IF( IFST.GT.1 ) THEN
 304          IF( A( IFST, IFST-1 ).NE.ZERO )
 305      $      IFST = IFST - 1
 306       END IF
 307       NBF = 1
 308       IF( IFST.LT.N ) THEN
 309          IF( A( IFST+1, IFST ).NE.ZERO )
 310      $      NBF = 2
 311       END IF
 312 *
 313 *     Determine the first row of the final block
 314 *     and find out if it is 1-by-1 or 2-by-2.
 315 *
 316       IF( ILST.GT.1 ) THEN
 317          IF( A( ILST, ILST-1 ).NE.ZERO )
 318      $      ILST = ILST - 1
 319       END IF
 320       NBL = 1
 321       IF( ILST.LT.N ) THEN
 322          IF( A( ILST+1, ILST ).NE.ZERO )
 323      $      NBL = 2
 324       END IF
 325       IF( IFST.EQ.ILST )
 326      $   RETURN
 327 *
 328       IF( IFST.LT.ILST ) THEN
 329 *
 330 *        Update ILST.
 331 *
 332          IF( NBF.EQ.2 .AND. NBL.EQ.1 )
 333      $      ILST = ILST - 1
 334          IF( NBF.EQ.1 .AND. NBL.EQ.2 )
 335      $      ILST = ILST + 1
 336 *
 337          HERE = IFST
 338 *
 339    10    CONTINUE
 340 *
 341 *        Swap with next one below.
 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+NBF+1.LE.N ) THEN
 349                IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
 350      $            NBNEXT = 2
 351             END IF
 352             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
 353      $                   LDZ, HERE, NBF, NBNEXT, WORK, LWORK, 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( A( 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+3.LE.N ) THEN
 374                IF( A( HERE+3, HERE+2 ).NE.ZERO )
 375      $            NBNEXT = 2
 376             END IF
 377             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
 378      $                   LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, 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.
 386 *
 387                CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
 388      $                      LDZ, HERE, 1, 1, WORK, LWORK, INFO )
 389                IF( INFO.NE.0 ) THEN
 390                   ILST = HERE
 391                   RETURN
 392                END IF
 393                HERE = HERE + 1
 394 *
 395             ELSE
 396 *
 397 *              Recompute NBNEXT in case of 2-by-2 split.
 398 *
 399                IF( A( HERE+2, HERE+1 ).EQ.ZERO )
 400      $            NBNEXT = 1
 401                IF( NBNEXT.EQ.2 ) THEN
 402 *
 403 *                 2-by-2 block did not split.
 404 *
 405                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
 406      $                         Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
 407      $                         INFO )
 408                   IF( INFO.NE.0 ) THEN
 409                      ILST = HERE
 410                      RETURN
 411                   END IF
 412                   HERE = HERE + 2
 413                ELSE
 414 *
 415 *                 2-by-2 block did split.
 416 *
 417                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
 418      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
 419                   IF( INFO.NE.0 ) THEN
 420                      ILST = HERE
 421                      RETURN
 422                   END IF
 423                   HERE = HERE + 1
 424                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
 425      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
 426                   IF( INFO.NE.0 ) THEN
 427                      ILST = HERE
 428                      RETURN
 429                   END IF
 430                   HERE = HERE + 1
 431                END IF
 432 *
 433             END IF
 434          END IF
 435          IF( HERE.LT.ILST )
 436      $      GO TO 10
 437       ELSE
 438          HERE = IFST
 439 *
 440    20    CONTINUE
 441 *
 442 *        Swap with next one below.
 443 *
 444          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
 445 *
 446 *           Current block either 1-by-1 or 2-by-2.
 447 *
 448             NBNEXT = 1
 449             IF( HERE.GE.3 ) THEN
 450                IF( A( HERE-1, HERE-2 ).NE.ZERO )
 451      $            NBNEXT = 2
 452             END IF
 453             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
 454      $                   LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
 455      $                   INFO )
 456             IF( INFO.NE.0 ) THEN
 457                ILST = HERE
 458                RETURN
 459             END IF
 460             HERE = HERE - NBNEXT
 461 *
 462 *           Test if 2-by-2 block breaks into two 1-by-1 blocks.
 463 *
 464             IF( NBF.EQ.2 ) THEN
 465                IF( A( HERE+1, HERE ).EQ.ZERO )
 466      $            NBF = 3
 467             END IF
 468 *
 469          ELSE
 470 *
 471 *           Current block consists of two 1-by-1 blocks, each of which
 472 *           must be swapped individually.
 473 *
 474             NBNEXT = 1
 475             IF( HERE.GE.3 ) THEN
 476                IF( A( HERE-1, HERE-2 ).NE.ZERO )
 477      $            NBNEXT = 2
 478             END IF
 479             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
 480      $                   LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
 481      $                   INFO )
 482             IF( INFO.NE.0 ) THEN
 483                ILST = HERE
 484                RETURN
 485             END IF
 486             IF( NBNEXT.EQ.1 ) THEN
 487 *
 488 *              Swap two 1-by-1 blocks.
 489 *
 490                CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
 491      $                      LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
 492                IF( INFO.NE.0 ) THEN
 493                   ILST = HERE
 494                   RETURN
 495                END IF
 496                HERE = HERE - 1
 497             ELSE
 498 *
 499 *             Recompute NBNEXT in case of 2-by-2 split.
 500 *
 501                IF( A( HERE, HERE-1 ).EQ.ZERO )
 502      $            NBNEXT = 1
 503                IF( NBNEXT.EQ.2 ) THEN
 504 *
 505 *                 2-by-2 block did not split.
 506 *
 507                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
 508      $                         Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
 509                   IF( INFO.NE.0 ) THEN
 510                      ILST = HERE
 511                      RETURN
 512                   END IF
 513                   HERE = HERE - 2
 514                ELSE
 515 *
 516 *                 2-by-2 block did split.
 517 *
 518                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
 519      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
 520                   IF( INFO.NE.0 ) THEN
 521                      ILST = HERE
 522                      RETURN
 523                   END IF
 524                   HERE = HERE - 1
 525                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
 526      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
 527                   IF( INFO.NE.0 ) THEN
 528                      ILST = HERE
 529                      RETURN
 530                   END IF
 531                   HERE = HERE - 1
 532                END IF
 533             END IF
 534          END IF
 535          IF( HERE.GT.ILST )
 536      $      GO TO 20
 537       END IF
 538       ILST = HERE
 539       WORK( 1 ) = LWMIN
 540       RETURN
 541 *
 542 *     End of DTGEXC
 543 *
 544       END