LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
strsen.f
Go to the documentation of this file.
00001 *> \brief \b STRSEN
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download STRSEN + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strsen.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strsen.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strsen.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE STRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
00022 *                          M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          COMPQ, JOB
00026 *       INTEGER            INFO, LDQ, LDT, LIWORK, LWORK, M, N
00027 *       REAL               S, SEP
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       LOGICAL            SELECT( * )
00031 *       INTEGER            IWORK( * )
00032 *       REAL               Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
00033 *      $                   WR( * )
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> STRSEN reorders the real Schur factorization of a real matrix
00043 *> A = Q*T*Q**T, so that a selected cluster of eigenvalues appears in
00044 *> the leading diagonal blocks of the upper quasi-triangular matrix T,
00045 *> and the leading columns of Q form an orthonormal basis of the
00046 *> corresponding right invariant subspace.
00047 *>
00048 *> Optionally the routine computes the reciprocal condition numbers of
00049 *> the cluster of eigenvalues and/or the invariant subspace.
00050 *>
00051 *> T must be in Schur canonical form (as returned by SHSEQR), that is,
00052 *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
00053 *> 2-by-2 diagonal block has its diagonal elements equal and its
00054 *> off-diagonal elements of opposite sign.
00055 *> \endverbatim
00056 *
00057 *  Arguments:
00058 *  ==========
00059 *
00060 *> \param[in] JOB
00061 *> \verbatim
00062 *>          JOB is CHARACTER*1
00063 *>          Specifies whether condition numbers are required for the
00064 *>          cluster of eigenvalues (S) or the invariant subspace (SEP):
00065 *>          = 'N': none;
00066 *>          = 'E': for eigenvalues only (S);
00067 *>          = 'V': for invariant subspace only (SEP);
00068 *>          = 'B': for both eigenvalues and invariant subspace (S and
00069 *>                 SEP).
00070 *> \endverbatim
00071 *>
00072 *> \param[in] COMPQ
00073 *> \verbatim
00074 *>          COMPQ is CHARACTER*1
00075 *>          = 'V': update the matrix Q of Schur vectors;
00076 *>          = 'N': do not update Q.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] SELECT
00080 *> \verbatim
00081 *>          SELECT is LOGICAL array, dimension (N)
00082 *>          SELECT specifies the eigenvalues in the selected cluster. To
00083 *>          select a real eigenvalue w(j), SELECT(j) must be set to
00084 *>          .TRUE.. To select a complex conjugate pair of eigenvalues
00085 *>          w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
00086 *>          either SELECT(j) or SELECT(j+1) or both must be set to
00087 *>          .TRUE.; a complex conjugate pair of eigenvalues must be
00088 *>          either both included in the cluster or both excluded.
00089 *> \endverbatim
00090 *>
00091 *> \param[in] N
00092 *> \verbatim
00093 *>          N is INTEGER
00094 *>          The order of the matrix T. N >= 0.
00095 *> \endverbatim
00096 *>
00097 *> \param[in,out] T
00098 *> \verbatim
00099 *>          T is REAL array, dimension (LDT,N)
00100 *>          On entry, the upper quasi-triangular matrix T, in Schur
00101 *>          canonical form.
00102 *>          On exit, T is overwritten by the reordered matrix T, again in
00103 *>          Schur canonical form, with the selected eigenvalues in the
00104 *>          leading diagonal blocks.
00105 *> \endverbatim
00106 *>
00107 *> \param[in] LDT
00108 *> \verbatim
00109 *>          LDT is INTEGER
00110 *>          The leading dimension of the array T. LDT >= max(1,N).
00111 *> \endverbatim
00112 *>
00113 *> \param[in,out] Q
00114 *> \verbatim
00115 *>          Q is REAL array, dimension (LDQ,N)
00116 *>          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
00117 *>          On exit, if COMPQ = 'V', Q has been postmultiplied by the
00118 *>          orthogonal transformation matrix which reorders T; the
00119 *>          leading M columns of Q form an orthonormal basis for the
00120 *>          specified invariant subspace.
00121 *>          If COMPQ = 'N', Q is not referenced.
00122 *> \endverbatim
00123 *>
00124 *> \param[in] LDQ
00125 *> \verbatim
00126 *>          LDQ is INTEGER
00127 *>          The leading dimension of the array Q.
00128 *>          LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
00129 *> \endverbatim
00130 *>
00131 *> \param[out] WR
00132 *> \verbatim
00133 *>          WR is REAL array, dimension (N)
00134 *> \endverbatim
00135 *>
00136 *> \param[out] WI
00137 *> \verbatim
00138 *>          WI is REAL array, dimension (N)
00139 *>
00140 *>          The real and imaginary parts, respectively, of the reordered
00141 *>          eigenvalues of T. The eigenvalues are stored in the same
00142 *>          order as on the diagonal of T, with WR(i) = T(i,i) and, if
00143 *>          T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and
00144 *>          WI(i+1) = -WI(i). Note that if a complex eigenvalue is
00145 *>          sufficiently ill-conditioned, then its value may differ
00146 *>          significantly from its value before reordering.
00147 *> \endverbatim
00148 *>
00149 *> \param[out] M
00150 *> \verbatim
00151 *>          M is INTEGER
00152 *>          The dimension of the specified invariant subspace.
00153 *>          0 < = M <= N.
00154 *> \endverbatim
00155 *>
00156 *> \param[out] S
00157 *> \verbatim
00158 *>          S is REAL
00159 *>          If JOB = 'E' or 'B', S is a lower bound on the reciprocal
00160 *>          condition number for the selected cluster of eigenvalues.
00161 *>          S cannot underestimate the true reciprocal condition number
00162 *>          by more than a factor of sqrt(N). If M = 0 or N, S = 1.
00163 *>          If JOB = 'N' or 'V', S is not referenced.
00164 *> \endverbatim
00165 *>
00166 *> \param[out] SEP
00167 *> \verbatim
00168 *>          SEP is REAL
00169 *>          If JOB = 'V' or 'B', SEP is the estimated reciprocal
00170 *>          condition number of the specified invariant subspace. If
00171 *>          M = 0 or N, SEP = norm(T).
00172 *>          If JOB = 'N' or 'E', SEP is not referenced.
00173 *> \endverbatim
00174 *>
00175 *> \param[out] WORK
00176 *> \verbatim
00177 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00178 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00179 *> \endverbatim
00180 *>
00181 *> \param[in] LWORK
00182 *> \verbatim
00183 *>          LWORK is INTEGER
00184 *>          The dimension of the array WORK.
00185 *>          If JOB = 'N', LWORK >= max(1,N);
00186 *>          if JOB = 'E', LWORK >= max(1,M*(N-M));
00187 *>          if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
00188 *>
00189 *>          If LWORK = -1, then a workspace query is assumed; the routine
00190 *>          only calculates the optimal size of the WORK array, returns
00191 *>          this value as the first entry of the WORK array, and no error
00192 *>          message related to LWORK is issued by XERBLA.
00193 *> \endverbatim
00194 *>
00195 *> \param[out] IWORK
00196 *> \verbatim
00197 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
00198 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
00199 *> \endverbatim
00200 *>
00201 *> \param[in] LIWORK
00202 *> \verbatim
00203 *>          LIWORK is INTEGER
00204 *>          The dimension of the array IWORK.
00205 *>          If JOB = 'N' or 'E', LIWORK >= 1;
00206 *>          if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)).
00207 *>
00208 *>          If LIWORK = -1, then a workspace query is assumed; the
00209 *>          routine only calculates the optimal size of the IWORK array,
00210 *>          returns this value as the first entry of the IWORK array, and
00211 *>          no error message related to LIWORK is issued by XERBLA.
00212 *> \endverbatim
00213 *>
00214 *> \param[out] INFO
00215 *> \verbatim
00216 *>          INFO is INTEGER
00217 *>          = 0: successful exit
00218 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00219 *>          = 1: reordering of T failed because some eigenvalues are too
00220 *>               close to separate (the problem is very ill-conditioned);
00221 *>               T may have been partially reordered, and WR and WI
00222 *>               contain the eigenvalues in the same order as in T; S and
00223 *>               SEP (if requested) are set to zero.
00224 *> \endverbatim
00225 *
00226 *  Authors:
00227 *  ========
00228 *
00229 *> \author Univ. of Tennessee 
00230 *> \author Univ. of California Berkeley 
00231 *> \author Univ. of Colorado Denver 
00232 *> \author NAG Ltd. 
00233 *
00234 *> \date April 2012
00235 *
00236 *> \ingroup realOTHERcomputational
00237 *
00238 *> \par Further Details:
00239 *  =====================
00240 *>
00241 *> \verbatim
00242 *>
00243 *>  STRSEN first collects the selected eigenvalues by computing an
00244 *>  orthogonal transformation Z to move them to the top left corner of T.
00245 *>  In other words, the selected eigenvalues are the eigenvalues of T11
00246 *>  in:
00247 *>
00248 *>          Z**T * T * Z = ( T11 T12 ) n1
00249 *>                         (  0  T22 ) n2
00250 *>                            n1  n2
00251 *>
00252 *>  where N = n1+n2 and Z**T means the transpose of Z. The first n1 columns
00253 *>  of Z span the specified invariant subspace of T.
00254 *>
00255 *>  If T has been obtained from the real Schur factorization of a matrix
00256 *>  A = Q*T*Q**T, then the reordered real Schur factorization of A is given
00257 *>  by A = (Q*Z)*(Z**T*T*Z)*(Q*Z)**T, and the first n1 columns of Q*Z span
00258 *>  the corresponding invariant subspace of A.
00259 *>
00260 *>  The reciprocal condition number of the average of the eigenvalues of
00261 *>  T11 may be returned in S. S lies between 0 (very badly conditioned)
00262 *>  and 1 (very well conditioned). It is computed as follows. First we
00263 *>  compute R so that
00264 *>
00265 *>                         P = ( I  R ) n1
00266 *>                             ( 0  0 ) n2
00267 *>                               n1 n2
00268 *>
00269 *>  is the projector on the invariant subspace associated with T11.
00270 *>  R is the solution of the Sylvester equation:
00271 *>
00272 *>                        T11*R - R*T22 = T12.
00273 *>
00274 *>  Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
00275 *>  the two-norm of M. Then S is computed as the lower bound
00276 *>
00277 *>                      (1 + F-norm(R)**2)**(-1/2)
00278 *>
00279 *>  on the reciprocal of 2-norm(P), the true reciprocal condition number.
00280 *>  S cannot underestimate 1 / 2-norm(P) by more than a factor of
00281 *>  sqrt(N).
00282 *>
00283 *>  An approximate error bound for the computed average of the
00284 *>  eigenvalues of T11 is
00285 *>
00286 *>                         EPS * norm(T) / S
00287 *>
00288 *>  where EPS is the machine precision.
00289 *>
00290 *>  The reciprocal condition number of the right invariant subspace
00291 *>  spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
00292 *>  SEP is defined as the separation of T11 and T22:
00293 *>
00294 *>                     sep( T11, T22 ) = sigma-min( C )
00295 *>
00296 *>  where sigma-min(C) is the smallest singular value of the
00297 *>  n1*n2-by-n1*n2 matrix
00298 *>
00299 *>     C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
00300 *>
00301 *>  I(m) is an m by m identity matrix, and kprod denotes the Kronecker
00302 *>  product. We estimate sigma-min(C) by the reciprocal of an estimate of
00303 *>  the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
00304 *>  cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
00305 *>
00306 *>  When SEP is small, small changes in T can cause large changes in
00307 *>  the invariant subspace. An approximate bound on the maximum angular
00308 *>  error in the computed right invariant subspace is
00309 *>
00310 *>                      EPS * norm(T) / SEP
00311 *> \endverbatim
00312 *>
00313 *  =====================================================================
00314       SUBROUTINE STRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
00315      $                   M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
00316 *
00317 *  -- LAPACK computational routine (version 3.4.1) --
00318 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00319 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00320 *     April 2012
00321 *
00322 *     .. Scalar Arguments ..
00323       CHARACTER          COMPQ, JOB
00324       INTEGER            INFO, LDQ, LDT, LIWORK, LWORK, M, N
00325       REAL               S, SEP
00326 *     ..
00327 *     .. Array Arguments ..
00328       LOGICAL            SELECT( * )
00329       INTEGER            IWORK( * )
00330       REAL               Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
00331      $                   WR( * )
00332 *     ..
00333 *
00334 *  =====================================================================
00335 *
00336 *     .. Parameters ..
00337       REAL               ZERO, ONE
00338       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00339 *     ..
00340 *     .. Local Scalars ..
00341       LOGICAL            LQUERY, PAIR, SWAP, WANTBH, WANTQ, WANTS,
00342      $                    WANTSP
00343       INTEGER            IERR, K, KASE, KK, KS, LIWMIN, LWMIN, N1, N2,
00344      $                   NN
00345       REAL               EST, RNORM, SCALE
00346 *     ..
00347 *     .. Local Arrays ..
00348       INTEGER            ISAVE( 3 )
00349 *     ..
00350 *     .. External Functions ..
00351       LOGICAL            LSAME
00352       REAL               SLANGE
00353       EXTERNAL           LSAME, SLANGE
00354 *     ..
00355 *     .. External Subroutines ..
00356       EXTERNAL           SLACN2, SLACPY, STREXC, STRSYL, XERBLA
00357 *     ..
00358 *     .. Intrinsic Functions ..
00359       INTRINSIC          ABS, MAX, SQRT
00360 *     ..
00361 *     .. Executable Statements ..
00362 *
00363 *     Decode and test the input parameters
00364 *
00365       WANTBH = LSAME( JOB, 'B' )
00366       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00367       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
00368       WANTQ = LSAME( COMPQ, 'V' )
00369 *
00370       INFO = 0
00371       LQUERY = ( LWORK.EQ.-1 )
00372       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
00373      $     THEN
00374          INFO = -1
00375       ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
00376          INFO = -2
00377       ELSE IF( N.LT.0 ) THEN
00378          INFO = -4
00379       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00380          INFO = -6
00381       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
00382          INFO = -8
00383       ELSE
00384 *
00385 *        Set M to the dimension of the specified invariant subspace,
00386 *        and test LWORK and LIWORK.
00387 *
00388          M = 0
00389          PAIR = .FALSE.
00390          DO 10 K = 1, N
00391             IF( PAIR ) THEN
00392                PAIR = .FALSE.
00393             ELSE
00394                IF( K.LT.N ) THEN
00395                   IF( T( K+1, K ).EQ.ZERO ) THEN
00396                      IF( SELECT( K ) )
00397      $                  M = M + 1
00398                   ELSE
00399                      PAIR = .TRUE.
00400                      IF( SELECT( K ) .OR. SELECT( K+1 ) )
00401      $                  M = M + 2
00402                   END IF
00403                ELSE
00404                   IF( SELECT( N ) )
00405      $               M = M + 1
00406                END IF
00407             END IF
00408    10    CONTINUE
00409 *
00410          N1 = M
00411          N2 = N - M
00412          NN = N1*N2
00413 *
00414          IF(  WANTSP ) THEN
00415             LWMIN = MAX( 1, 2*NN )
00416             LIWMIN = MAX( 1, NN )
00417          ELSE IF( LSAME( JOB, 'N' ) ) THEN
00418             LWMIN = MAX( 1, N )
00419             LIWMIN = 1
00420          ELSE IF( LSAME( JOB, 'E' ) ) THEN
00421             LWMIN = MAX( 1, NN )
00422             LIWMIN = 1
00423          END IF
00424 *
00425          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00426             INFO = -15
00427          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
00428             INFO = -17
00429          END IF
00430       END IF
00431 *
00432       IF( INFO.EQ.0 ) THEN
00433          WORK( 1 ) = LWMIN
00434          IWORK( 1 ) = LIWMIN
00435       END IF
00436 *
00437       IF( INFO.NE.0 ) THEN
00438          CALL XERBLA( 'STRSEN', -INFO )
00439          RETURN
00440       ELSE IF( LQUERY ) THEN
00441          RETURN
00442       END IF
00443 *
00444 *     Quick return if possible.
00445 *
00446       IF( M.EQ.N .OR. M.EQ.0 ) THEN
00447          IF( WANTS )
00448      $      S = ONE
00449          IF( WANTSP )
00450      $      SEP = SLANGE( '1', N, N, T, LDT, WORK )
00451          GO TO 40
00452       END IF
00453 *
00454 *     Collect the selected blocks at the top-left corner of T.
00455 *
00456       KS = 0
00457       PAIR = .FALSE.
00458       DO 20 K = 1, N
00459          IF( PAIR ) THEN
00460             PAIR = .FALSE.
00461          ELSE
00462             SWAP = SELECT( K )
00463             IF( K.LT.N ) THEN
00464                IF( T( K+1, K ).NE.ZERO ) THEN
00465                   PAIR = .TRUE.
00466                   SWAP = SWAP .OR. SELECT( K+1 )
00467                END IF
00468             END IF
00469             IF( SWAP ) THEN
00470                KS = KS + 1
00471 *
00472 *              Swap the K-th block to position KS.
00473 *
00474                IERR = 0
00475                KK = K
00476                IF( K.NE.KS )
00477      $            CALL STREXC( COMPQ, N, T, LDT, Q, LDQ, KK, KS, WORK,
00478      $                         IERR )
00479                IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
00480 *
00481 *                 Blocks too close to swap: exit.
00482 *
00483                   INFO = 1
00484                   IF( WANTS )
00485      $               S = ZERO
00486                   IF( WANTSP )
00487      $               SEP = ZERO
00488                   GO TO 40
00489                END IF
00490                IF( PAIR )
00491      $            KS = KS + 1
00492             END IF
00493          END IF
00494    20 CONTINUE
00495 *
00496       IF( WANTS ) THEN
00497 *
00498 *        Solve Sylvester equation for R:
00499 *
00500 *           T11*R - R*T22 = scale*T12
00501 *
00502          CALL SLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
00503          CALL STRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
00504      $                LDT, WORK, N1, SCALE, IERR )
00505 *
00506 *        Estimate the reciprocal of the condition number of the cluster
00507 *        of eigenvalues.
00508 *
00509          RNORM = SLANGE( 'F', N1, N2, WORK, N1, WORK )
00510          IF( RNORM.EQ.ZERO ) THEN
00511             S = ONE
00512          ELSE
00513             S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
00514      $          SQRT( RNORM ) )
00515          END IF
00516       END IF
00517 *
00518       IF( WANTSP ) THEN
00519 *
00520 *        Estimate sep(T11,T22).
00521 *
00522          EST = ZERO
00523          KASE = 0
00524    30    CONTINUE
00525          CALL SLACN2( NN, WORK( NN+1 ), WORK, IWORK, EST, KASE, ISAVE )
00526          IF( KASE.NE.0 ) THEN
00527             IF( KASE.EQ.1 ) THEN
00528 *
00529 *              Solve  T11*R - R*T22 = scale*X.
00530 *
00531                CALL STRSYL( 'N', 'N', -1, N1, N2, T, LDT,
00532      $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
00533      $                      IERR )
00534             ELSE
00535 *
00536 *              Solve T11**T*R - R*T22**T = scale*X.
00537 *
00538                CALL STRSYL( 'T', 'T', -1, N1, N2, T, LDT,
00539      $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
00540      $                      IERR )
00541             END IF
00542             GO TO 30
00543          END IF
00544 *
00545          SEP = SCALE / EST
00546       END IF
00547 *
00548    40 CONTINUE
00549 *
00550 *     Store the output eigenvalues in WR and WI.
00551 *
00552       DO 50 K = 1, N
00553          WR( K ) = T( K, K )
00554          WI( K ) = ZERO
00555    50 CONTINUE
00556       DO 60 K = 1, N - 1
00557          IF( T( K+1, K ).NE.ZERO ) THEN
00558             WI( K ) = SQRT( ABS( T( K, K+1 ) ) )*
00559      $                SQRT( ABS( T( K+1, K ) ) )
00560             WI( K+1 ) = -WI( K )
00561          END IF
00562    60 CONTINUE
00563 *
00564       WORK( 1 ) = LWMIN
00565       IWORK( 1 ) = LIWMIN
00566 *
00567       RETURN
00568 *
00569 *     End of STRSEN
00570 *
00571       END
 All Files Functions