Return to ztrsen.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief \b ZTRSEN
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZTRSEN + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsen.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsen.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsen.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZTRSEN( 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: *       DOUBLE PRECISION   S, SEP
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       LOGICAL            SELECT( * )
   31: *       COMPLEX*16         Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
   32: *       ..
   33: *  
   34: *
   35: *> \par Purpose:
   36: *  =============
   37: *>
   38: *> \verbatim
   39: *>
   40: *> ZTRSEN 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*16 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*16 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*16 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 DOUBLE PRECISION
  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 DOUBLE PRECISION
  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*16 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 complex16OTHERcomputational
  188: *
  189: *> \par Further Details:
  190: *  =====================
  191: *>
  192: *> \verbatim
  193: *>
  194: *>  ZTRSEN 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 ZTRSEN( 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:       DOUBLE PRECISION   S, SEP
  276: *     ..
  277: *     .. Array Arguments ..
  278:       LOGICAL            SELECT( * )
  279:       COMPLEX*16         Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
  280: *     ..
  281: *
  282: *  =====================================================================
  283: *
  284: *     .. Parameters ..
  285:       DOUBLE PRECISION   ZERO, ONE
  286:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  287: *     ..
  288: *     .. Local Scalars ..
  289:       LOGICAL            LQUERY, WANTBH, WANTQ, WANTS, WANTSP
  290:       INTEGER            IERR, K, KASE, KS, LWMIN, N1, N2, NN
  291:       DOUBLE PRECISION   EST, RNORM, SCALE
  292: *     ..
  293: *     .. Local Arrays ..
  294:       INTEGER            ISAVE( 3 )
  295:       DOUBLE PRECISION   RWORK( 1 )
  296: *     ..
  297: *     .. External Functions ..
  298:       LOGICAL            LSAME
  299:       DOUBLE PRECISION   ZLANGE
  300:       EXTERNAL           LSAME, ZLANGE
  301: *     ..
  302: *     .. External Subroutines ..
  303:       EXTERNAL           XERBLA, ZLACN2, ZLACPY, ZTREXC, ZTRSYL
  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( 'ZTRSEN', -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 = ZLANGE( '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 ZTREXC( 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 ZLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
  397:          CALL ZTRSYL( '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 = ZLANGE( '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 ZLACN2( 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 ZTRSYL( '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 ZTRSYL( '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 ZTRSEN
  455: *
  456:       END