SRC/ctrsen.f

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