LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
stgexc.f
Go to the documentation of this file.
00001 *> \brief \b STGEXC
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download STGEXC + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgexc.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgexc.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgexc.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE STGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00022 *                          LDZ, IFST, ILST, WORK, LWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       LOGICAL            WANTQ, WANTZ
00026 *       INTEGER            IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       REAL               A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00030 *      $                   WORK( * ), Z( LDZ, * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> STGEXC reorders the generalized real Schur decomposition of a real
00040 *> matrix pair (A,B) using an orthogonal equivalence transformation
00041 *>
00042 *>                (A, B) = Q * (A, B) * Z**T,
00043 *>
00044 *> so that the diagonal block of (A, B) with row index IFST is moved
00045 *> to row ILST.
00046 *>
00047 *> (A, B) must be in generalized real Schur canonical form (as returned
00048 *> by SGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
00049 *> diagonal blocks. B is upper triangular.
00050 *>
00051 *> Optionally, the matrices Q and Z of generalized Schur vectors are
00052 *> updated.
00053 *>
00054 *>        Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
00055 *>        Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
00056 *>
00057 *> \endverbatim
00058 *
00059 *  Arguments:
00060 *  ==========
00061 *
00062 *> \param[in] WANTQ
00063 *> \verbatim
00064 *>          WANTQ is LOGICAL
00065 *>          .TRUE. : update the left transformation matrix Q;
00066 *>          .FALSE.: do not update Q.
00067 *> \endverbatim
00068 *>
00069 *> \param[in] WANTZ
00070 *> \verbatim
00071 *>          WANTZ is LOGICAL
00072 *>          .TRUE. : update the right transformation matrix Z;
00073 *>          .FALSE.: do not update Z.
00074 *> \endverbatim
00075 *>
00076 *> \param[in] N
00077 *> \verbatim
00078 *>          N is INTEGER
00079 *>          The order of the matrices A and B. N >= 0.
00080 *> \endverbatim
00081 *>
00082 *> \param[in,out] A
00083 *> \verbatim
00084 *>          A is REAL array, dimension (LDA,N)
00085 *>          On entry, the matrix A in generalized real Schur canonical
00086 *>          form.
00087 *>          On exit, the updated matrix A, again in generalized
00088 *>          real Schur canonical form.
00089 *> \endverbatim
00090 *>
00091 *> \param[in] LDA
00092 *> \verbatim
00093 *>          LDA is INTEGER
00094 *>          The leading dimension of the array A. LDA >= max(1,N).
00095 *> \endverbatim
00096 *>
00097 *> \param[in,out] B
00098 *> \verbatim
00099 *>          B is REAL array, dimension (LDB,N)
00100 *>          On entry, the matrix B in generalized real Schur canonical
00101 *>          form (A,B).
00102 *>          On exit, the updated matrix B, again in generalized
00103 *>          real Schur canonical form (A,B).
00104 *> \endverbatim
00105 *>
00106 *> \param[in] LDB
00107 *> \verbatim
00108 *>          LDB is INTEGER
00109 *>          The leading dimension of the array B. LDB >= max(1,N).
00110 *> \endverbatim
00111 *>
00112 *> \param[in,out] Q
00113 *> \verbatim
00114 *>          Q is REAL array, dimension (LDZ,N)
00115 *>          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
00116 *>          On exit, the updated matrix Q.
00117 *>          If WANTQ = .FALSE., Q is not referenced.
00118 *> \endverbatim
00119 *>
00120 *> \param[in] LDQ
00121 *> \verbatim
00122 *>          LDQ is INTEGER
00123 *>          The leading dimension of the array Q. LDQ >= 1.
00124 *>          If WANTQ = .TRUE., LDQ >= N.
00125 *> \endverbatim
00126 *>
00127 *> \param[in,out] Z
00128 *> \verbatim
00129 *>          Z is REAL array, dimension (LDZ,N)
00130 *>          On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
00131 *>          On exit, the updated matrix Z.
00132 *>          If WANTZ = .FALSE., Z is not referenced.
00133 *> \endverbatim
00134 *>
00135 *> \param[in] LDZ
00136 *> \verbatim
00137 *>          LDZ is INTEGER
00138 *>          The leading dimension of the array Z. LDZ >= 1.
00139 *>          If WANTZ = .TRUE., LDZ >= N.
00140 *> \endverbatim
00141 *>
00142 *> \param[in,out] IFST
00143 *> \verbatim
00144 *>          IFST is INTEGER
00145 *> \endverbatim
00146 *>
00147 *> \param[in,out] ILST
00148 *> \verbatim
00149 *>          ILST is INTEGER
00150 *>          Specify the reordering of the diagonal blocks of (A, B).
00151 *>          The block with row index IFST is moved to row ILST, by a
00152 *>          sequence of swapping between adjacent blocks.
00153 *>          On exit, if IFST pointed on entry to the second row of
00154 *>          a 2-by-2 block, it is changed to point to the first row;
00155 *>          ILST always points to the first row of the block in its
00156 *>          final position (which may differ from its input value by
00157 *>          +1 or -1). 1 <= IFST, ILST <= N.
00158 *> \endverbatim
00159 *>
00160 *> \param[out] WORK
00161 *> \verbatim
00162 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00163 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00164 *> \endverbatim
00165 *>
00166 *> \param[in] LWORK
00167 *> \verbatim
00168 *>          LWORK is INTEGER
00169 *>          The dimension of the array WORK.
00170 *>          LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
00171 *>
00172 *>          If LWORK = -1, then a workspace query is assumed; the routine
00173 *>          only calculates the optimal size of the WORK array, returns
00174 *>          this value as the first entry of the WORK array, and no error
00175 *>          message related to LWORK is issued by XERBLA.
00176 *> \endverbatim
00177 *>
00178 *> \param[out] INFO
00179 *> \verbatim
00180 *>          INFO is INTEGER
00181 *>           =0:  successful exit.
00182 *>           <0:  if INFO = -i, the i-th argument had an illegal value.
00183 *>           =1:  The transformed matrix pair (A, B) would be too far
00184 *>                from generalized Schur form; the problem is ill-
00185 *>                conditioned. (A, B) may have been partially reordered,
00186 *>                and ILST points to the first row of the current
00187 *>                position of the block being moved.
00188 *> \endverbatim
00189 *
00190 *  Authors:
00191 *  ========
00192 *
00193 *> \author Univ. of Tennessee 
00194 *> \author Univ. of California Berkeley 
00195 *> \author Univ. of Colorado Denver 
00196 *> \author NAG Ltd. 
00197 *
00198 *> \date November 2011
00199 *
00200 *> \ingroup realGEcomputational
00201 *
00202 *> \par Contributors:
00203 *  ==================
00204 *>
00205 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00206 *>     Umea University, S-901 87 Umea, Sweden.
00207 *
00208 *> \par References:
00209 *  ================
00210 *>
00211 *> \verbatim
00212 *>
00213 *>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
00214 *>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
00215 *>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
00216 *>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
00217 *> \endverbatim
00218 *>
00219 *  =====================================================================
00220       SUBROUTINE STGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00221      $                   LDZ, IFST, ILST, WORK, LWORK, INFO )
00222 *
00223 *  -- LAPACK computational routine (version 3.4.0) --
00224 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00225 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00226 *     November 2011
00227 *
00228 *     .. Scalar Arguments ..
00229       LOGICAL            WANTQ, WANTZ
00230       INTEGER            IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
00231 *     ..
00232 *     .. Array Arguments ..
00233       REAL               A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00234      $                   WORK( * ), Z( LDZ, * )
00235 *     ..
00236 *
00237 *  =====================================================================
00238 *
00239 *     .. Parameters ..
00240       REAL               ZERO
00241       PARAMETER          ( ZERO = 0.0E+0 )
00242 *     ..
00243 *     .. Local Scalars ..
00244       LOGICAL            LQUERY
00245       INTEGER            HERE, LWMIN, NBF, NBL, NBNEXT
00246 *     ..
00247 *     .. External Subroutines ..
00248       EXTERNAL           STGEX2, XERBLA
00249 *     ..
00250 *     .. Intrinsic Functions ..
00251       INTRINSIC          MAX
00252 *     ..
00253 *     .. Executable Statements ..
00254 *
00255 *     Decode and test input arguments.
00256 *
00257       INFO = 0
00258       LQUERY = ( LWORK.EQ.-1 )
00259       IF( N.LT.0 ) THEN
00260          INFO = -3
00261       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00262          INFO = -5
00263       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00264          INFO = -7
00265       ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
00266          INFO = -9
00267       ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
00268          INFO = -11
00269       ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
00270          INFO = -12
00271       ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
00272          INFO = -13
00273       END IF
00274 *
00275       IF( INFO.EQ.0 ) THEN
00276          IF( N.LE.1 ) THEN
00277             LWMIN = 1
00278          ELSE
00279             LWMIN = 4*N + 16
00280          END IF
00281          WORK(1) = LWMIN
00282 *
00283          IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
00284             INFO = -15
00285          END IF
00286       END IF
00287 *
00288       IF( INFO.NE.0 ) THEN
00289          CALL XERBLA( 'STGEXC', -INFO )
00290          RETURN
00291       ELSE IF( LQUERY ) THEN
00292          RETURN
00293       END IF
00294 *
00295 *     Quick return if possible
00296 *
00297       IF( N.LE.1 )
00298      $   RETURN
00299 *
00300 *     Determine the first row of the specified block and find out
00301 *     if it is 1-by-1 or 2-by-2.
00302 *
00303       IF( IFST.GT.1 ) THEN
00304          IF( A( IFST, IFST-1 ).NE.ZERO )
00305      $      IFST = IFST - 1
00306       END IF
00307       NBF = 1
00308       IF( IFST.LT.N ) THEN
00309          IF( A( IFST+1, IFST ).NE.ZERO )
00310      $      NBF = 2
00311       END IF
00312 *
00313 *     Determine the first row of the final block
00314 *     and find out if it is 1-by-1 or 2-by-2.
00315 *
00316       IF( ILST.GT.1 ) THEN
00317          IF( A( ILST, ILST-1 ).NE.ZERO )
00318      $      ILST = ILST - 1
00319       END IF
00320       NBL = 1
00321       IF( ILST.LT.N ) THEN
00322          IF( A( ILST+1, ILST ).NE.ZERO )
00323      $      NBL = 2
00324       END IF
00325       IF( IFST.EQ.ILST )
00326      $   RETURN
00327 *
00328       IF( IFST.LT.ILST ) THEN
00329 *
00330 *        Update ILST.
00331 *
00332          IF( NBF.EQ.2 .AND. NBL.EQ.1 )
00333      $      ILST = ILST - 1
00334          IF( NBF.EQ.1 .AND. NBL.EQ.2 )
00335      $      ILST = ILST + 1
00336 *
00337          HERE = IFST
00338 *
00339    10    CONTINUE
00340 *
00341 *        Swap with next one below.
00342 *
00343          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
00344 *
00345 *           Current block either 1-by-1 or 2-by-2.
00346 *
00347             NBNEXT = 1
00348             IF( HERE+NBF+1.LE.N ) THEN
00349                IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
00350      $            NBNEXT = 2
00351             END IF
00352             CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00353      $                   LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO )
00354             IF( INFO.NE.0 ) THEN
00355                ILST = HERE
00356                RETURN
00357             END IF
00358             HERE = HERE + NBNEXT
00359 *
00360 *           Test if 2-by-2 block breaks into two 1-by-1 blocks.
00361 *
00362             IF( NBF.EQ.2 ) THEN
00363                IF( A( HERE+1, HERE ).EQ.ZERO )
00364      $            NBF = 3
00365             END IF
00366 *
00367          ELSE
00368 *
00369 *           Current block consists of two 1-by-1 blocks, each of which
00370 *           must be swapped individually.
00371 *
00372             NBNEXT = 1
00373             IF( HERE+3.LE.N ) THEN
00374                IF( A( HERE+3, HERE+2 ).NE.ZERO )
00375      $            NBNEXT = 2
00376             END IF
00377             CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00378      $                   LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
00379             IF( INFO.NE.0 ) THEN
00380                ILST = HERE
00381                RETURN
00382             END IF
00383             IF( NBNEXT.EQ.1 ) THEN
00384 *
00385 *              Swap two 1-by-1 blocks.
00386 *
00387                CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00388      $                      LDZ, HERE, 1, 1, WORK, LWORK, INFO )
00389                IF( INFO.NE.0 ) THEN
00390                   ILST = HERE
00391                   RETURN
00392                END IF
00393                HERE = HERE + 1
00394 *
00395             ELSE
00396 *
00397 *              Recompute NBNEXT in case of 2-by-2 split.
00398 *
00399                IF( A( HERE+2, HERE+1 ).EQ.ZERO )
00400      $            NBNEXT = 1
00401                IF( NBNEXT.EQ.2 ) THEN
00402 *
00403 *                 2-by-2 block did not split.
00404 *
00405                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00406      $                         Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
00407      $                         INFO )
00408                   IF( INFO.NE.0 ) THEN
00409                      ILST = HERE
00410                      RETURN
00411                   END IF
00412                   HERE = HERE + 2
00413                ELSE
00414 *
00415 *                 2-by-2 block did split.
00416 *
00417                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00418      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
00419                   IF( INFO.NE.0 ) THEN
00420                      ILST = HERE
00421                      RETURN
00422                   END IF
00423                   HERE = HERE + 1
00424                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00425      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
00426                   IF( INFO.NE.0 ) THEN
00427                      ILST = HERE
00428                      RETURN
00429                   END IF
00430                   HERE = HERE + 1
00431                END IF
00432 *
00433             END IF
00434          END IF
00435          IF( HERE.LT.ILST )
00436      $      GO TO 10
00437       ELSE
00438          HERE = IFST
00439 *
00440    20    CONTINUE
00441 *
00442 *        Swap with next one below.
00443 *
00444          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
00445 *
00446 *           Current block either 1-by-1 or 2-by-2.
00447 *
00448             NBNEXT = 1
00449             IF( HERE.GE.3 ) THEN
00450                IF( A( HERE-1, HERE-2 ).NE.ZERO )
00451      $            NBNEXT = 2
00452             END IF
00453             CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00454      $                   LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
00455      $                   INFO )
00456             IF( INFO.NE.0 ) THEN
00457                ILST = HERE
00458                RETURN
00459             END IF
00460             HERE = HERE - NBNEXT
00461 *
00462 *           Test if 2-by-2 block breaks into two 1-by-1 blocks.
00463 *
00464             IF( NBF.EQ.2 ) THEN
00465                IF( A( HERE+1, HERE ).EQ.ZERO )
00466      $            NBF = 3
00467             END IF
00468 *
00469          ELSE
00470 *
00471 *           Current block consists of two 1-by-1 blocks, each of which
00472 *           must be swapped individually.
00473 *
00474             NBNEXT = 1
00475             IF( HERE.GE.3 ) THEN
00476                IF( A( HERE-1, HERE-2 ).NE.ZERO )
00477      $            NBNEXT = 2
00478             END IF
00479             CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00480      $                   LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
00481      $                   INFO )
00482             IF( INFO.NE.0 ) THEN
00483                ILST = HERE
00484                RETURN
00485             END IF
00486             IF( NBNEXT.EQ.1 ) THEN
00487 *
00488 *              Swap two 1-by-1 blocks.
00489 *
00490                CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00491      $                      LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
00492                IF( INFO.NE.0 ) THEN
00493                   ILST = HERE
00494                   RETURN
00495                END IF
00496                HERE = HERE - 1
00497             ELSE
00498 *
00499 *             Recompute NBNEXT in case of 2-by-2 split.
00500 *
00501                IF( A( HERE, HERE-1 ).EQ.ZERO )
00502      $            NBNEXT = 1
00503                IF( NBNEXT.EQ.2 ) THEN
00504 *
00505 *                 2-by-2 block did not split.
00506 *
00507                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00508      $                         Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
00509                   IF( INFO.NE.0 ) THEN
00510                      ILST = HERE
00511                      RETURN
00512                   END IF
00513                   HERE = HERE - 2
00514                ELSE
00515 *
00516 *                 2-by-2 block did split.
00517 *
00518                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00519      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
00520                   IF( INFO.NE.0 ) THEN
00521                      ILST = HERE
00522                      RETURN
00523                   END IF
00524                   HERE = HERE - 1
00525                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00526      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
00527                   IF( INFO.NE.0 ) THEN
00528                      ILST = HERE
00529                      RETURN
00530                   END IF
00531                   HERE = HERE - 1
00532                END IF
00533             END IF
00534          END IF
00535          IF( HERE.GT.ILST )
00536      $      GO TO 20
00537       END IF
00538       ILST = HERE
00539       WORK( 1 ) = LWMIN
00540       RETURN
00541 *
00542 *     End of STGEXC
00543 *
00544       END
 All Files Functions