LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ctrsna.f
Go to the documentation of this file.
00001 *> \brief \b CTRSNA
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CTRSNA + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrsna.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrsna.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrsna.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00022 *                          LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
00023 *                          INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          HOWMNY, JOB
00027 *       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       LOGICAL            SELECT( * )
00031 *       REAL               RWORK( * ), S( * ), SEP( * )
00032 *       COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00033 *      $                   WORK( LDWORK, * )
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> CTRSNA estimates reciprocal condition numbers for specified
00043 *> eigenvalues and/or right eigenvectors of a complex upper triangular
00044 *> matrix T (or of any matrix Q*T*Q**H with Q unitary).
00045 *> \endverbatim
00046 *
00047 *  Arguments:
00048 *  ==========
00049 *
00050 *> \param[in] JOB
00051 *> \verbatim
00052 *>          JOB is CHARACTER*1
00053 *>          Specifies whether condition numbers are required for
00054 *>          eigenvalues (S) or eigenvectors (SEP):
00055 *>          = 'E': for eigenvalues only (S);
00056 *>          = 'V': for eigenvectors only (SEP);
00057 *>          = 'B': for both eigenvalues and eigenvectors (S and SEP).
00058 *> \endverbatim
00059 *>
00060 *> \param[in] HOWMNY
00061 *> \verbatim
00062 *>          HOWMNY is CHARACTER*1
00063 *>          = 'A': compute condition numbers for all eigenpairs;
00064 *>          = 'S': compute condition numbers for selected eigenpairs
00065 *>                 specified by the array SELECT.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] SELECT
00069 *> \verbatim
00070 *>          SELECT is LOGICAL array, dimension (N)
00071 *>          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
00072 *>          condition numbers are required. To select condition numbers
00073 *>          for the j-th eigenpair, SELECT(j) must be set to .TRUE..
00074 *>          If HOWMNY = 'A', SELECT is not referenced.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] N
00078 *> \verbatim
00079 *>          N is INTEGER
00080 *>          The order of the matrix T. N >= 0.
00081 *> \endverbatim
00082 *>
00083 *> \param[in] T
00084 *> \verbatim
00085 *>          T is COMPLEX array, dimension (LDT,N)
00086 *>          The upper triangular matrix T.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] LDT
00090 *> \verbatim
00091 *>          LDT is INTEGER
00092 *>          The leading dimension of the array T. LDT >= max(1,N).
00093 *> \endverbatim
00094 *>
00095 *> \param[in] VL
00096 *> \verbatim
00097 *>          VL is COMPLEX array, dimension (LDVL,M)
00098 *>          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
00099 *>          (or of any Q*T*Q**H with Q unitary), corresponding to the
00100 *>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
00101 *>          must be stored in consecutive columns of VL, as returned by
00102 *>          CHSEIN or CTREVC.
00103 *>          If JOB = 'V', VL is not referenced.
00104 *> \endverbatim
00105 *>
00106 *> \param[in] LDVL
00107 *> \verbatim
00108 *>          LDVL is INTEGER
00109 *>          The leading dimension of the array VL.
00110 *>          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
00111 *> \endverbatim
00112 *>
00113 *> \param[in] VR
00114 *> \verbatim
00115 *>          VR is COMPLEX array, dimension (LDVR,M)
00116 *>          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
00117 *>          (or of any Q*T*Q**H with Q unitary), corresponding to the
00118 *>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
00119 *>          must be stored in consecutive columns of VR, as returned by
00120 *>          CHSEIN or CTREVC.
00121 *>          If JOB = 'V', VR is not referenced.
00122 *> \endverbatim
00123 *>
00124 *> \param[in] LDVR
00125 *> \verbatim
00126 *>          LDVR is INTEGER
00127 *>          The leading dimension of the array VR.
00128 *>          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
00129 *> \endverbatim
00130 *>
00131 *> \param[out] S
00132 *> \verbatim
00133 *>          S is REAL array, dimension (MM)
00134 *>          If JOB = 'E' or 'B', the reciprocal condition numbers of the
00135 *>          selected eigenvalues, stored in consecutive elements of the
00136 *>          array. Thus S(j), SEP(j), and the j-th columns of VL and VR
00137 *>          all correspond to the same eigenpair (but not in general the
00138 *>          j-th eigenpair, unless all eigenpairs are selected).
00139 *>          If JOB = 'V', S is not referenced.
00140 *> \endverbatim
00141 *>
00142 *> \param[out] SEP
00143 *> \verbatim
00144 *>          SEP is REAL array, dimension (MM)
00145 *>          If JOB = 'V' or 'B', the estimated reciprocal condition
00146 *>          numbers of the selected eigenvectors, stored in consecutive
00147 *>          elements of the array.
00148 *>          If JOB = 'E', SEP is not referenced.
00149 *> \endverbatim
00150 *>
00151 *> \param[in] MM
00152 *> \verbatim
00153 *>          MM is INTEGER
00154 *>          The number of elements in the arrays S (if JOB = 'E' or 'B')
00155 *>           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
00156 *> \endverbatim
00157 *>
00158 *> \param[out] M
00159 *> \verbatim
00160 *>          M is INTEGER
00161 *>          The number of elements of the arrays S and/or SEP actually
00162 *>          used to store the estimated condition numbers.
00163 *>          If HOWMNY = 'A', M is set to N.
00164 *> \endverbatim
00165 *>
00166 *> \param[out] WORK
00167 *> \verbatim
00168 *>          WORK is COMPLEX array, dimension (LDWORK,N+6)
00169 *>          If JOB = 'E', WORK is not referenced.
00170 *> \endverbatim
00171 *>
00172 *> \param[in] LDWORK
00173 *> \verbatim
00174 *>          LDWORK is INTEGER
00175 *>          The leading dimension of the array WORK.
00176 *>          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
00177 *> \endverbatim
00178 *>
00179 *> \param[out] RWORK
00180 *> \verbatim
00181 *>          RWORK is REAL array, dimension (N)
00182 *>          If JOB = 'E', RWORK is not referenced.
00183 *> \endverbatim
00184 *>
00185 *> \param[out] INFO
00186 *> \verbatim
00187 *>          INFO is INTEGER
00188 *>          = 0: successful exit
00189 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00190 *> \endverbatim
00191 *
00192 *  Authors:
00193 *  ========
00194 *
00195 *> \author Univ. of Tennessee 
00196 *> \author Univ. of California Berkeley 
00197 *> \author Univ. of Colorado Denver 
00198 *> \author NAG Ltd. 
00199 *
00200 *> \date November 2011
00201 *
00202 *> \ingroup complexOTHERcomputational
00203 *
00204 *> \par Further Details:
00205 *  =====================
00206 *>
00207 *> \verbatim
00208 *>
00209 *>  The reciprocal of the condition number of an eigenvalue lambda is
00210 *>  defined as
00211 *>
00212 *>          S(lambda) = |v**H*u| / (norm(u)*norm(v))
00213 *>
00214 *>  where u and v are the right and left eigenvectors of T corresponding
00215 *>  to lambda; v**H denotes the conjugate transpose of v, and norm(u)
00216 *>  denotes the Euclidean norm. These reciprocal condition numbers always
00217 *>  lie between zero (very badly conditioned) and one (very well
00218 *>  conditioned). If n = 1, S(lambda) is defined to be 1.
00219 *>
00220 *>  An approximate error bound for a computed eigenvalue W(i) is given by
00221 *>
00222 *>                      EPS * norm(T) / S(i)
00223 *>
00224 *>  where EPS is the machine precision.
00225 *>
00226 *>  The reciprocal of the condition number of the right eigenvector u
00227 *>  corresponding to lambda is defined as follows. Suppose
00228 *>
00229 *>              T = ( lambda  c  )
00230 *>                  (   0    T22 )
00231 *>
00232 *>  Then the reciprocal condition number is
00233 *>
00234 *>          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
00235 *>
00236 *>  where sigma-min denotes the smallest singular value. We approximate
00237 *>  the smallest singular value by the reciprocal of an estimate of the
00238 *>  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
00239 *>  defined to be abs(T(1,1)).
00240 *>
00241 *>  An approximate error bound for a computed right eigenvector VR(i)
00242 *>  is given by
00243 *>
00244 *>                      EPS * norm(T) / SEP(i)
00245 *> \endverbatim
00246 *>
00247 *  =====================================================================
00248       SUBROUTINE CTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00249      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
00250      $                   INFO )
00251 *
00252 *  -- LAPACK computational routine (version 3.4.0) --
00253 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00254 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00255 *     November 2011
00256 *
00257 *     .. Scalar Arguments ..
00258       CHARACTER          HOWMNY, JOB
00259       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
00260 *     ..
00261 *     .. Array Arguments ..
00262       LOGICAL            SELECT( * )
00263       REAL               RWORK( * ), S( * ), SEP( * )
00264       COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00265      $                   WORK( LDWORK, * )
00266 *     ..
00267 *
00268 *  =====================================================================
00269 *
00270 *     .. Parameters ..
00271       REAL               ZERO, ONE
00272       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0+0 )
00273 *     ..
00274 *     .. Local Scalars ..
00275       LOGICAL            SOMCON, WANTBH, WANTS, WANTSP
00276       CHARACTER          NORMIN
00277       INTEGER            I, IERR, IX, J, K, KASE, KS
00278       REAL               BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
00279      $                   XNORM
00280       COMPLEX            CDUM, PROD
00281 *     ..
00282 *     .. Local Arrays ..
00283       INTEGER            ISAVE( 3 )
00284       COMPLEX            DUMMY( 1 )
00285 *     ..
00286 *     .. External Functions ..
00287       LOGICAL            LSAME
00288       INTEGER            ICAMAX
00289       REAL               SCNRM2, SLAMCH
00290       COMPLEX            CDOTC
00291       EXTERNAL           LSAME, ICAMAX, SCNRM2, SLAMCH, CDOTC
00292 *     ..
00293 *     .. External Subroutines ..
00294       EXTERNAL           CLACN2, CLACPY, CLATRS, CSRSCL, CTREXC, SLABAD,
00295      $                   XERBLA
00296 *     ..
00297 *     .. Intrinsic Functions ..
00298       INTRINSIC          ABS, AIMAG, MAX, REAL
00299 *     ..
00300 *     .. Statement Functions ..
00301       REAL               CABS1
00302 *     ..
00303 *     .. Statement Function definitions ..
00304       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
00305 *     ..
00306 *     .. Executable Statements ..
00307 *
00308 *     Decode and test the input parameters
00309 *
00310       WANTBH = LSAME( JOB, 'B' )
00311       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00312       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
00313 *
00314       SOMCON = LSAME( HOWMNY, 'S' )
00315 *
00316 *     Set M to the number of eigenpairs for which condition numbers are
00317 *     to be computed.
00318 *
00319       IF( SOMCON ) THEN
00320          M = 0
00321          DO 10 J = 1, N
00322             IF( SELECT( J ) )
00323      $         M = M + 1
00324    10    CONTINUE
00325       ELSE
00326          M = N
00327       END IF
00328 *
00329       INFO = 0
00330       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
00331          INFO = -1
00332       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
00333          INFO = -2
00334       ELSE IF( N.LT.0 ) THEN
00335          INFO = -4
00336       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00337          INFO = -6
00338       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
00339          INFO = -8
00340       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
00341          INFO = -10
00342       ELSE IF( MM.LT.M ) THEN
00343          INFO = -13
00344       ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
00345          INFO = -16
00346       END IF
00347       IF( INFO.NE.0 ) THEN
00348          CALL XERBLA( 'CTRSNA', -INFO )
00349          RETURN
00350       END IF
00351 *
00352 *     Quick return if possible
00353 *
00354       IF( N.EQ.0 )
00355      $   RETURN
00356 *
00357       IF( N.EQ.1 ) THEN
00358          IF( SOMCON ) THEN
00359             IF( .NOT.SELECT( 1 ) )
00360      $         RETURN
00361          END IF
00362          IF( WANTS )
00363      $      S( 1 ) = ONE
00364          IF( WANTSP )
00365      $      SEP( 1 ) = ABS( T( 1, 1 ) )
00366          RETURN
00367       END IF
00368 *
00369 *     Get machine constants
00370 *
00371       EPS = SLAMCH( 'P' )
00372       SMLNUM = SLAMCH( 'S' ) / EPS
00373       BIGNUM = ONE / SMLNUM
00374       CALL SLABAD( SMLNUM, BIGNUM )
00375 *
00376       KS = 1
00377       DO 50 K = 1, N
00378 *
00379          IF( SOMCON ) THEN
00380             IF( .NOT.SELECT( K ) )
00381      $         GO TO 50
00382          END IF
00383 *
00384          IF( WANTS ) THEN
00385 *
00386 *           Compute the reciprocal condition number of the k-th
00387 *           eigenvalue.
00388 *
00389             PROD = CDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
00390             RNRM = SCNRM2( N, VR( 1, KS ), 1 )
00391             LNRM = SCNRM2( N, VL( 1, KS ), 1 )
00392             S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
00393 *
00394          END IF
00395 *
00396          IF( WANTSP ) THEN
00397 *
00398 *           Estimate the reciprocal condition number of the k-th
00399 *           eigenvector.
00400 *
00401 *           Copy the matrix T to the array WORK and swap the k-th
00402 *           diagonal element to the (1,1) position.
00403 *
00404             CALL CLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
00405             CALL CTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
00406 *
00407 *           Form  C = T22 - lambda*I in WORK(2:N,2:N).
00408 *
00409             DO 20 I = 2, N
00410                WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
00411    20       CONTINUE
00412 *
00413 *           Estimate a lower bound for the 1-norm of inv(C**H). The 1st
00414 *           and (N+1)th columns of WORK are used to store work vectors.
00415 *
00416             SEP( KS ) = ZERO
00417             EST = ZERO
00418             KASE = 0
00419             NORMIN = 'N'
00420    30       CONTINUE
00421             CALL CLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
00422 *
00423             IF( KASE.NE.0 ) THEN
00424                IF( KASE.EQ.1 ) THEN
00425 *
00426 *                 Solve C**H*x = scale*b
00427 *
00428                   CALL CLATRS( 'Upper', 'Conjugate transpose',
00429      $                         'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
00430      $                         LDWORK, WORK, SCALE, RWORK, IERR )
00431                ELSE
00432 *
00433 *                 Solve C*x = scale*b
00434 *
00435                   CALL CLATRS( 'Upper', 'No transpose', 'Nonunit',
00436      $                         NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
00437      $                         SCALE, RWORK, IERR )
00438                END IF
00439                NORMIN = 'Y'
00440                IF( SCALE.NE.ONE ) THEN
00441 *
00442 *                 Multiply by 1/SCALE if doing so will not cause
00443 *                 overflow.
00444 *
00445                   IX = ICAMAX( N-1, WORK, 1 )
00446                   XNORM = CABS1( WORK( IX, 1 ) )
00447                   IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
00448      $               GO TO 40
00449                   CALL CSRSCL( N, SCALE, WORK, 1 )
00450                END IF
00451                GO TO 30
00452             END IF
00453 *
00454             SEP( KS ) = ONE / MAX( EST, SMLNUM )
00455          END IF
00456 *
00457    40    CONTINUE
00458          KS = KS + 1
00459    50 CONTINUE
00460       RETURN
00461 *
00462 *     End of CTRSNA
00463 *
00464       END
 All Files Functions