LAPACK 3.3.1
Linear Algebra PACKage

dtrsen.f

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