Return to ztgsna.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence Lapack.

    1: *> \brief \b ZTGSNA
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download ZTGSNA + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsna.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsna.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsna.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
   22: *                          LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
   23: *                          IWORK, INFO )
   24: * 
   25: *       .. Scalar Arguments ..
   26: *       CHARACTER          HOWMNY, JOB
   27: *       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       LOGICAL            SELECT( * )
   31: *       INTEGER            IWORK( * )
   32: *       DOUBLE PRECISION   DIF( * ), S( * )
   33: *       COMPLEX*16         A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
   34: *      $                   VR( LDVR, * ), WORK( * )
   35: *       ..
   36: *  
   37: *
   38: *> \par Purpose:
   39: *  =============
   40: *>
   41: *> \verbatim
   42: *>
   43: *> ZTGSNA estimates reciprocal condition numbers for specified
   44: *> eigenvalues and/or eigenvectors of a matrix pair (A, B).
   45: *>
   46: *> (A, B) must be in generalized Schur canonical form, that is, A and
   47: *> B are both upper triangular.
   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
   57: *>          eigenvalues (S) or eigenvectors (DIF):
   58: *>          = 'E': for eigenvalues only (S);
   59: *>          = 'V': for eigenvectors only (DIF);
   60: *>          = 'B': for both eigenvalues and eigenvectors (S and DIF).
   61: *> \endverbatim
   62: *>
   63: *> \param[in] HOWMNY
   64: *> \verbatim
   65: *>          HOWMNY is CHARACTER*1
   66: *>          = 'A': compute condition numbers for all eigenpairs;
   67: *>          = 'S': compute condition numbers for selected eigenpairs
   68: *>                 specified by the array SELECT.
   69: *> \endverbatim
   70: *>
   71: *> \param[in] SELECT
   72: *> \verbatim
   73: *>          SELECT is LOGICAL array, dimension (N)
   74: *>          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
   75: *>          condition numbers are required. To select condition numbers
   76: *>          for the corresponding j-th eigenvalue and/or eigenvector,
   77: *>          SELECT(j) must be set to .TRUE..
   78: *>          If HOWMNY = 'A', SELECT is not referenced.
   79: *> \endverbatim
   80: *>
   81: *> \param[in] N
   82: *> \verbatim
   83: *>          N is INTEGER
   84: *>          The order of the square matrix pair (A, B). N >= 0.
   85: *> \endverbatim
   86: *>
   87: *> \param[in] A
   88: *> \verbatim
   89: *>          A is COMPLEX*16 array, dimension (LDA,N)
   90: *>          The upper triangular matrix A in the pair (A,B).
   91: *> \endverbatim
   92: *>
   93: *> \param[in] LDA
   94: *> \verbatim
   95: *>          LDA is INTEGER
   96: *>          The leading dimension of the array A. LDA >= max(1,N).
   97: *> \endverbatim
   98: *>
   99: *> \param[in] B
  100: *> \verbatim
  101: *>          B is COMPLEX*16 array, dimension (LDB,N)
  102: *>          The upper triangular matrix B in the pair (A, B).
  103: *> \endverbatim
  104: *>
  105: *> \param[in] LDB
  106: *> \verbatim
  107: *>          LDB is INTEGER
  108: *>          The leading dimension of the array B. LDB >= max(1,N).
  109: *> \endverbatim
  110: *>
  111: *> \param[in] VL
  112: *> \verbatim
  113: *>          VL is COMPLEX*16 array, dimension (LDVL,M)
  114: *>          IF JOB = 'E' or 'B', VL must contain left eigenvectors of
  115: *>          (A, B), corresponding to the eigenpairs specified by HOWMNY
  116: *>          and SELECT.  The eigenvectors must be stored in consecutive
  117: *>          columns of VL, as returned by ZTGEVC.
  118: *>          If JOB = 'V', VL is not referenced.
  119: *> \endverbatim
  120: *>
  121: *> \param[in] LDVL
  122: *> \verbatim
  123: *>          LDVL is INTEGER
  124: *>          The leading dimension of the array VL. LDVL >= 1; and
  125: *>          If JOB = 'E' or 'B', LDVL >= N.
  126: *> \endverbatim
  127: *>
  128: *> \param[in] VR
  129: *> \verbatim
  130: *>          VR is COMPLEX*16 array, dimension (LDVR,M)
  131: *>          IF JOB = 'E' or 'B', VR must contain right eigenvectors of
  132: *>          (A, B), corresponding to the eigenpairs specified by HOWMNY
  133: *>          and SELECT.  The eigenvectors must be stored in consecutive
  134: *>          columns of VR, as returned by ZTGEVC.
  135: *>          If JOB = 'V', VR is not referenced.
  136: *> \endverbatim
  137: *>
  138: *> \param[in] LDVR
  139: *> \verbatim
  140: *>          LDVR is INTEGER
  141: *>          The leading dimension of the array VR. LDVR >= 1;
  142: *>          If JOB = 'E' or 'B', LDVR >= N.
  143: *> \endverbatim
  144: *>
  145: *> \param[out] S
  146: *> \verbatim
  147: *>          S is DOUBLE PRECISION array, dimension (MM)
  148: *>          If JOB = 'E' or 'B', the reciprocal condition numbers of the
  149: *>          selected eigenvalues, stored in consecutive elements of the
  150: *>          array.
  151: *>          If JOB = 'V', S is not referenced.
  152: *> \endverbatim
  153: *>
  154: *> \param[out] DIF
  155: *> \verbatim
  156: *>          DIF is DOUBLE PRECISION array, dimension (MM)
  157: *>          If JOB = 'V' or 'B', the estimated reciprocal condition
  158: *>          numbers of the selected eigenvectors, stored in consecutive
  159: *>          elements of the array.
  160: *>          If the eigenvalues cannot be reordered to compute DIF(j),
  161: *>          DIF(j) is set to 0; this can only occur when the true value
  162: *>          would be very small anyway.
  163: *>          For each eigenvalue/vector specified by SELECT, DIF stores
  164: *>          a Frobenius norm-based estimate of Difl.
  165: *>          If JOB = 'E', DIF is not referenced.
  166: *> \endverbatim
  167: *>
  168: *> \param[in] MM
  169: *> \verbatim
  170: *>          MM is INTEGER
  171: *>          The number of elements in the arrays S and DIF. MM >= M.
  172: *> \endverbatim
  173: *>
  174: *> \param[out] M
  175: *> \verbatim
  176: *>          M is INTEGER
  177: *>          The number of elements of the arrays S and DIF used to store
  178: *>          the specified condition numbers; for each selected eigenvalue
  179: *>          one element is used. If HOWMNY = 'A', M is set to N.
  180: *> \endverbatim
  181: *>
  182: *> \param[out] WORK
  183: *> \verbatim
  184: *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
  185: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  186: *> \endverbatim
  187: *>
  188: *> \param[in] LWORK
  189: *> \verbatim
  190: *>          LWORK is INTEGER
  191: *>          The dimension of the array WORK. LWORK >= max(1,N).
  192: *>          If JOB = 'V' or 'B', LWORK >= max(1,2*N*N).
  193: *> \endverbatim
  194: *>
  195: *> \param[out] IWORK
  196: *> \verbatim
  197: *>          IWORK is INTEGER array, dimension (N+2)
  198: *>          If JOB = 'E', IWORK is not referenced.
  199: *> \endverbatim
  200: *>
  201: *> \param[out] INFO
  202: *> \verbatim
  203: *>          INFO is INTEGER
  204: *>          = 0: Successful exit
  205: *>          < 0: If INFO = -i, the i-th argument had an illegal value
  206: *> \endverbatim
  207: *
  208: *  Authors:
  209: *  ========
  210: *
  211: *> \author Univ. of Tennessee 
  212: *> \author Univ. of California Berkeley 
  213: *> \author Univ. of Colorado Denver 
  214: *> \author NAG Ltd. 
  215: *
  216: *> \date November 2011
  217: *
  218: *> \ingroup complex16OTHERcomputational
  219: *
  220: *> \par Further Details:
  221: *  =====================
  222: *>
  223: *> \verbatim
  224: *>
  225: *>  The reciprocal of the condition number of the i-th generalized
  226: *>  eigenvalue w = (a, b) is defined as
  227: *>
  228: *>          S(I) = (|v**HAu|**2 + |v**HBu|**2)**(1/2) / (norm(u)*norm(v))
  229: *>
  230: *>  where u and v are the right and left eigenvectors of (A, B)
  231: *>  corresponding to w; |z| denotes the absolute value of the complex
  232: *>  number, and norm(u) denotes the 2-norm of the vector u. The pair
  233: *>  (a, b) corresponds to an eigenvalue w = a/b (= v**HAu/v**HBu) of the
  234: *>  matrix pair (A, B). If both a and b equal zero, then (A,B) is
  235: *>  singular and S(I) = -1 is returned.
  236: *>
  237: *>  An approximate error bound on the chordal distance between the i-th
  238: *>  computed generalized eigenvalue w and the corresponding exact
  239: *>  eigenvalue lambda is
  240: *>
  241: *>          chord(w, lambda) <=   EPS * norm(A, B) / S(I),
  242: *>
  243: *>  where EPS is the machine precision.
  244: *>
  245: *>  The reciprocal of the condition number of the right eigenvector u
  246: *>  and left eigenvector v corresponding to the generalized eigenvalue w
  247: *>  is defined as follows. Suppose
  248: *>
  249: *>                   (A, B) = ( a   *  ) ( b  *  )  1
  250: *>                            ( 0  A22 ),( 0 B22 )  n-1
  251: *>                              1  n-1     1 n-1
  252: *>
  253: *>  Then the reciprocal condition number DIF(I) is
  254: *>
  255: *>          Difl[(a, b), (A22, B22)]  = sigma-min( Zl )
  256: *>
  257: *>  where sigma-min(Zl) denotes the smallest singular value of
  258: *>
  259: *>         Zl = [ kron(a, In-1) -kron(1, A22) ]
  260: *>              [ kron(b, In-1) -kron(1, B22) ].
  261: *>
  262: *>  Here In-1 is the identity matrix of size n-1 and X**H is the conjugate
  263: *>  transpose of X. kron(X, Y) is the Kronecker product between the
  264: *>  matrices X and Y.
  265: *>
  266: *>  We approximate the smallest singular value of Zl with an upper
  267: *>  bound. This is done by ZLATDF.
  268: *>
  269: *>  An approximate error bound for a computed eigenvector VL(i) or
  270: *>  VR(i) is given by
  271: *>
  272: *>                      EPS * norm(A, B) / DIF(i).
  273: *>
  274: *>  See ref. [2-3] for more details and further references.
  275: *> \endverbatim
  276: *
  277: *> \par Contributors:
  278: *  ==================
  279: *>
  280: *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
  281: *>     Umea University, S-901 87 Umea, Sweden.
  282: *
  283: *> \par References:
  284: *  ================
  285: *>
  286: *> \verbatim
  287: *>
  288: *>  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
  289: *>      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
  290: *>      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
  291: *>      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
  292: *>
  293: *>  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
  294: *>      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
  295: *>      Estimation: Theory, Algorithms and Software, Report
  296: *>      UMINF - 94.04, Department of Computing Science, Umea University,
  297: *>      S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
  298: *>      To appear in Numerical Algorithms, 1996.
  299: *>
  300: *>  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
  301: *>      for Solving the Generalized Sylvester Equation and Estimating the
  302: *>      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
  303: *>      Department of Computing Science, Umea University, S-901 87 Umea,
  304: *>      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
  305: *>      Note 75.
  306: *>      To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
  307: *> \endverbatim
  308: *>
  309: *  =====================================================================
  310:       SUBROUTINE ZTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
  311:      $                   LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
  312:      $                   IWORK, INFO )
  313: *
  314: *  -- LAPACK computational routine (version 3.4.0) --
  315: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  316: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  317: *     November 2011
  318: *
  319: *     .. Scalar Arguments ..
  320:       CHARACTER          HOWMNY, JOB
  321:       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
  322: *     ..
  323: *     .. Array Arguments ..
  324:       LOGICAL            SELECT( * )
  325:       INTEGER            IWORK( * )
  326:       DOUBLE PRECISION   DIF( * ), S( * )
  327:       COMPLEX*16         A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
  328:      $                   VR( LDVR, * ), WORK( * )
  329: *     ..
  330: *
  331: *  =====================================================================
  332: *
  333: *     .. Parameters ..
  334:       DOUBLE PRECISION   ZERO, ONE
  335:       INTEGER            IDIFJB
  336:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, IDIFJB = 3 )
  337: *     ..
  338: *     .. Local Scalars ..
  339:       LOGICAL            LQUERY, SOMCON, WANTBH, WANTDF, WANTS
  340:       INTEGER            I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2
  341:       DOUBLE PRECISION   BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM
  342:       COMPLEX*16         YHAX, YHBX
  343: *     ..
  344: *     .. Local Arrays ..
  345:       COMPLEX*16         DUMMY( 1 ), DUMMY1( 1 )
  346: *     ..
  347: *     .. External Functions ..
  348:       LOGICAL            LSAME
  349:       DOUBLE PRECISION   DLAMCH, DLAPY2, DZNRM2
  350:       COMPLEX*16         ZDOTC
  351:       EXTERNAL           LSAME, DLAMCH, DLAPY2, DZNRM2, ZDOTC
  352: *     ..
  353: *     .. External Subroutines ..
  354:       EXTERNAL           DLABAD, XERBLA, ZGEMV, ZLACPY, ZTGEXC, ZTGSYL
  355: *     ..
  356: *     .. Intrinsic Functions ..
  357:       INTRINSIC          ABS, DCMPLX, MAX
  358: *     ..
  359: *     .. Executable Statements ..
  360: *
  361: *     Decode and test the input parameters
  362: *
  363:       WANTBH = LSAME( JOB, 'B' )
  364:       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
  365:       WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
  366: *
  367:       SOMCON = LSAME( HOWMNY, 'S' )
  368: *
  369:       INFO = 0
  370:       LQUERY = ( LWORK.EQ.-1 )
  371: *
  372:       IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
  373:          INFO = -1
  374:       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
  375:          INFO = -2
  376:       ELSE IF( N.LT.0 ) THEN
  377:          INFO = -4
  378:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  379:          INFO = -6
  380:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  381:          INFO = -8
  382:       ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
  383:          INFO = -10
  384:       ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
  385:          INFO = -12
  386:       ELSE
  387: *
  388: *        Set M to the number of eigenpairs for which condition numbers
  389: *        are required, and test MM.
  390: *
  391:          IF( SOMCON ) THEN
  392:             M = 0
  393:             DO 10 K = 1, N
  394:                IF( SELECT( K ) )
  395:      $            M = M + 1
  396:    10       CONTINUE
  397:          ELSE
  398:             M = N
  399:          END IF
  400: *
  401:          IF( N.EQ.0 ) THEN
  402:             LWMIN = 1
  403:          ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
  404:             LWMIN = 2*N*N
  405:          ELSE
  406:             LWMIN = N
  407:          END IF
  408:          WORK( 1 ) = LWMIN
  409: *
  410:          IF( MM.LT.M ) THEN
  411:             INFO = -15
  412:          ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
  413:             INFO = -18
  414:          END IF
  415:       END IF
  416: *
  417:       IF( INFO.NE.0 ) THEN
  418:          CALL XERBLA( 'ZTGSNA', -INFO )
  419:          RETURN
  420:       ELSE IF( LQUERY ) THEN
  421:          RETURN
  422:       END IF
  423: *
  424: *     Quick return if possible
  425: *
  426:       IF( N.EQ.0 )
  427:      $   RETURN
  428: *
  429: *     Get machine constants
  430: *
  431:       EPS = DLAMCH( 'P' )
  432:       SMLNUM = DLAMCH( 'S' ) / EPS
  433:       BIGNUM = ONE / SMLNUM
  434:       CALL DLABAD( SMLNUM, BIGNUM )
  435:       KS = 0
  436:       DO 20 K = 1, N
  437: *
  438: *        Determine whether condition numbers are required for the k-th
  439: *        eigenpair.
  440: *
  441:          IF( SOMCON ) THEN
  442:             IF( .NOT.SELECT( K ) )
  443:      $         GO TO 20
  444:          END IF
  445: *
  446:          KS = KS + 1
  447: *
  448:          IF( WANTS ) THEN
  449: *
  450: *           Compute the reciprocal condition number of the k-th
  451: *           eigenvalue.
  452: *
  453:             RNRM = DZNRM2( N, VR( 1, KS ), 1 )
  454:             LNRM = DZNRM2( N, VL( 1, KS ), 1 )
  455:             CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), A, LDA,
  456:      $                  VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
  457:             YHAX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
  458:             CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), B, LDB,
  459:      $                  VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
  460:             YHBX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
  461:             COND = DLAPY2( ABS( YHAX ), ABS( YHBX ) )
  462:             IF( COND.EQ.ZERO ) THEN
  463:                S( KS ) = -ONE
  464:             ELSE
  465:                S( KS ) = COND / ( RNRM*LNRM )
  466:             END IF
  467:          END IF
  468: *
  469:          IF( WANTDF ) THEN
  470:             IF( N.EQ.1 ) THEN
  471:                DIF( KS ) = DLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) )
  472:             ELSE
  473: *
  474: *              Estimate the reciprocal condition number of the k-th
  475: *              eigenvectors.
  476: *
  477: *              Copy the matrix (A, B) to the array WORK and move the
  478: *              (k,k)th pair to the (1,1) position.
  479: *
  480:                CALL ZLACPY( 'Full', N, N, A, LDA, WORK, N )
  481:                CALL ZLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
  482:                IFST = K
  483:                ILST = 1
  484: *
  485:                CALL ZTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ),
  486:      $                      N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR )
  487: *
  488:                IF( IERR.GT.0 ) THEN
  489: *
  490: *                 Ill-conditioned problem - swap rejected.
  491: *
  492:                   DIF( KS ) = ZERO
  493:                ELSE
  494: *
  495: *                 Reordering successful, solve generalized Sylvester
  496: *                 equation for R and L,
  497: *                            A22 * R - L * A11 = A12
  498: *                            B22 * R - L * B11 = B12,
  499: *                 and compute estimate of Difl[(A11,B11), (A22, B22)].
  500: *
  501:                   N1 = 1
  502:                   N2 = N - N1
  503:                   I = N*N + 1
  504:                   CALL ZTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ),
  505:      $                         N, WORK, N, WORK( N1+1 ), N,
  506:      $                         WORK( N*N1+N1+I ), N, WORK( I ), N,
  507:      $                         WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY,
  508:      $                         1, IWORK, IERR )
  509:                END IF
  510:             END IF
  511:          END IF
  512: *
  513:    20 CONTINUE
  514:       WORK( 1 ) = LWMIN
  515:       RETURN
  516: *
  517: *     End of ZTGSNA
  518: *
  519:       END