dtrsen.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 484 行 · 第 1/3 页

</span><span class="comment">*</span><span class="comment">          = 1: reordering of T failed because some eigenvalues are too
</span><span class="comment">*</span><span class="comment">               close to separate (the problem is very ill-conditioned);
</span><span class="comment">*</span><span class="comment">               T may have been partially reordered, and WR and WI
</span><span class="comment">*</span><span class="comment">               contain the eigenvalues in the same order as in T; S and
</span><span class="comment">*</span><span class="comment">               SEP (if requested) are set to zero.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="DTRSEN.153"></a><a href="dtrsen.f.html#DTRSEN.1">DTRSEN</a> first collects the selected eigenvalues by computing an
</span><span class="comment">*</span><span class="comment">  orthogonal transformation Z to move them to the top left corner of T.
</span><span class="comment">*</span><span class="comment">  In other words, the selected eigenvalues are the eigenvalues of T11
</span><span class="comment">*</span><span class="comment">  in:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                Z'*T*Z = ( T11 T12 ) n1
</span><span class="comment">*</span><span class="comment">                         (  0  T22 ) n2
</span><span class="comment">*</span><span class="comment">                            n1  n2
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where N = n1+n2 and Z' means the transpose of Z. The first n1 columns
</span><span class="comment">*</span><span class="comment">  of Z span the specified invariant subspace of T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  If T has been obtained from the real Schur factorization of a matrix
</span><span class="comment">*</span><span class="comment">  A = Q*T*Q', then the reordered real Schur factorization of A is given
</span><span class="comment">*</span><span class="comment">  by A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span
</span><span class="comment">*</span><span class="comment">  the corresponding invariant subspace of A.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The reciprocal condition number of the average of the eigenvalues of
</span><span class="comment">*</span><span class="comment">  T11 may be returned in S. S lies between 0 (very badly conditioned)
</span><span class="comment">*</span><span class="comment">  and 1 (very well conditioned). It is computed as follows. First we
</span><span class="comment">*</span><span class="comment">  compute R so that
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                         P = ( I  R ) n1
</span><span class="comment">*</span><span class="comment">                             ( 0  0 ) n2
</span><span class="comment">*</span><span class="comment">                               n1 n2
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  is the projector on the invariant subspace associated with T11.
</span><span class="comment">*</span><span class="comment">  R is the solution of the Sylvester equation:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                        T11*R - R*T22 = T12.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
</span><span class="comment">*</span><span class="comment">  the two-norm of M. Then S is computed as the lower bound
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                      (1 + F-norm(R)**2)**(-1/2)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  on the reciprocal of 2-norm(P), the true reciprocal condition number.
</span><span class="comment">*</span><span class="comment">  S cannot underestimate 1 / 2-norm(P) by more than a factor of
</span><span class="comment">*</span><span class="comment">  sqrt(N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  An approximate error bound for the computed average of the
</span><span class="comment">*</span><span class="comment">  eigenvalues of T11 is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                         EPS * norm(T) / S
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where EPS is the machine precision.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The reciprocal condition number of the right invariant subspace
</span><span class="comment">*</span><span class="comment">  spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
</span><span class="comment">*</span><span class="comment">  SEP is defined as the separation of T11 and T22:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                     sep( T11, T22 ) = sigma-min( C )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where sigma-min(C) is the smallest singular value of the
</span><span class="comment">*</span><span class="comment">  n1*n2-by-n1*n2 matrix
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  I(m) is an m by m identity matrix, and kprod denotes the Kronecker
</span><span class="comment">*</span><span class="comment">  product. We estimate sigma-min(C) by the reciprocal of an estimate of
</span><span class="comment">*</span><span class="comment">  the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
</span><span class="comment">*</span><span class="comment">  cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  When SEP is small, small changes in T can cause large changes in
</span><span class="comment">*</span><span class="comment">  the invariant subspace. An approximate bound on the maximum angular
</span><span class="comment">*</span><span class="comment">  error in the computed right invariant subspace is
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                      EPS * norm(T) / SEP
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY, PAIR, SWAP, WANTBH, WANTQ, WANTS,
     $                   WANTSP
      INTEGER            IERR, K, KASE, KK, KS, LIWMIN, LWMIN, N1, N2,
     $                   NN
      DOUBLE PRECISION   EST, RNORM, SCALE
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            ISAVE( 3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.239"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      DOUBLE PRECISION   <a name="DLANGE.240"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
      EXTERNAL           <a name="LSAME.241"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="DLANGE.241"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DLACN2.244"></a><a href="dlacn2.f.html#DLACN2.1">DLACN2</a>, <a name="DLACPY.244"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>, <a name="DTREXC.244"></a><a href="dtrexc.f.html#DTREXC.1">DTREXC</a>, <a name="DTRSYL.244"></a><a href="dtrsyl.f.html#DTRSYL.1">DTRSYL</a>, <a name="XERBLA.244"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Decode and test the input parameters
</span><span class="comment">*</span><span class="comment">
</span>      WANTBH = <a name="LSAME.253"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> )
      WANTS = <a name="LSAME.254"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'E'</span> ) .OR. WANTBH
      WANTSP = <a name="LSAME.255"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'V'</span> ) .OR. WANTBH
      WANTQ = <a name="LSAME.256"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'V'</span> )
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      LQUERY = ( LWORK.EQ.-1 )
      IF( .NOT.<a name="LSAME.260"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'N'</span> ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
     $     THEN
         INFO = -1
      ELSE IF( .NOT.<a name="LSAME.263"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( COMPQ, <span class="string">'N'</span> ) .AND. .NOT.WANTQ ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
         INFO = -8
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Set M to the dimension of the specified invariant subspace,
</span><span class="comment">*</span><span class="comment">        and test LWORK and LIWORK.
</span><span class="comment">*</span><span class="comment">
</span>         M = 0
         PAIR = .FALSE.
         DO 10 K = 1, N
            IF( PAIR ) THEN
               PAIR = .FALSE.
            ELSE
               IF( K.LT.N ) THEN
                  IF( T( K+1, K ).EQ.ZERO ) THEN
                     IF( SELECT( K ) )
     $                  M = M + 1
                  ELSE
                     PAIR = .TRUE.
                     IF( SELECT( K ) .OR. SELECT( K+1 ) )
     $                  M = M + 2
                  END IF
               ELSE
                  IF( SELECT( N ) )
     $               M = M + 1
               END IF
            END IF
   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         N1 = M
         N2 = N - M
         NN = N1*N2
<span class="comment">*</span><span class="comment">
</span>         IF( WANTSP ) THEN
            LWMIN = MAX( 1, 2*NN )
            LIWMIN = MAX( 1, NN )
         ELSE IF( <a name="LSAME.305"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'N'</span> ) ) THEN