LAPACK 3.3.1
Linear Algebra PACKage

ztrsna.f

Go to the documentation of this file.
00001       SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00002      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
00003      $                   INFO )
00004 *
00005 *  -- LAPACK routine (version 3.3.1) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *  -- April 2011                                                      --
00009 *
00010 *     Modified to call ZLACN2 in place of ZLACON, 10 Feb 03, SJH.
00011 *
00012 *     .. Scalar Arguments ..
00013       CHARACTER          HOWMNY, JOB
00014       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
00015 *     ..
00016 *     .. Array Arguments ..
00017       LOGICAL            SELECT( * )
00018       DOUBLE PRECISION   RWORK( * ), S( * ), SEP( * )
00019       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00020      $                   WORK( LDWORK, * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  ZTRSNA estimates reciprocal condition numbers for specified
00027 *  eigenvalues and/or right eigenvectors of a complex upper triangular
00028 *  matrix T (or of any matrix Q*T*Q**H with Q unitary).
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  JOB     (input) CHARACTER*1
00034 *          Specifies whether condition numbers are required for
00035 *          eigenvalues (S) or eigenvectors (SEP):
00036 *          = 'E': for eigenvalues only (S);
00037 *          = 'V': for eigenvectors only (SEP);
00038 *          = 'B': for both eigenvalues and eigenvectors (S and SEP).
00039 *
00040 *  HOWMNY  (input) CHARACTER*1
00041 *          = 'A': compute condition numbers for all eigenpairs;
00042 *          = 'S': compute condition numbers for selected eigenpairs
00043 *                 specified by the array SELECT.
00044 *
00045 *  SELECT  (input) LOGICAL array, dimension (N)
00046 *          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
00047 *          condition numbers are required. To select condition numbers
00048 *          for the j-th eigenpair, SELECT(j) must be set to .TRUE..
00049 *          If HOWMNY = 'A', SELECT is not referenced.
00050 *
00051 *  N       (input) INTEGER
00052 *          The order of the matrix T. N >= 0.
00053 *
00054 *  T       (input) COMPLEX*16 array, dimension (LDT,N)
00055 *          The upper triangular matrix T.
00056 *
00057 *  LDT     (input) INTEGER
00058 *          The leading dimension of the array T. LDT >= max(1,N).
00059 *
00060 *  VL      (input) COMPLEX*16 array, dimension (LDVL,M)
00061 *          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
00062 *          (or of any Q*T*Q**H with Q unitary), corresponding to the
00063 *          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
00064 *          must be stored in consecutive columns of VL, as returned by
00065 *          ZHSEIN or ZTREVC.
00066 *          If JOB = 'V', VL is not referenced.
00067 *
00068 *  LDVL    (input) INTEGER
00069 *          The leading dimension of the array VL.
00070 *          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
00071 *
00072 *  VR      (input) COMPLEX*16 array, dimension (LDVR,M)
00073 *          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
00074 *          (or of any Q*T*Q**H with Q unitary), corresponding to the
00075 *          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
00076 *          must be stored in consecutive columns of VR, as returned by
00077 *          ZHSEIN or ZTREVC.
00078 *          If JOB = 'V', VR is not referenced.
00079 *
00080 *  LDVR    (input) INTEGER
00081 *          The leading dimension of the array VR.
00082 *          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
00083 *
00084 *  S       (output) DOUBLE PRECISION array, dimension (MM)
00085 *          If JOB = 'E' or 'B', the reciprocal condition numbers of the
00086 *          selected eigenvalues, stored in consecutive elements of the
00087 *          array. Thus S(j), SEP(j), and the j-th columns of VL and VR
00088 *          all correspond to the same eigenpair (but not in general the
00089 *          j-th eigenpair, unless all eigenpairs are selected).
00090 *          If JOB = 'V', S is not referenced.
00091 *
00092 *  SEP     (output) DOUBLE PRECISION array, dimension (MM)
00093 *          If JOB = 'V' or 'B', the estimated reciprocal condition
00094 *          numbers of the selected eigenvectors, stored in consecutive
00095 *          elements of the array.
00096 *          If JOB = 'E', SEP is not referenced.
00097 *
00098 *  MM      (input) INTEGER
00099 *          The number of elements in the arrays S (if JOB = 'E' or 'B')
00100 *           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
00101 *
00102 *  M       (output) INTEGER
00103 *          The number of elements of the arrays S and/or SEP actually
00104 *          used to store the estimated condition numbers.
00105 *          If HOWMNY = 'A', M is set to N.
00106 *
00107 *  WORK    (workspace) COMPLEX*16 array, dimension (LDWORK,N+6)
00108 *          If JOB = 'E', WORK is not referenced.
00109 *
00110 *  LDWORK  (input) INTEGER
00111 *          The leading dimension of the array WORK.
00112 *          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
00113 *
00114 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00115 *          If JOB = 'E', RWORK is not referenced.
00116 *
00117 *  INFO    (output) INTEGER
00118 *          = 0: successful exit
00119 *          < 0: if INFO = -i, the i-th argument had an illegal value
00120 *
00121 *  Further Details
00122 *  ===============
00123 *
00124 *  The reciprocal of the condition number of an eigenvalue lambda is
00125 *  defined as
00126 *
00127 *          S(lambda) = |v**H*u| / (norm(u)*norm(v))
00128 *
00129 *  where u and v are the right and left eigenvectors of T corresponding
00130 *  to lambda; v**H denotes the conjugate transpose of v, and norm(u)
00131 *  denotes the Euclidean norm. These reciprocal condition numbers always
00132 *  lie between zero (very badly conditioned) and one (very well
00133 *  conditioned). If n = 1, S(lambda) is defined to be 1.
00134 *
00135 *  An approximate error bound for a computed eigenvalue W(i) is given by
00136 *
00137 *                      EPS * norm(T) / S(i)
00138 *
00139 *  where EPS is the machine precision.
00140 *
00141 *  The reciprocal of the condition number of the right eigenvector u
00142 *  corresponding to lambda is defined as follows. Suppose
00143 *
00144 *              T = ( lambda  c  )
00145 *                  (   0    T22 )
00146 *
00147 *  Then the reciprocal condition number is
00148 *
00149 *          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
00150 *
00151 *  where sigma-min denotes the smallest singular value. We approximate
00152 *  the smallest singular value by the reciprocal of an estimate of the
00153 *  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
00154 *  defined to be abs(T(1,1)).
00155 *
00156 *  An approximate error bound for a computed right eigenvector VR(i)
00157 *  is given by
00158 *
00159 *                      EPS * norm(T) / SEP(i)
00160 *
00161 *  =====================================================================
00162 *
00163 *     .. Parameters ..
00164       DOUBLE PRECISION   ZERO, ONE
00165       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D0+0 )
00166 *     ..
00167 *     .. Local Scalars ..
00168       LOGICAL            SOMCON, WANTBH, WANTS, WANTSP
00169       CHARACTER          NORMIN
00170       INTEGER            I, IERR, IX, J, K, KASE, KS
00171       DOUBLE PRECISION   BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
00172      $                   XNORM
00173       COMPLEX*16         CDUM, PROD
00174 *     ..
00175 *     .. Local Arrays ..
00176       INTEGER            ISAVE( 3 )
00177       COMPLEX*16         DUMMY( 1 )
00178 *     ..
00179 *     .. External Functions ..
00180       LOGICAL            LSAME
00181       INTEGER            IZAMAX
00182       DOUBLE PRECISION   DLAMCH, DZNRM2
00183       COMPLEX*16         ZDOTC
00184       EXTERNAL           LSAME, IZAMAX, DLAMCH, DZNRM2, ZDOTC
00185 *     ..
00186 *     .. External Subroutines ..
00187       EXTERNAL           XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC
00188 *     ..
00189 *     .. Intrinsic Functions ..
00190       INTRINSIC          ABS, DBLE, DIMAG, MAX
00191 *     ..
00192 *     .. Statement Functions ..
00193       DOUBLE PRECISION   CABS1
00194 *     ..
00195 *     .. Statement Function definitions ..
00196       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
00197 *     ..
00198 *     .. Executable Statements ..
00199 *
00200 *     Decode and test the input parameters
00201 *
00202       WANTBH = LSAME( JOB, 'B' )
00203       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
00204       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
00205 *
00206       SOMCON = LSAME( HOWMNY, 'S' )
00207 *
00208 *     Set M to the number of eigenpairs for which condition numbers are
00209 *     to be computed.
00210 *
00211       IF( SOMCON ) THEN
00212          M = 0
00213          DO 10 J = 1, N
00214             IF( SELECT( J ) )
00215      $         M = M + 1
00216    10    CONTINUE
00217       ELSE
00218          M = N
00219       END IF
00220 *
00221       INFO = 0
00222       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
00223          INFO = -1
00224       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
00225          INFO = -2
00226       ELSE IF( N.LT.0 ) THEN
00227          INFO = -4
00228       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00229          INFO = -6
00230       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
00231          INFO = -8
00232       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
00233          INFO = -10
00234       ELSE IF( MM.LT.M ) THEN
00235          INFO = -13
00236       ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
00237          INFO = -16
00238       END IF
00239       IF( INFO.NE.0 ) THEN
00240          CALL XERBLA( 'ZTRSNA', -INFO )
00241          RETURN
00242       END IF
00243 *
00244 *     Quick return if possible
00245 *
00246       IF( N.EQ.0 )
00247      $   RETURN
00248 *
00249       IF( N.EQ.1 ) THEN
00250          IF( SOMCON ) THEN
00251             IF( .NOT.SELECT( 1 ) )
00252      $         RETURN
00253          END IF
00254          IF( WANTS )
00255      $      S( 1 ) = ONE
00256          IF( WANTSP )
00257      $      SEP( 1 ) = ABS( T( 1, 1 ) )
00258          RETURN
00259       END IF
00260 *
00261 *     Get machine constants
00262 *
00263       EPS = DLAMCH( 'P' )
00264       SMLNUM = DLAMCH( 'S' ) / EPS
00265       BIGNUM = ONE / SMLNUM
00266       CALL DLABAD( SMLNUM, BIGNUM )
00267 *
00268       KS = 1
00269       DO 50 K = 1, N
00270 *
00271          IF( SOMCON ) THEN
00272             IF( .NOT.SELECT( K ) )
00273      $         GO TO 50
00274          END IF
00275 *
00276          IF( WANTS ) THEN
00277 *
00278 *           Compute the reciprocal condition number of the k-th
00279 *           eigenvalue.
00280 *
00281             PROD = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
00282             RNRM = DZNRM2( N, VR( 1, KS ), 1 )
00283             LNRM = DZNRM2( N, VL( 1, KS ), 1 )
00284             S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
00285 *
00286          END IF
00287 *
00288          IF( WANTSP ) THEN
00289 *
00290 *           Estimate the reciprocal condition number of the k-th
00291 *           eigenvector.
00292 *
00293 *           Copy the matrix T to the array WORK and swap the k-th
00294 *           diagonal element to the (1,1) position.
00295 *
00296             CALL ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
00297             CALL ZTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
00298 *
00299 *           Form  C = T22 - lambda*I in WORK(2:N,2:N).
00300 *
00301             DO 20 I = 2, N
00302                WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
00303    20       CONTINUE
00304 *
00305 *           Estimate a lower bound for the 1-norm of inv(C**H). The 1st
00306 *           and (N+1)th columns of WORK are used to store work vectors.
00307 *
00308             SEP( KS ) = ZERO
00309             EST = ZERO
00310             KASE = 0
00311             NORMIN = 'N'
00312    30       CONTINUE
00313             CALL ZLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
00314 *
00315             IF( KASE.NE.0 ) THEN
00316                IF( KASE.EQ.1 ) THEN
00317 *
00318 *                 Solve C**H*x = scale*b
00319 *
00320                   CALL ZLATRS( 'Upper', 'Conjugate transpose',
00321      $                         'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
00322      $                         LDWORK, WORK, SCALE, RWORK, IERR )
00323                ELSE
00324 *
00325 *                 Solve C*x = scale*b
00326 *
00327                   CALL ZLATRS( 'Upper', 'No transpose', 'Nonunit',
00328      $                         NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
00329      $                         SCALE, RWORK, IERR )
00330                END IF
00331                NORMIN = 'Y'
00332                IF( SCALE.NE.ONE ) THEN
00333 *
00334 *                 Multiply by 1/SCALE if doing so will not cause
00335 *                 overflow.
00336 *
00337                   IX = IZAMAX( N-1, WORK, 1 )
00338                   XNORM = CABS1( WORK( IX, 1 ) )
00339                   IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
00340      $               GO TO 40
00341                   CALL ZDRSCL( N, SCALE, WORK, 1 )
00342                END IF
00343                GO TO 30
00344             END IF
00345 *
00346             SEP( KS ) = ONE / MAX( EST, SMLNUM )
00347          END IF
00348 *
00349    40    CONTINUE
00350          KS = KS + 1
00351    50 CONTINUE
00352       RETURN
00353 *
00354 *     End of ZTRSNA
00355 *
00356       END
 All Files Functions