LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ztrsen.f
Go to the documentation of this file.
00001 *> \brief \b ZTRSEN
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZTRSEN + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsen.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsen.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsen.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
00022 *                          SEP, WORK, LWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          COMPQ, JOB
00026 *       INTEGER            INFO, LDQ, LDT, LWORK, M, N
00027 *       DOUBLE PRECISION   S, SEP
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       LOGICAL            SELECT( * )
00031 *       COMPLEX*16         Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> ZTRSEN reorders the Schur factorization of a complex matrix
00041 *> A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in
00042 *> the leading positions on the diagonal of the upper triangular matrix
00043 *> T, and the leading columns of Q form an orthonormal basis of the
00044 *> corresponding right invariant subspace.
00045 *>
00046 *> Optionally the routine computes the reciprocal condition numbers of
00047 *> the cluster of eigenvalues and/or the invariant subspace.
00048 *> \endverbatim
00049 *
00050 *  Arguments:
00051 *  ==========
00052 *
00053 *> \param[in] JOB
00054 *> \verbatim
00055 *>          JOB is CHARACTER*1
00056 *>          Specifies whether condition numbers are required for the
00057 *>          cluster of eigenvalues (S) or the invariant subspace (SEP):
00058 *>          = 'N': none;
00059 *>          = 'E': for eigenvalues only (S);
00060 *>          = 'V': for invariant subspace only (SEP);
00061 *>          = 'B': for both eigenvalues and invariant subspace (S and
00062 *>                 SEP).
00063 *> \endverbatim
00064 *>
00065 *> \param[in] COMPQ
00066 *> \verbatim
00067 *>          COMPQ is CHARACTER*1
00068 *>          = 'V': update the matrix Q of Schur vectors;
00069 *>          = 'N': do not update Q.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] SELECT
00073 *> \verbatim
00074 *>          SELECT is LOGICAL array, dimension (N)
00075 *>          SELECT specifies the eigenvalues in the selected cluster. To
00076 *>          select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
00077 *> \endverbatim
00078 *>
00079 *> \param[in] N
00080 *> \verbatim
00081 *>          N is INTEGER
00082 *>          The order of the matrix T. N >= 0.
00083 *> \endverbatim
00084 *>
00085 *> \param[in,out] T
00086 *> \verbatim
00087 *>          T is COMPLEX*16 array, dimension (LDT,N)
00088 *>          On entry, the upper triangular matrix T.
00089 *>          On exit, T is overwritten by the reordered matrix T, with the
00090 *>          selected eigenvalues as the leading diagonal elements.
00091 *> \endverbatim
00092 *>
00093 *> \param[in] LDT
00094 *> \verbatim
00095 *>          LDT is INTEGER
00096 *>          The leading dimension of the array T. LDT >= max(1,N).
00097 *> \endverbatim
00098 *>
00099 *> \param[in,out] Q
00100 *> \verbatim
00101 *>          Q is COMPLEX*16 array, dimension (LDQ,N)
00102 *>          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
00103 *>          On exit, if COMPQ = 'V', Q has been postmultiplied by the
00104 *>          unitary transformation matrix which reorders T; the leading M
00105 *>          columns of Q form an orthonormal basis for the specified
00106 *>          invariant subspace.
00107 *>          If COMPQ = 'N', Q is not referenced.
00108 *> \endverbatim
00109 *>
00110 *> \param[in] LDQ
00111 *> \verbatim
00112 *>          LDQ is INTEGER
00113 *>          The leading dimension of the array Q.
00114 *>          LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
00115 *> \endverbatim
00116 *>
00117 *> \param[out] W
00118 *> \verbatim
00119 *>          W is COMPLEX*16 array, dimension (N)
00120 *>          The reordered eigenvalues of T, in the same order as they
00121 *>          appear on the diagonal of T.
00122 *> \endverbatim
00123 *>
00124 *> \param[out] M
00125 *> \verbatim
00126 *>          M is INTEGER
00127 *>          The dimension of the specified invariant subspace.
00128 *>          0 <= M <= N.
00129 *> \endverbatim
00130 *>
00131 *> \param[out] S
00132 *> \verbatim
00133 *>          S is DOUBLE PRECISION
00134 *>          If JOB = 'E' or 'B', S is a lower bound on the reciprocal
00135 *>          condition number for the selected cluster of eigenvalues.
00136 *>          S cannot underestimate the true reciprocal condition number
00137 *>          by more than a factor of sqrt(N). If M = 0 or N, S = 1.
00138 *>          If JOB = 'N' or 'V', S is not referenced.
00139 *> \endverbatim
00140 *>
00141 *> \param[out] SEP
00142 *> \verbatim
00143 *>          SEP is DOUBLE PRECISION
00144 *>          If JOB = 'V' or 'B', SEP is the estimated reciprocal
00145 *>          condition number of the specified invariant subspace. If
00146 *>          M = 0 or N, SEP = norm(T).
00147 *>          If JOB = 'N' or 'E', SEP is not referenced.
00148 *> \endverbatim
00149 *>
00150 *> \param[out] WORK
00151 *> \verbatim
00152 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
00153 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00154 *> \endverbatim
00155 *>
00156 *> \param[in] LWORK
00157 *> \verbatim
00158 *>          LWORK is INTEGER
00159 *>          The dimension of the array WORK.
00160 *>          If JOB = 'N', LWORK >= 1;
00161 *>          if JOB = 'E', LWORK = max(1,M*(N-M));
00162 *>          if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).
00163 *>
00164 *>          If LWORK = -1, then a workspace query is assumed; the routine
00165 *>          only calculates the optimal size of the WORK array, returns
00166 *>          this value as the first entry of the WORK array, and no error
00167 *>          message related to LWORK is issued by XERBLA.
00168 *> \endverbatim
00169 *>
00170 *> \param[out] INFO
00171 *> \verbatim
00172 *>          INFO is INTEGER
00173 *>          = 0:  successful exit
00174 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00175 *> \endverbatim
00176 *
00177 *  Authors:
00178 *  ========
00179 *
00180 *> \author Univ. of Tennessee 
00181 *> \author Univ. of California Berkeley 
00182 *> \author Univ. of Colorado Denver 
00183 *> \author NAG Ltd. 
00184 *
00185 *> \date November 2011
00186 *
00187 *> \ingroup complex16OTHERcomputational
00188 *
00189 *> \par Further Details:
00190 *  =====================
00191 *>
00192 *> \verbatim
00193 *>
00194 *>  ZTRSEN first collects the selected eigenvalues by computing a unitary
00195 *>  transformation Z to move them to the top left corner of T. In other
00196 *>  words, the selected eigenvalues are the eigenvalues of T11 in:
00197 *>
00198 *>          Z**H * T * Z = ( T11 T12 ) n1
00199 *>                         (  0  T22 ) n2
00200 *>                            n1  n2
00201 *>
00202 *>  where N = n1+n2. The first
00203 *>  n1 columns of Z span the specified invariant subspace of T.
00204 *>
00205 *>  If T has been obtained from the Schur factorization of a matrix
00206 *>  A = Q*T*Q**H, then the reordered Schur factorization of A is given by
00207 *>  A = (Q*Z)*(Z**H*T*Z)*(Q*Z)**H, and the first n1 columns of Q*Z span the
00208 *>  corresponding invariant subspace of A.
00209 *>
00210 *>  The reciprocal condition number of the average of the eigenvalues of
00211 *>  T11 may be returned in S. S lies between 0 (very badly conditioned)
00212 *>  and 1 (very well conditioned). It is computed as follows. First we
00213 *>  compute R so that
00214 *>
00215 *>                         P = ( I  R ) n1
00216 *>                             ( 0  0 ) n2
00217 *>                               n1 n2
00218 *>
00219 *>  is the projector on the invariant subspace associated with T11.
00220 *>  R is the solution of the Sylvester equation:
00221 *>
00222 *>                        T11*R - R*T22 = T12.
00223 *>
00224 *>  Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
00225 *>  the two-norm of M. Then S is computed as the lower bound
00226 *>
00227 *>                      (1 + F-norm(R)**2)**(-1/2)
00228 *>
00229 *>  on the reciprocal of 2-norm(P), the true reciprocal condition number.
00230 *>  S cannot underestimate 1 / 2-norm(P) by more than a factor of
00231 *>  sqrt(N).
00232 *>
00233 *>  An approximate error bound for the computed average of the
00234 *>  eigenvalues of T11 is
00235 *>
00236 *>                         EPS * norm(T) / S
00237 *>
00238 *>  where EPS is the machine precision.
00239 *>
00240 *>  The reciprocal condition number of the right invariant subspace
00241 *>  spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
00242 *>  SEP is defined as the separation of T11 and T22:
00243 *>
00244 *>                     sep( T11, T22 ) = sigma-min( C )
00245 *>
00246 *>  where sigma-min(C) is the smallest singular value of the
00247 *>  n1*n2-by-n1*n2 matrix
00248 *>
00249 *>     C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
00250 *>
00251 *>  I(m) is an m by m identity matrix, and kprod denotes the Kronecker
00252 *>  product. We estimate sigma-min(C) by the reciprocal of an estimate of
00253 *>  the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
00254 *>  cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
00255 *>
00256 *>  When SEP is small, small changes in T can cause large changes in
00257 *>  the invariant subspace. An approximate bound on the maximum angular
00258 *>  error in the computed right invariant subspace is
00259 *>
00260 *>                      EPS * norm(T) / SEP
00261 *> \endverbatim
00262 *>
00263 *  =====================================================================
00264       SUBROUTINE ZTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S,
00265      $                   SEP, WORK, LWORK, INFO )
00266 *
00267 *  -- LAPACK computational routine (version 3.4.0) --
00268 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00269 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00270 *     November 2011
00271 *
00272 *     .. Scalar Arguments ..
00273       CHARACTER          COMPQ, JOB
00274       INTEGER            INFO, LDQ, LDT, LWORK, M, N
00275       DOUBLE PRECISION   S, SEP
00276 *     ..
00277 *     .. Array Arguments ..
00278       LOGICAL            SELECT( * )
00279       COMPLEX*16         Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
00280 *     ..
00281 *
00282 *  =====================================================================
00283 *
00284 *     .. Parameters ..
00285       DOUBLE PRECISION   ZERO, ONE
00286       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00287 *     ..
00288 *     .. Local Scalars ..
00289       LOGICAL            LQUERY, WANTBH, WANTQ, WANTS, WANTSP
00290       INTEGER            IERR, K, KASE, KS, LWMIN, N1, N2, NN
00291       DOUBLE PRECISION   EST, RNORM, SCALE
00292 *     ..
00293 *     .. Local Arrays ..
00294       INTEGER            ISAVE( 3 )
00295       DOUBLE PRECISION   RWORK( 1 )
00296 *     ..
00297 *     .. External Functions ..
00298       LOGICAL            LSAME
00299       DOUBLE PRECISION   ZLANGE
00300       EXTERNAL           LSAME, ZLANGE
00301 *     ..
00302 *     .. External Subroutines ..
00303       EXTERNAL           XERBLA, ZLACN2, ZLACPY, ZTREXC, ZTRSYL
00304 *     ..
00305 *     .. Intrinsic Functions ..
00306       INTRINSIC          MAX, SQRT
00307 *     ..
00308 *     .. Executable Statements ..
00309 *
00310 *     Decode and test the input parameters.
00311 *
00312       WANTBH = LSAME( JOB, 'B' )
00313       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00314       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
00315       WANTQ = LSAME( COMPQ, 'V' )
00316 *
00317 *     Set M to the number of selected eigenvalues.
00318 *
00319       M = 0
00320       DO 10 K = 1, N
00321          IF( SELECT( K ) )
00322      $      M = M + 1
00323    10 CONTINUE
00324 *
00325       N1 = M
00326       N2 = N - M
00327       NN = N1*N2
00328 *
00329       INFO = 0
00330       LQUERY = ( LWORK.EQ.-1 )
00331 *
00332       IF( WANTSP ) THEN
00333          LWMIN = MAX( 1, 2*NN )
00334       ELSE IF( LSAME( JOB, 'N' ) ) THEN
00335          LWMIN = 1
00336       ELSE IF( LSAME( JOB, 'E' ) ) THEN
00337          LWMIN = MAX( 1, NN )
00338       END IF
00339 *
00340       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
00341      $     THEN
00342          INFO = -1
00343       ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
00344          INFO = -2
00345       ELSE IF( N.LT.0 ) THEN
00346          INFO = -4
00347       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00348          INFO = -6
00349       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
00350          INFO = -8
00351       ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00352          INFO = -14
00353       END IF
00354 *
00355       IF( INFO.EQ.0 ) THEN
00356          WORK( 1 ) = LWMIN
00357       END IF
00358 *
00359       IF( INFO.NE.0 ) THEN
00360          CALL XERBLA( 'ZTRSEN', -INFO )
00361          RETURN
00362       ELSE IF( LQUERY ) THEN
00363          RETURN
00364       END IF
00365 *
00366 *     Quick return if possible
00367 *
00368       IF( M.EQ.N .OR. M.EQ.0 ) THEN
00369          IF( WANTS )
00370      $      S = ONE
00371          IF( WANTSP )
00372      $      SEP = ZLANGE( '1', N, N, T, LDT, RWORK )
00373          GO TO 40
00374       END IF
00375 *
00376 *     Collect the selected eigenvalues at the top left corner of T.
00377 *
00378       KS = 0
00379       DO 20 K = 1, N
00380          IF( SELECT( K ) ) THEN
00381             KS = KS + 1
00382 *
00383 *           Swap the K-th eigenvalue to position KS.
00384 *
00385             IF( K.NE.KS )
00386      $         CALL ZTREXC( COMPQ, N, T, LDT, Q, LDQ, K, KS, IERR )
00387          END IF
00388    20 CONTINUE
00389 *
00390       IF( WANTS ) THEN
00391 *
00392 *        Solve the Sylvester equation for R:
00393 *
00394 *           T11*R - R*T22 = scale*T12
00395 *
00396          CALL ZLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
00397          CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
00398      $                LDT, WORK, N1, SCALE, IERR )
00399 *
00400 *        Estimate the reciprocal of the condition number of the cluster
00401 *        of eigenvalues.
00402 *
00403          RNORM = ZLANGE( 'F', N1, N2, WORK, N1, RWORK )
00404          IF( RNORM.EQ.ZERO ) THEN
00405             S = ONE
00406          ELSE
00407             S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
00408      $          SQRT( RNORM ) )
00409          END IF
00410       END IF
00411 *
00412       IF( WANTSP ) THEN
00413 *
00414 *        Estimate sep(T11,T22).
00415 *
00416          EST = ZERO
00417          KASE = 0
00418    30    CONTINUE
00419          CALL ZLACN2( NN, WORK( NN+1 ), WORK, EST, KASE, ISAVE )
00420          IF( KASE.NE.0 ) THEN
00421             IF( KASE.EQ.1 ) THEN
00422 *
00423 *              Solve T11*R - R*T22 = scale*X.
00424 *
00425                CALL ZTRSYL( 'N', 'N', -1, N1, N2, T, LDT,
00426      $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
00427      $                      IERR )
00428             ELSE
00429 *
00430 *              Solve T11**H*R - R*T22**H = scale*X.
00431 *
00432                CALL ZTRSYL( 'C', 'C', -1, N1, N2, T, LDT,
00433      $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
00434      $                      IERR )
00435             END IF
00436             GO TO 30
00437          END IF
00438 *
00439          SEP = SCALE / EST
00440       END IF
00441 *
00442    40 CONTINUE
00443 *
00444 *     Copy reordered eigenvalues to W.
00445 *
00446       DO 50 K = 1, N
00447          W( K ) = T( K, K )
00448    50 CONTINUE
00449 *
00450       WORK( 1 ) = LWMIN
00451 *
00452       RETURN
00453 *
00454 *     End of ZTRSEN
00455 *
00456       END
 All Files Functions