SRC/strsen.f

   1 *> \brief \b STRSEN
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download STRSEN + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strsen.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strsen.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strsen.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE STRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
  22 *                          M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          COMPQ, JOB
  26 *       INTEGER            INFO, LDQ, LDT, LIWORK, LWORK, M, N
  27 *       REAL               S, SEP
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       LOGICAL            SELECT( * )
  31 *       INTEGER            IWORK( * )
  32 *       REAL               Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
  33 *      $                   WR( * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> STRSEN reorders the real Schur factorization of a real matrix
  43 *> A = Q*T*Q**T, so that a selected cluster of eigenvalues appears in
  44 *> the leading diagonal blocks of the upper quasi-triangular matrix T,
  45 *> and the leading columns of Q form an orthonormal basis of the
  46 *> corresponding right invariant subspace.
  47 *>
  48 *> Optionally the routine computes the reciprocal condition numbers of
  49 *> the cluster of eigenvalues and/or the invariant subspace.
  50 *>
  51 *> T must be in Schur canonical form (as returned by SHSEQR), that is,
  52 *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
  53 *> 2-by-2 diagonal block has its diagonal elements equal and its
  54 *> off-diagonal elements of opposite sign.
  55 *> \endverbatim
  56 *
  57 *  Arguments:
  58 *  ==========
  59 *
  60 *> \param[in] JOB
  61 *> \verbatim
  62 *>          JOB is CHARACTER*1
  63 *>          Specifies whether condition numbers are required for the
  64 *>          cluster of eigenvalues (S) or the invariant subspace (SEP):
  65 *>          = 'N': none;
  66 *>          = 'E': for eigenvalues only (S);
  67 *>          = 'V': for invariant subspace only (SEP);
  68 *>          = 'B': for both eigenvalues and invariant subspace (S and
  69 *>                 SEP).
  70 *> \endverbatim
  71 *>
  72 *> \param[in] COMPQ
  73 *> \verbatim
  74 *>          COMPQ is CHARACTER*1
  75 *>          = 'V': update the matrix Q of Schur vectors;
  76 *>          = 'N': do not update Q.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] SELECT
  80 *> \verbatim
  81 *>          SELECT is LOGICAL array, dimension (N)
  82 *>          SELECT specifies the eigenvalues in the selected cluster. To
  83 *>          select a real eigenvalue w(j), SELECT(j) must be set to
  84 *>          .TRUE.. To select a complex conjugate pair of eigenvalues
  85 *>          w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
  86 *>          either SELECT(j) or SELECT(j+1) or both must be set to
  87 *>          .TRUE.; a complex conjugate pair of eigenvalues must be
  88 *>          either both included in the cluster or both excluded.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] N
  92 *> \verbatim
  93 *>          N is INTEGER
  94 *>          The order of the matrix T. N >= 0.
  95 *> \endverbatim
  96 *>
  97 *> \param[in,out] T
  98 *> \verbatim
  99 *>          T is REAL array, dimension (LDT,N)
 100 *>          On entry, the upper quasi-triangular matrix T, in Schur
 101 *>          canonical form.
 102 *>          On exit, T is overwritten by the reordered matrix T, again in
 103 *>          Schur canonical form, with the selected eigenvalues in the
 104 *>          leading diagonal blocks.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] LDT
 108 *> \verbatim
 109 *>          LDT is INTEGER
 110 *>          The leading dimension of the array T. LDT >= max(1,N).
 111 *> \endverbatim
 112 *>
 113 *> \param[in,out] Q
 114 *> \verbatim
 115 *>          Q is REAL array, dimension (LDQ,N)
 116 *>          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
 117 *>          On exit, if COMPQ = 'V', Q has been postmultiplied by the
 118 *>          orthogonal transformation matrix which reorders T; the
 119 *>          leading M columns of Q form an orthonormal basis for the
 120 *>          specified invariant subspace.
 121 *>          If COMPQ = 'N', Q is not referenced.
 122 *> \endverbatim
 123 *>
 124 *> \param[in] LDQ
 125 *> \verbatim
 126 *>          LDQ is INTEGER
 127 *>          The leading dimension of the array Q.
 128 *>          LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
 129 *> \endverbatim
 130 *>
 131 *> \param[out] WR
 132 *> \verbatim
 133 *>          WR is REAL array, dimension (N)
 134 *> \endverbatim
 135 *>
 136 *> \param[out] WI
 137 *> \verbatim
 138 *>          WI is REAL array, dimension (N)
 139 *>
 140 *>          The real and imaginary parts, respectively, of the reordered
 141 *>          eigenvalues of T. The eigenvalues are stored in the same
 142 *>          order as on the diagonal of T, with WR(i) = T(i,i) and, if
 143 *>          T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and
 144 *>          WI(i+1) = -WI(i). Note that if a complex eigenvalue is
 145 *>          sufficiently ill-conditioned, then its value may differ
 146 *>          significantly from its value before reordering.
 147 *> \endverbatim
 148 *>
 149 *> \param[out] M
 150 *> \verbatim
 151 *>          M is INTEGER
 152 *>          The dimension of the specified invariant subspace.
 153 *>          0 < = M <= N.
 154 *> \endverbatim
 155 *>
 156 *> \param[out] S
 157 *> \verbatim
 158 *>          S is REAL
 159 *>          If JOB = 'E' or 'B', S is a lower bound on the reciprocal
 160 *>          condition number for the selected cluster of eigenvalues.
 161 *>          S cannot underestimate the true reciprocal condition number
 162 *>          by more than a factor of sqrt(N). If M = 0 or N, S = 1.
 163 *>          If JOB = 'N' or 'V', S is not referenced.
 164 *> \endverbatim
 165 *>
 166 *> \param[out] SEP
 167 *> \verbatim
 168 *>          SEP is REAL
 169 *>          If JOB = 'V' or 'B', SEP is the estimated reciprocal
 170 *>          condition number of the specified invariant subspace. If
 171 *>          M = 0 or N, SEP = norm(T).
 172 *>          If JOB = 'N' or 'E', SEP is not referenced.
 173 *> \endverbatim
 174 *>
 175 *> \param[out] WORK
 176 *> \verbatim
 177 *>          WORK is REAL array, dimension (MAX(1,LWORK))
 178 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 179 *> \endverbatim
 180 *>
 181 *> \param[in] LWORK
 182 *> \verbatim
 183 *>          LWORK is INTEGER
 184 *>          The dimension of the array WORK.
 185 *>          If JOB = 'N', LWORK >= max(1,N);
 186 *>          if JOB = 'E', LWORK >= max(1,M*(N-M));
 187 *>          if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
 188 *>
 189 *>          If LWORK = -1, then a workspace query is assumed; the routine
 190 *>          only calculates the optimal size of the WORK array, returns
 191 *>          this value as the first entry of the WORK array, and no error
 192 *>          message related to LWORK is issued by XERBLA.
 193 *> \endverbatim
 194 *>
 195 *> \param[out] IWORK
 196 *> \verbatim
 197 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
 198 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
 199 *> \endverbatim
 200 *>
 201 *> \param[in] LIWORK
 202 *> \verbatim
 203 *>          LIWORK is INTEGER
 204 *>          The dimension of the array IWORK.
 205 *>          If JOB = 'N' or 'E', LIWORK >= 1;
 206 *>          if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)).
 207 *>
 208 *>          If LIWORK = -1, then a workspace query is assumed; the
 209 *>          routine only calculates the optimal size of the IWORK array,
 210 *>          returns this value as the first entry of the IWORK array, and
 211 *>          no error message related to LIWORK is issued by XERBLA.
 212 *> \endverbatim
 213 *>
 214 *> \param[out] INFO
 215 *> \verbatim
 216 *>          INFO is INTEGER
 217 *>          = 0: successful exit
 218 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 219 *>          = 1: reordering of T failed because some eigenvalues are too
 220 *>               close to separate (the problem is very ill-conditioned);
 221 *>               T may have been partially reordered, and WR and WI
 222 *>               contain the eigenvalues in the same order as in T; S and
 223 *>               SEP (if requested) are set to zero.
 224 *> \endverbatim
 225 *
 226 *  Authors:
 227 *  ========
 228 *
 229 *> \author Univ. of Tennessee
 230 *> \author Univ. of California Berkeley
 231 *> \author Univ. of Colorado Denver
 232 *> \author NAG Ltd.
 233 *
 234 *> \date April 2012
 235 *
 236 *> \ingroup realOTHERcomputational
 237 *
 238 *> \par Further Details:
 239 *  =====================
 240 *>
 241 *> \verbatim
 242 *>
 243 *>  STRSEN first collects the selected eigenvalues by computing an
 244 *>  orthogonal transformation Z to move them to the top left corner of T.
 245 *>  In other words, the selected eigenvalues are the eigenvalues of T11
 246 *>  in:
 247 *>
 248 *>          Z**T * T * Z = ( T11 T12 ) n1
 249 *>                         (  0  T22 ) n2
 250 *>                            n1  n2
 251 *>
 252 *>  where N = n1+n2 and Z**T means the transpose of Z. The first n1 columns
 253 *>  of Z span the specified invariant subspace of T.
 254 *>
 255 *>  If T has been obtained from the real Schur factorization of a matrix
 256 *>  A = Q*T*Q**T, then the reordered real Schur factorization of A is given
 257 *>  by A = (Q*Z)*(Z**T*T*Z)*(Q*Z)**T, and the first n1 columns of Q*Z span
 258 *>  the corresponding invariant subspace of A.
 259 *>
 260 *>  The reciprocal condition number of the average of the eigenvalues of
 261 *>  T11 may be returned in S. S lies between 0 (very badly conditioned)
 262 *>  and 1 (very well conditioned). It is computed as follows. First we
 263 *>  compute R so that
 264 *>
 265 *>                         P = ( I  R ) n1
 266 *>                             ( 0  0 ) n2
 267 *>                               n1 n2
 268 *>
 269 *>  is the projector on the invariant subspace associated with T11.
 270 *>  R is the solution of the Sylvester equation:
 271 *>
 272 *>                        T11*R - R*T22 = T12.
 273 *>
 274 *>  Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
 275 *>  the two-norm of M. Then S is computed as the lower bound
 276 *>
 277 *>                      (1 + F-norm(R)**2)**(-1/2)
 278 *>
 279 *>  on the reciprocal of 2-norm(P), the true reciprocal condition number.
 280 *>  S cannot underestimate 1 / 2-norm(P) by more than a factor of
 281 *>  sqrt(N).
 282 *>
 283 *>  An approximate error bound for the computed average of the
 284 *>  eigenvalues of T11 is
 285 *>
 286 *>                         EPS * norm(T) / S
 287 *>
 288 *>  where EPS is the machine precision.
 289 *>
 290 *>  The reciprocal condition number of the right invariant subspace
 291 *>  spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
 292 *>  SEP is defined as the separation of T11 and T22:
 293 *>
 294 *>                     sep( T11, T22 ) = sigma-min( C )
 295 *>
 296 *>  where sigma-min(C) is the smallest singular value of the
 297 *>  n1*n2-by-n1*n2 matrix
 298 *>
 299 *>     C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
 300 *>
 301 *>  I(m) is an m by m identity matrix, and kprod denotes the Kronecker
 302 *>  product. We estimate sigma-min(C) by the reciprocal of an estimate of
 303 *>  the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
 304 *>  cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
 305 *>
 306 *>  When SEP is small, small changes in T can cause large changes in
 307 *>  the invariant subspace. An approximate bound on the maximum angular
 308 *>  error in the computed right invariant subspace is
 309 *>
 310 *>                      EPS * norm(T) / SEP
 311 *> \endverbatim
 312 *>
 313 *  =====================================================================
 314       SUBROUTINE STRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
 315      $                   M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
 316 *
 317 *  -- LAPACK computational routine (version 3.4.1) --
 318 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 319 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 320 *     April 2012
 321 *
 322 *     .. Scalar Arguments ..
 323       CHARACTER          COMPQ, JOB
 324       INTEGER            INFO, LDQ, LDT, LIWORK, LWORK, M, N
 325       REAL               S, SEP
 326 *     ..
 327 *     .. Array Arguments ..
 328       LOGICAL            SELECT( * )
 329       INTEGER            IWORK( * )
 330       REAL               Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
 331      $                   WR( * )
 332 *     ..
 333 *
 334 *  =====================================================================
 335 *
 336 *     .. Parameters ..
 337       REAL               ZERO, ONE
 338       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 339 *     ..
 340 *     .. Local Scalars ..
 341       LOGICAL            LQUERY, PAIR, SWAP, WANTBH, WANTQ, WANTS,
 342      $                    WANTSP
 343       INTEGER            IERR, K, KASE, KK, KS, LIWMIN, LWMIN, N1, N2,
 344      $                   NN
 345       REAL               EST, RNORM, SCALE
 346 *     ..
 347 *     .. Local Arrays ..
 348       INTEGER            ISAVE( 3 )
 349 *     ..
 350 *     .. External Functions ..
 351       LOGICAL            LSAME
 352       REAL               SLANGE
 353       EXTERNAL           LSAME, SLANGE
 354 *     ..
 355 *     .. External Subroutines ..
 356       EXTERNAL           SLACN2, SLACPY, STREXC, STRSYL, XERBLA
 357 *     ..
 358 *     .. Intrinsic Functions ..
 359       INTRINSIC          ABS, MAX, SQRT
 360 *     ..
 361 *     .. Executable Statements ..
 362 *
 363 *     Decode and test the input parameters
 364 *
 365       WANTBH = LSAME( JOB, 'B' )
 366       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
 367       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
 368       WANTQ = LSAME( COMPQ, 'V' )
 369 *
 370       INFO = 0
 371       LQUERY = ( LWORK.EQ.-1 )
 372       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
 373      $     THEN
 374          INFO = -1
 375       ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
 376          INFO = -2
 377       ELSE IF( N.LT.0 ) THEN
 378          INFO = -4
 379       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
 380          INFO = -6
 381       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
 382          INFO = -8
 383       ELSE
 384 *
 385 *        Set M to the dimension of the specified invariant subspace,
 386 *        and test LWORK and LIWORK.
 387 *
 388          M = 0
 389          PAIR = .FALSE.
 390          DO 10 K = 1, N
 391             IF( PAIR ) THEN
 392                PAIR = .FALSE.
 393             ELSE
 394                IF( K.LT.N ) THEN
 395                   IF( T( K+1, K ).EQ.ZERO ) THEN
 396                      IF( SELECT( K ) )
 397      $                  M = M + 1
 398                   ELSE
 399                      PAIR = .TRUE.
 400                      IF( SELECT( K ) .OR. SELECT( K+1 ) )
 401      $                  M = M + 2
 402                   END IF
 403                ELSE
 404                   IF( SELECT( N ) )
 405      $               M = M + 1
 406                END IF
 407             END IF
 408    10    CONTINUE
 409 *
 410          N1 = M
 411          N2 = N - M
 412          NN = N1*N2
 413 *
 414          IF(  WANTSP ) THEN
 415             LWMIN = MAX( 1, 2*NN )
 416             LIWMIN = MAX( 1, NN )
 417          ELSE IF( LSAME( JOB, 'N' ) ) THEN
 418             LWMIN = MAX( 1, N )
 419             LIWMIN = 1
 420          ELSE IF( LSAME( JOB, 'E' ) ) THEN
 421             LWMIN = MAX( 1, NN )
 422             LIWMIN = 1
 423          END IF
 424 *
 425          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 426             INFO = -15
 427          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
 428             INFO = -17
 429          END IF
 430       END IF
 431 *
 432       IF( INFO.EQ.0 ) THEN
 433          WORK( 1 ) = LWMIN
 434          IWORK( 1 ) = LIWMIN
 435       END IF
 436 *
 437       IF( INFO.NE.0 ) THEN
 438          CALL XERBLA( 'STRSEN', -INFO )
 439          RETURN
 440       ELSE IF( LQUERY ) THEN
 441          RETURN
 442       END IF
 443 *
 444 *     Quick return if possible.
 445 *
 446       IF( M.EQ.N .OR. M.EQ.0 ) THEN
 447          IF( WANTS )
 448      $      S = ONE
 449          IF( WANTSP )
 450      $      SEP = SLANGE( '1', N, N, T, LDT, WORK )
 451          GO TO 40
 452       END IF
 453 *
 454 *     Collect the selected blocks at the top-left corner of T.
 455 *
 456       KS = 0
 457       PAIR = .FALSE.
 458       DO 20 K = 1, N
 459          IF( PAIR ) THEN
 460             PAIR = .FALSE.
 461          ELSE
 462             SWAP = SELECT( K )
 463             IF( K.LT.N ) THEN
 464                IF( T( K+1, K ).NE.ZERO ) THEN
 465                   PAIR = .TRUE.
 466                   SWAP = SWAP .OR. SELECT( K+1 )
 467                END IF
 468             END IF
 469             IF( SWAP ) THEN
 470                KS = KS + 1
 471 *
 472 *              Swap the K-th block to position KS.
 473 *
 474                IERR = 0
 475                KK = K
 476                IF( K.NE.KS )
 477      $            CALL STREXC( COMPQ, N, T, LDT, Q, LDQ, KK, KS, WORK,
 478      $                         IERR )
 479                IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
 480 *
 481 *                 Blocks too close to swap: exit.
 482 *
 483                   INFO = 1
 484                   IF( WANTS )
 485      $               S = ZERO
 486                   IF( WANTSP )
 487      $               SEP = ZERO
 488                   GO TO 40
 489                END IF
 490                IF( PAIR )
 491      $            KS = KS + 1
 492             END IF
 493          END IF
 494    20 CONTINUE
 495 *
 496       IF( WANTS ) THEN
 497 *
 498 *        Solve Sylvester equation for R:
 499 *
 500 *           T11*R - R*T22 = scale*T12
 501 *
 502          CALL SLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
 503          CALL STRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
 504      $                LDT, WORK, N1, SCALE, IERR )
 505 *
 506 *        Estimate the reciprocal of the condition number of the cluster
 507 *        of eigenvalues.
 508 *
 509          RNORM = SLANGE( 'F', N1, N2, WORK, N1, WORK )
 510          IF( RNORM.EQ.ZERO ) THEN
 511             S = ONE
 512          ELSE
 513             S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
 514      $          SQRT( RNORM ) )
 515          END IF
 516       END IF
 517 *
 518       IF( WANTSP ) THEN
 519 *
 520 *        Estimate sep(T11,T22).
 521 *
 522          EST = ZERO
 523          KASE = 0
 524    30    CONTINUE
 525          CALL SLACN2( NN, WORK( NN+1 ), WORK, IWORK, EST, KASE, ISAVE )
 526          IF( KASE.NE.0 ) THEN
 527             IF( KASE.EQ.1 ) THEN
 528 *
 529 *              Solve  T11*R - R*T22 = scale*X.
 530 *
 531                CALL STRSYL( 'N', 'N', -1, N1, N2, T, LDT,
 532      $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
 533      $                      IERR )
 534             ELSE
 535 *
 536 *              Solve T11**T*R - R*T22**T = scale*X.
 537 *
 538                CALL STRSYL( 'T', 'T', -1, N1, N2, T, LDT,
 539      $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
 540      $                      IERR )
 541             END IF
 542             GO TO 30
 543          END IF
 544 *
 545          SEP = SCALE / EST
 546       END IF
 547 *
 548    40 CONTINUE
 549 *
 550 *     Store the output eigenvalues in WR and WI.
 551 *
 552       DO 50 K = 1, N
 553          WR( K ) = T( K, K )
 554          WI( K ) = ZERO
 555    50 CONTINUE
 556       DO 60 K = 1, N - 1
 557          IF( T( K+1, K ).NE.ZERO ) THEN
 558             WI( K ) = SQRT( ABS( T( K, K+1 ) ) )*
 559      $                SQRT( ABS( T( K+1, K ) ) )
 560             WI( K+1 ) = -WI( K )
 561          END IF
 562    60 CONTINUE
 563 *
 564       WORK( 1 ) = LWMIN
 565       IWORK( 1 ) = LIWMIN
 566 *
 567       RETURN
 568 *
 569 *     End of STRSEN
 570 *
 571       END