LAPACK 3.3.1
Linear Algebra PACKage

ctrsen.f

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