LAPACK: SRC/stgexc.f Source File

Go to the documentation of this file.
00001       SUBROUTINE STGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00002      $                   LDZ, IFST, ILST, WORK, LWORK, 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       LOGICAL            WANTQ, WANTZ
00011       INTEGER            IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
00015      $                   WORK( * ), Z( LDZ, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  STGEXC reorders the generalized real Schur decomposition of a real
00022 *  matrix pair (A,B) using an orthogonal equivalence transformation
00023 *
00024 *                 (A, B) = Q * (A, B) * Z',
00025 *
00026 *  so that the diagonal block of (A, B) with row index IFST is moved
00027 *  to row ILST.
00028 *
00029 *  (A, B) must be in generalized real Schur canonical form (as returned
00030 *  by SGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
00031 *  diagonal blocks. B is upper triangular.
00032 *
00033 *  Optionally, the matrices Q and Z of generalized Schur vectors are
00034 *  updated.
00035 *
00036 *         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
00037 *         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
00038 *
00039 *
00040 *  Arguments
00041 *  =========
00042 *
00043 *  WANTQ   (input) LOGICAL
00044 *          .TRUE. : update the left transformation matrix Q;
00045 *          .FALSE.: do not update Q.
00046 *
00047 *  WANTZ   (input) LOGICAL
00048 *          .TRUE. : update the right transformation matrix Z;
00049 *          .FALSE.: do not update Z.
00050 *
00051 *  N       (input) INTEGER
00052 *          The order of the matrices A and B. N >= 0.
00053 *
00054 *  A       (input/output) REAL array, dimension (LDA,N)
00055 *          On entry, the matrix A in generalized real Schur canonical
00056 *          form.
00057 *          On exit, the updated matrix A, again in generalized
00058 *          real Schur canonical form.
00059 *
00060 *  LDA     (input)  INTEGER
00061 *          The leading dimension of the array A. LDA >= max(1,N).
00062 *
00063 *  B       (input/output) REAL array, dimension (LDB,N)
00064 *          On entry, the matrix B in generalized real Schur canonical
00065 *          form (A,B).
00066 *          On exit, the updated matrix B, again in generalized
00067 *          real Schur canonical form (A,B).
00068 *
00069 *  LDB     (input)  INTEGER
00070 *          The leading dimension of the array B. LDB >= max(1,N).
00071 *
00072 *  Q       (input/output) REAL array, dimension (LDZ,N)
00073 *          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
00074 *          On exit, the updated matrix Q.
00075 *          If WANTQ = .FALSE., Q is not referenced.
00076 *
00077 *  LDQ     (input) INTEGER
00078 *          The leading dimension of the array Q. LDQ >= 1.
00079 *          If WANTQ = .TRUE., LDQ >= N.
00080 *
00081 *  Z       (input/output) REAL array, dimension (LDZ,N)
00082 *          On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
00083 *          On exit, the updated matrix Z.
00084 *          If WANTZ = .FALSE., Z is not referenced.
00085 *
00086 *  LDZ     (input) INTEGER
00087 *          The leading dimension of the array Z. LDZ >= 1.
00088 *          If WANTZ = .TRUE., LDZ >= N.
00089 *
00090 *  IFST    (input/output) INTEGER
00091 *  ILST    (input/output) INTEGER
00092 *          Specify the reordering of the diagonal blocks of (A, B).
00093 *          The block with row index IFST is moved to row ILST, by a
00094 *          sequence of swapping between adjacent blocks.
00095 *          On exit, if IFST pointed on entry to the second row of
00096 *          a 2-by-2 block, it is changed to point to the first row;
00097 *          ILST always points to the first row of the block in its
00098 *          final position (which may differ from its input value by
00099 *          +1 or -1). 1 <= IFST, ILST <= N.
00100 *
00101 *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
00102 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00103 *
00104 *  LWORK   (input) INTEGER
00105 *          The dimension of the array WORK.
00106 *          LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
00107 *
00108 *          If LWORK = -1, then a workspace query is assumed; the routine
00109 *          only calculates the optimal size of the WORK array, returns
00110 *          this value as the first entry of the WORK array, and no error
00111 *          message related to LWORK is issued by XERBLA.
00112 *
00113 *  INFO    (output) INTEGER
00114 *           =0:  successful exit.
00115 *           <0:  if INFO = -i, the i-th argument had an illegal value.
00116 *           =1:  The transformed matrix pair (A, B) would be too far
00117 *                from generalized Schur form; the problem is ill-
00118 *                conditioned. (A, B) may have been partially reordered,
00119 *                and ILST points to the first row of the current
00120 *                position of the block being moved.
00121 *
00122 *  Further Details
00123 *  ===============
00124 *
00125 *  Based on contributions by
00126 *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00127 *     Umea University, S-901 87 Umea, Sweden.
00128 *
00129 *  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
00130 *      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
00131 *      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
00132 *      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
00133 *
00134 *  =====================================================================
00135 *
00136 *     .. Parameters ..
00137       REAL               ZERO
00138       PARAMETER          ( ZERO = 0.0E+0 )
00139 *     ..
00140 *     .. Local Scalars ..
00141       LOGICAL            LQUERY
00142       INTEGER            HERE, LWMIN, NBF, NBL, NBNEXT
00143 *     ..
00144 *     .. External Subroutines ..
00145       EXTERNAL           STGEX2, XERBLA
00146 *     ..
00147 *     .. Intrinsic Functions ..
00148       INTRINSIC          MAX
00149 *     ..
00150 *     .. Executable Statements ..
00151 *
00152 *     Decode and test input arguments.
00153 *
00154       INFO = 0
00155       LQUERY = ( LWORK.EQ.-1 )
00156       IF( N.LT.0 ) THEN
00157          INFO = -3
00158       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00159          INFO = -5
00160       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00161          INFO = -7
00162       ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
00163          INFO = -9
00164       ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
00165          INFO = -11
00166       ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
00167          INFO = -12
00168       ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
00169          INFO = -13
00170       END IF
00171 *
00172       IF( INFO.EQ.0 ) THEN
00173          IF( N.LE.1 ) THEN
00174             LWMIN = 1
00175          ELSE
00176             LWMIN = 4*N + 16
00177          END IF
00178          WORK(1) = LWMIN
00179 *
00180          IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
00181             INFO = -15
00182          END IF
00183       END IF
00184 *
00185       IF( INFO.NE.0 ) THEN
00186          CALL XERBLA( 'STGEXC', -INFO )
00187          RETURN
00188       ELSE IF( LQUERY ) THEN
00189          RETURN
00190       END IF
00191 *
00192 *     Quick return if possible
00193 *
00194       IF( N.LE.1 )
00195      $   RETURN
00196 *
00197 *     Determine the first row of the specified block and find out
00198 *     if it is 1-by-1 or 2-by-2.
00199 *
00200       IF( IFST.GT.1 ) THEN
00201          IF( A( IFST, IFST-1 ).NE.ZERO )
00202      $      IFST = IFST - 1
00203       END IF
00204       NBF = 1
00205       IF( IFST.LT.N ) THEN
00206          IF( A( IFST+1, IFST ).NE.ZERO )
00207      $      NBF = 2
00208       END IF
00209 *
00210 *     Determine the first row of the final block
00211 *     and find out if it is 1-by-1 or 2-by-2.
00212 *
00213       IF( ILST.GT.1 ) THEN
00214          IF( A( ILST, ILST-1 ).NE.ZERO )
00215      $      ILST = ILST - 1
00216       END IF
00217       NBL = 1
00218       IF( ILST.LT.N ) THEN
00219          IF( A( ILST+1, ILST ).NE.ZERO )
00220      $      NBL = 2
00221       END IF
00222       IF( IFST.EQ.ILST )
00223      $   RETURN
00224 *
00225       IF( IFST.LT.ILST ) THEN
00226 *
00227 *        Update ILST.
00228 *
00229          IF( NBF.EQ.2 .AND. NBL.EQ.1 )
00230      $      ILST = ILST - 1
00231          IF( NBF.EQ.1 .AND. NBL.EQ.2 )
00232      $      ILST = ILST + 1
00233 *
00234          HERE = IFST
00235 *
00236    10    CONTINUE
00237 *
00238 *        Swap with next one below.
00239 *
00240          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
00241 *
00242 *           Current block either 1-by-1 or 2-by-2.
00243 *
00244             NBNEXT = 1
00245             IF( HERE+NBF+1.LE.N ) THEN
00246                IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
00247      $            NBNEXT = 2
00248             END IF
00249             CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00250      $                   LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO )
00251             IF( INFO.NE.0 ) THEN
00252                ILST = HERE
00253                RETURN
00254             END IF
00255             HERE = HERE + NBNEXT
00256 *
00257 *           Test if 2-by-2 block breaks into two 1-by-1 blocks.
00258 *
00259             IF( NBF.EQ.2 ) THEN
00260                IF( A( HERE+1, HERE ).EQ.ZERO )
00261      $            NBF = 3
00262             END IF
00263 *
00264          ELSE
00265 *
00266 *           Current block consists of two 1-by-1 blocks, each of which
00267 *           must be swapped individually.
00268 *
00269             NBNEXT = 1
00270             IF( HERE+3.LE.N ) THEN
00271                IF( A( HERE+3, HERE+2 ).NE.ZERO )
00272      $            NBNEXT = 2
00273             END IF
00274             CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00275      $                   LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
00276             IF( INFO.NE.0 ) THEN
00277                ILST = HERE
00278                RETURN
00279             END IF
00280             IF( NBNEXT.EQ.1 ) THEN
00281 *
00282 *              Swap two 1-by-1 blocks.
00283 *
00284                CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00285      $                      LDZ, HERE, 1, 1, WORK, LWORK, INFO )
00286                IF( INFO.NE.0 ) THEN
00287                   ILST = HERE
00288                   RETURN
00289                END IF
00290                HERE = HERE + 1
00291 *
00292             ELSE
00293 *
00294 *              Recompute NBNEXT in case of 2-by-2 split.
00295 *
00296                IF( A( HERE+2, HERE+1 ).EQ.ZERO )
00297      $            NBNEXT = 1
00298                IF( NBNEXT.EQ.2 ) THEN
00299 *
00300 *                 2-by-2 block did not split.
00301 *
00302                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00303      $                         Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
00304      $                         INFO )
00305                   IF( INFO.NE.0 ) THEN
00306                      ILST = HERE
00307                      RETURN
00308                   END IF
00309                   HERE = HERE + 2
00310                ELSE
00311 *
00312 *                 2-by-2 block did split.
00313 *
00314                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00315      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
00316                   IF( INFO.NE.0 ) THEN
00317                      ILST = HERE
00318                      RETURN
00319                   END IF
00320                   HERE = HERE + 1
00321                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00322      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
00323                   IF( INFO.NE.0 ) THEN
00324                      ILST = HERE
00325                      RETURN
00326                   END IF
00327                   HERE = HERE + 1
00328                END IF
00329 *
00330             END IF
00331          END IF
00332          IF( HERE.LT.ILST )
00333      $      GO TO 10
00334       ELSE
00335          HERE = IFST
00336 *
00337    20    CONTINUE
00338 *
00339 *        Swap with next one below.
00340 *
00341          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
00342 *
00343 *           Current block either 1-by-1 or 2-by-2.
00344 *
00345             NBNEXT = 1
00346             IF( HERE.GE.3 ) THEN
00347                IF( A( HERE-1, HERE-2 ).NE.ZERO )
00348      $            NBNEXT = 2
00349             END IF
00350             CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00351      $                   LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
00352      $                   INFO )
00353             IF( INFO.NE.0 ) THEN
00354                ILST = HERE
00355                RETURN
00356             END IF
00357             HERE = HERE - NBNEXT
00358 *
00359 *           Test if 2-by-2 block breaks into two 1-by-1 blocks.
00360 *
00361             IF( NBF.EQ.2 ) THEN
00362                IF( A( HERE+1, HERE ).EQ.ZERO )
00363      $            NBF = 3
00364             END IF
00365 *
00366          ELSE
00367 *
00368 *           Current block consists of two 1-by-1 blocks, each of which
00369 *           must be swapped individually.
00370 *
00371             NBNEXT = 1
00372             IF( HERE.GE.3 ) THEN
00373                IF( A( HERE-1, HERE-2 ).NE.ZERO )
00374      $            NBNEXT = 2
00375             END IF
00376             CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
00377      $                   LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
00378      $                   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, NBNEXT, 1, WORK, LWORK, INFO )
00389                IF( INFO.NE.0 ) THEN
00390                   ILST = HERE
00391                   RETURN
00392                END IF
00393                HERE = HERE - 1
00394             ELSE
00395 *
00396 *             Recompute NBNEXT in case of 2-by-2 split.
00397 *
00398                IF( A( HERE, HERE-1 ).EQ.ZERO )
00399      $            NBNEXT = 1
00400                IF( NBNEXT.EQ.2 ) THEN
00401 *
00402 *                 2-by-2 block did not split.
00403 *
00404                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00405      $                         Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
00406                   IF( INFO.NE.0 ) THEN
00407                      ILST = HERE
00408                      RETURN
00409                   END IF
00410                   HERE = HERE - 2
00411                ELSE
00412 *
00413 *                 2-by-2 block did split.
00414 *
00415                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00416      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
00417                   IF( INFO.NE.0 ) THEN
00418                      ILST = HERE
00419                      RETURN
00420                   END IF
00421                   HERE = HERE - 1
00422                   CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
00423      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
00424                   IF( INFO.NE.0 ) THEN
00425                      ILST = HERE
00426                      RETURN
00427                   END IF
00428                   HERE = HERE - 1
00429                END IF
00430             END IF
00431          END IF
00432          IF( HERE.GT.ILST )
00433      $      GO TO 20
00434       END IF
00435       ILST = HERE
00436       WORK( 1 ) = LWMIN
00437       RETURN
00438 *
00439 *     End of STGEXC
00440 *
00441       END