LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ztrsna.f
Go to the documentation of this file.
00001 *> \brief \b ZTRSNA
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZTRSNA + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsna.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsna.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsna.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZTRSNA( 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 *       DOUBLE PRECISION   RWORK( * ), S( * ), SEP( * )
00032 *       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00033 *      $                   WORK( LDWORK, * )
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> ZTRSNA 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*16 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*16 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 *>          ZHSEIN or ZTREVC.
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*16 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 *>          ZHSEIN or ZTREVC.
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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 DOUBLE PRECISION 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 complex16OTHERcomputational
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 ZTRSNA( 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       DOUBLE PRECISION   RWORK( * ), S( * ), SEP( * )
00264       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00265      $                   WORK( LDWORK, * )
00266 *     ..
00267 *
00268 *  =====================================================================
00269 *
00270 *     .. Parameters ..
00271       DOUBLE PRECISION   ZERO, ONE
00272       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D0+0 )
00273 *     ..
00274 *     .. Local Scalars ..
00275       LOGICAL            SOMCON, WANTBH, WANTS, WANTSP
00276       CHARACTER          NORMIN
00277       INTEGER            I, IERR, IX, J, K, KASE, KS
00278       DOUBLE PRECISION   BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
00279      $                   XNORM
00280       COMPLEX*16         CDUM, PROD
00281 *     ..
00282 *     .. Local Arrays ..
00283       INTEGER            ISAVE( 3 )
00284       COMPLEX*16         DUMMY( 1 )
00285 *     ..
00286 *     .. External Functions ..
00287       LOGICAL            LSAME
00288       INTEGER            IZAMAX
00289       DOUBLE PRECISION   DLAMCH, DZNRM2
00290       COMPLEX*16         ZDOTC
00291       EXTERNAL           LSAME, IZAMAX, DLAMCH, DZNRM2, ZDOTC
00292 *     ..
00293 *     .. External Subroutines ..
00294       EXTERNAL           XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC
00295 *     ..
00296 *     .. Intrinsic Functions ..
00297       INTRINSIC          ABS, DBLE, DIMAG, MAX
00298 *     ..
00299 *     .. Statement Functions ..
00300       DOUBLE PRECISION   CABS1
00301 *     ..
00302 *     .. Statement Function definitions ..
00303       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
00304 *     ..
00305 *     .. Executable Statements ..
00306 *
00307 *     Decode and test the input parameters
00308 *
00309       WANTBH = LSAME( JOB, 'B' )
00310       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00311       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
00312 *
00313       SOMCON = LSAME( HOWMNY, 'S' )
00314 *
00315 *     Set M to the number of eigenpairs for which condition numbers are
00316 *     to be computed.
00317 *
00318       IF( SOMCON ) THEN
00319          M = 0
00320          DO 10 J = 1, N
00321             IF( SELECT( J ) )
00322      $         M = M + 1
00323    10    CONTINUE
00324       ELSE
00325          M = N
00326       END IF
00327 *
00328       INFO = 0
00329       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
00330          INFO = -1
00331       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
00332          INFO = -2
00333       ELSE IF( N.LT.0 ) THEN
00334          INFO = -4
00335       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00336          INFO = -6
00337       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
00338          INFO = -8
00339       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
00340          INFO = -10
00341       ELSE IF( MM.LT.M ) THEN
00342          INFO = -13
00343       ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
00344          INFO = -16
00345       END IF
00346       IF( INFO.NE.0 ) THEN
00347          CALL XERBLA( 'ZTRSNA', -INFO )
00348          RETURN
00349       END IF
00350 *
00351 *     Quick return if possible
00352 *
00353       IF( N.EQ.0 )
00354      $   RETURN
00355 *
00356       IF( N.EQ.1 ) THEN
00357          IF( SOMCON ) THEN
00358             IF( .NOT.SELECT( 1 ) )
00359      $         RETURN
00360          END IF
00361          IF( WANTS )
00362      $      S( 1 ) = ONE
00363          IF( WANTSP )
00364      $      SEP( 1 ) = ABS( T( 1, 1 ) )
00365          RETURN
00366       END IF
00367 *
00368 *     Get machine constants
00369 *
00370       EPS = DLAMCH( 'P' )
00371       SMLNUM = DLAMCH( 'S' ) / EPS
00372       BIGNUM = ONE / SMLNUM
00373       CALL DLABAD( SMLNUM, BIGNUM )
00374 *
00375       KS = 1
00376       DO 50 K = 1, N
00377 *
00378          IF( SOMCON ) THEN
00379             IF( .NOT.SELECT( K ) )
00380      $         GO TO 50
00381          END IF
00382 *
00383          IF( WANTS ) THEN
00384 *
00385 *           Compute the reciprocal condition number of the k-th
00386 *           eigenvalue.
00387 *
00388             PROD = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
00389             RNRM = DZNRM2( N, VR( 1, KS ), 1 )
00390             LNRM = DZNRM2( N, VL( 1, KS ), 1 )
00391             S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
00392 *
00393          END IF
00394 *
00395          IF( WANTSP ) THEN
00396 *
00397 *           Estimate the reciprocal condition number of the k-th
00398 *           eigenvector.
00399 *
00400 *           Copy the matrix T to the array WORK and swap the k-th
00401 *           diagonal element to the (1,1) position.
00402 *
00403             CALL ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
00404             CALL ZTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
00405 *
00406 *           Form  C = T22 - lambda*I in WORK(2:N,2:N).
00407 *
00408             DO 20 I = 2, N
00409                WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
00410    20       CONTINUE
00411 *
00412 *           Estimate a lower bound for the 1-norm of inv(C**H). The 1st
00413 *           and (N+1)th columns of WORK are used to store work vectors.
00414 *
00415             SEP( KS ) = ZERO
00416             EST = ZERO
00417             KASE = 0
00418             NORMIN = 'N'
00419    30       CONTINUE
00420             CALL ZLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
00421 *
00422             IF( KASE.NE.0 ) THEN
00423                IF( KASE.EQ.1 ) THEN
00424 *
00425 *                 Solve C**H*x = scale*b
00426 *
00427                   CALL ZLATRS( 'Upper', 'Conjugate transpose',
00428      $                         'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
00429      $                         LDWORK, WORK, SCALE, RWORK, IERR )
00430                ELSE
00431 *
00432 *                 Solve C*x = scale*b
00433 *
00434                   CALL ZLATRS( 'Upper', 'No transpose', 'Nonunit',
00435      $                         NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
00436      $                         SCALE, RWORK, IERR )
00437                END IF
00438                NORMIN = 'Y'
00439                IF( SCALE.NE.ONE ) THEN
00440 *
00441 *                 Multiply by 1/SCALE if doing so will not cause
00442 *                 overflow.
00443 *
00444                   IX = IZAMAX( N-1, WORK, 1 )
00445                   XNORM = CABS1( WORK( IX, 1 ) )
00446                   IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
00447      $               GO TO 40
00448                   CALL ZDRSCL( N, SCALE, WORK, 1 )
00449                END IF
00450                GO TO 30
00451             END IF
00452 *
00453             SEP( KS ) = ONE / MAX( EST, SMLNUM )
00454          END IF
00455 *
00456    40    CONTINUE
00457          KS = KS + 1
00458    50 CONTINUE
00459       RETURN
00460 *
00461 *     End of ZTRSNA
00462 *
00463       END
 All Files Functions