SRC/ztrsna.f

   1 *> \brief \b ZTRSNA
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZTRSNA + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsna.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsna.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsna.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
  22 *                          LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
  23 *                          INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          HOWMNY, JOB
  27 *       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       LOGICAL            SELECT( * )
  31 *       DOUBLE PRECISION   RWORK( * ), S( * ), SEP( * )
  32 *       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
  33 *      $                   WORK( LDWORK, * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> ZTRSNA estimates reciprocal condition numbers for specified
  43 *> eigenvalues and/or right eigenvectors of a complex upper triangular
  44 *> matrix T (or of any matrix Q*T*Q**H with Q unitary).
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] JOB
  51 *> \verbatim
  52 *>          JOB is CHARACTER*1
  53 *>          Specifies whether condition numbers are required for
  54 *>          eigenvalues (S) or eigenvectors (SEP):
  55 *>          = 'E': for eigenvalues only (S);
  56 *>          = 'V': for eigenvectors only (SEP);
  57 *>          = 'B': for both eigenvalues and eigenvectors (S and SEP).
  58 *> \endverbatim
  59 *>
  60 *> \param[in] HOWMNY
  61 *> \verbatim
  62 *>          HOWMNY is CHARACTER*1
  63 *>          = 'A': compute condition numbers for all eigenpairs;
  64 *>          = 'S': compute condition numbers for selected eigenpairs
  65 *>                 specified by the array SELECT.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] SELECT
  69 *> \verbatim
  70 *>          SELECT is LOGICAL array, dimension (N)
  71 *>          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
  72 *>          condition numbers are required. To select condition numbers
  73 *>          for the j-th eigenpair, SELECT(j) must be set to .TRUE..
  74 *>          If HOWMNY = 'A', SELECT is not referenced.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] N
  78 *> \verbatim
  79 *>          N is INTEGER
  80 *>          The order of the matrix T. N >= 0.
  81 *> \endverbatim
  82 *>
  83 *> \param[in] T
  84 *> \verbatim
  85 *>          T is COMPLEX*16 array, dimension (LDT,N)
  86 *>          The upper triangular matrix T.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] LDT
  90 *> \verbatim
  91 *>          LDT is INTEGER
  92 *>          The leading dimension of the array T. LDT >= max(1,N).
  93 *> \endverbatim
  94 *>
  95 *> \param[in] VL
  96 *> \verbatim
  97 *>          VL is COMPLEX*16 array, dimension (LDVL,M)
  98 *>          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
  99 *>          (or of any Q*T*Q**H with Q unitary), corresponding to the
 100 *>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
 101 *>          must be stored in consecutive columns of VL, as returned by
 102 *>          ZHSEIN or ZTREVC.
 103 *>          If JOB = 'V', VL is not referenced.
 104 *> \endverbatim
 105 *>
 106 *> \param[in] LDVL
 107 *> \verbatim
 108 *>          LDVL is INTEGER
 109 *>          The leading dimension of the array VL.
 110 *>          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] VR
 114 *> \verbatim
 115 *>          VR is COMPLEX*16 array, dimension (LDVR,M)
 116 *>          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
 117 *>          (or of any Q*T*Q**H with Q unitary), corresponding to the
 118 *>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
 119 *>          must be stored in consecutive columns of VR, as returned by
 120 *>          ZHSEIN or ZTREVC.
 121 *>          If JOB = 'V', VR is not referenced.
 122 *> \endverbatim
 123 *>
 124 *> \param[in] LDVR
 125 *> \verbatim
 126 *>          LDVR is INTEGER
 127 *>          The leading dimension of the array VR.
 128 *>          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
 129 *> \endverbatim
 130 *>
 131 *> \param[out] S
 132 *> \verbatim
 133 *>          S is DOUBLE PRECISION array, dimension (MM)
 134 *>          If JOB = 'E' or 'B', the reciprocal condition numbers of the
 135 *>          selected eigenvalues, stored in consecutive elements of the
 136 *>          array. Thus S(j), SEP(j), and the j-th columns of VL and VR
 137 *>          all correspond to the same eigenpair (but not in general the
 138 *>          j-th eigenpair, unless all eigenpairs are selected).
 139 *>          If JOB = 'V', S is not referenced.
 140 *> \endverbatim
 141 *>
 142 *> \param[out] SEP
 143 *> \verbatim
 144 *>          SEP is DOUBLE PRECISION array, dimension (MM)
 145 *>          If JOB = 'V' or 'B', the estimated reciprocal condition
 146 *>          numbers of the selected eigenvectors, stored in consecutive
 147 *>          elements of the array.
 148 *>          If JOB = 'E', SEP is not referenced.
 149 *> \endverbatim
 150 *>
 151 *> \param[in] MM
 152 *> \verbatim
 153 *>          MM is INTEGER
 154 *>          The number of elements in the arrays S (if JOB = 'E' or 'B')
 155 *>           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
 156 *> \endverbatim
 157 *>
 158 *> \param[out] M
 159 *> \verbatim
 160 *>          M is INTEGER
 161 *>          The number of elements of the arrays S and/or SEP actually
 162 *>          used to store the estimated condition numbers.
 163 *>          If HOWMNY = 'A', M is set to N.
 164 *> \endverbatim
 165 *>
 166 *> \param[out] WORK
 167 *> \verbatim
 168 *>          WORK is COMPLEX*16 array, dimension (LDWORK,N+6)
 169 *>          If JOB = 'E', WORK is not referenced.
 170 *> \endverbatim
 171 *>
 172 *> \param[in] LDWORK
 173 *> \verbatim
 174 *>          LDWORK is INTEGER
 175 *>          The leading dimension of the array WORK.
 176 *>          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
 177 *> \endverbatim
 178 *>
 179 *> \param[out] RWORK
 180 *> \verbatim
 181 *>          RWORK is DOUBLE PRECISION array, dimension (N)
 182 *>          If JOB = 'E', RWORK is not referenced.
 183 *> \endverbatim
 184 *>
 185 *> \param[out] INFO
 186 *> \verbatim
 187 *>          INFO is INTEGER
 188 *>          = 0: successful exit
 189 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 190 *> \endverbatim
 191 *
 192 *  Authors:
 193 *  ========
 194 *
 195 *> \author Univ. of Tennessee
 196 *> \author Univ. of California Berkeley
 197 *> \author Univ. of Colorado Denver
 198 *> \author NAG Ltd.
 199 *
 200 *> \date November 2011
 201 *
 202 *> \ingroup complex16OTHERcomputational
 203 *
 204 *> \par Further Details:
 205 *  =====================
 206 *>
 207 *> \verbatim
 208 *>
 209 *>  The reciprocal of the condition number of an eigenvalue lambda is
 210 *>  defined as
 211 *>
 212 *>          S(lambda) = |v**H*u| / (norm(u)*norm(v))
 213 *>
 214 *>  where u and v are the right and left eigenvectors of T corresponding
 215 *>  to lambda; v**H denotes the conjugate transpose of v, and norm(u)
 216 *>  denotes the Euclidean norm. These reciprocal condition numbers always
 217 *>  lie between zero (very badly conditioned) and one (very well
 218 *>  conditioned). If n = 1, S(lambda) is defined to be 1.
 219 *>
 220 *>  An approximate error bound for a computed eigenvalue W(i) is given by
 221 *>
 222 *>                      EPS * norm(T) / S(i)
 223 *>
 224 *>  where EPS is the machine precision.
 225 *>
 226 *>  The reciprocal of the condition number of the right eigenvector u
 227 *>  corresponding to lambda is defined as follows. Suppose
 228 *>
 229 *>              T = ( lambda  c  )
 230 *>                  (   0    T22 )
 231 *>
 232 *>  Then the reciprocal condition number is
 233 *>
 234 *>          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
 235 *>
 236 *>  where sigma-min denotes the smallest singular value. We approximate
 237 *>  the smallest singular value by the reciprocal of an estimate of the
 238 *>  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
 239 *>  defined to be abs(T(1,1)).
 240 *>
 241 *>  An approximate error bound for a computed right eigenvector VR(i)
 242 *>  is given by
 243 *>
 244 *>                      EPS * norm(T) / SEP(i)
 245 *> \endverbatim
 246 *>
 247 *  =====================================================================
 248       SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
 249      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
 250      $                   INFO )
 251 *
 252 *  -- LAPACK computational routine (version 3.4.0) --
 253 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 254 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 255 *     November 2011
 256 *
 257 *     .. Scalar Arguments ..
 258       CHARACTER          HOWMNY, JOB
 259       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
 260 *     ..
 261 *     .. Array Arguments ..
 262       LOGICAL            SELECT( * )
 263       DOUBLE PRECISION   RWORK( * ), S( * ), SEP( * )
 264       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
 265      $                   WORK( LDWORK, * )
 266 *     ..
 267 *
 268 *  =====================================================================
 269 *
 270 *     .. Parameters ..
 271       DOUBLE PRECISION   ZERO, ONE
 272       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D0+0 )
 273 *     ..
 274 *     .. Local Scalars ..
 275       LOGICAL            SOMCON, WANTBH, WANTS, WANTSP
 276       CHARACTER          NORMIN
 277       INTEGER            I, IERR, IX, J, K, KASE, KS
 278       DOUBLE PRECISION   BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
 279      $                   XNORM
 280       COMPLEX*16         CDUM, PROD
 281 *     ..
 282 *     .. Local Arrays ..
 283       INTEGER            ISAVE( 3 )
 284       COMPLEX*16         DUMMY( 1 )
 285 *     ..
 286 *     .. External Functions ..
 287       LOGICAL            LSAME
 288       INTEGER            IZAMAX
 289       DOUBLE PRECISION   DLAMCH, DZNRM2
 290       COMPLEX*16         ZDOTC
 291       EXTERNAL           LSAME, IZAMAX, DLAMCH, DZNRM2, ZDOTC
 292 *     ..
 293 *     .. External Subroutines ..
 294       EXTERNAL           XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC
 295 *     ..
 296 *     .. Intrinsic Functions ..
 297       INTRINSIC          ABS, DBLE, DIMAG, MAX
 298 *     ..
 299 *     .. Statement Functions ..
 300       DOUBLE PRECISION   CABS1
 301 *     ..
 302 *     .. Statement Function definitions ..
 303       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
 304 *     ..
 305 *     .. Executable Statements ..
 306 *
 307 *     Decode and test the input parameters
 308 *
 309       WANTBH = LSAME( JOB, 'B' )
 310       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
 311       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
 312 *
 313       SOMCON = LSAME( HOWMNY, 'S' )
 314 *
 315 *     Set M to the number of eigenpairs for which condition numbers are
 316 *     to be computed.
 317 *
 318       IF( SOMCON ) THEN
 319          M = 0
 320          DO 10 J = 1, N
 321             IF( SELECT( J ) )
 322      $         M = M + 1
 323    10    CONTINUE
 324       ELSE
 325          M = N
 326       END IF
 327 *
 328       INFO = 0
 329       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
 330          INFO = -1
 331       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
 332          INFO = -2
 333       ELSE IF( N.LT.0 ) THEN
 334          INFO = -4
 335       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
 336          INFO = -6
 337       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
 338          INFO = -8
 339       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
 340          INFO = -10
 341       ELSE IF( MM.LT.M ) THEN
 342          INFO = -13
 343       ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
 344          INFO = -16
 345       END IF
 346       IF( INFO.NE.0 ) THEN
 347          CALL XERBLA( 'ZTRSNA', -INFO )
 348          RETURN
 349       END IF
 350 *
 351 *     Quick return if possible
 352 *
 353       IF( N.EQ.0 )
 354      $   RETURN
 355 *
 356       IF( N.EQ.1 ) THEN
 357          IF( SOMCON ) THEN
 358             IF( .NOT.SELECT( 1 ) )
 359      $         RETURN
 360          END IF
 361          IF( WANTS )
 362      $      S( 1 ) = ONE
 363          IF( WANTSP )
 364      $      SEP( 1 ) = ABS( T( 1, 1 ) )
 365          RETURN
 366       END IF
 367 *
 368 *     Get machine constants
 369 *
 370       EPS = DLAMCH( 'P' )
 371       SMLNUM = DLAMCH( 'S' ) / EPS
 372       BIGNUM = ONE / SMLNUM
 373       CALL DLABAD( SMLNUM, BIGNUM )
 374 *
 375       KS = 1
 376       DO 50 K = 1, N
 377 *
 378          IF( SOMCON ) THEN
 379             IF( .NOT.SELECT( K ) )
 380      $         GO TO 50
 381          END IF
 382 *
 383          IF( WANTS ) THEN
 384 *
 385 *           Compute the reciprocal condition number of the k-th
 386 *           eigenvalue.
 387 *
 388             PROD = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
 389             RNRM = DZNRM2( N, VR( 1, KS ), 1 )
 390             LNRM = DZNRM2( N, VL( 1, KS ), 1 )
 391             S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
 392 *
 393          END IF
 394 *
 395          IF( WANTSP ) THEN
 396 *
 397 *           Estimate the reciprocal condition number of the k-th
 398 *           eigenvector.
 399 *
 400 *           Copy the matrix T to the array WORK and swap the k-th
 401 *           diagonal element to the (1,1) position.
 402 *
 403             CALL ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
 404             CALL ZTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
 405 *
 406 *           Form  C = T22 - lambda*I in WORK(2:N,2:N).
 407 *
 408             DO 20 I = 2, N
 409                WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
 410    20       CONTINUE
 411 *
 412 *           Estimate a lower bound for the 1-norm of inv(C**H). The 1st
 413 *           and (N+1)th columns of WORK are used to store work vectors.
 414 *
 415             SEP( KS ) = ZERO
 416             EST = ZERO
 417             KASE = 0
 418             NORMIN = 'N'
 419    30       CONTINUE
 420             CALL ZLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
 421 *
 422             IF( KASE.NE.0 ) THEN
 423                IF( KASE.EQ.1 ) THEN
 424 *
 425 *                 Solve C**H*x = scale*b
 426 *
 427                   CALL ZLATRS( 'Upper', 'Conjugate transpose',
 428      $                         'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
 429      $                         LDWORK, WORK, SCALE, RWORK, IERR )
 430                ELSE
 431 *
 432 *                 Solve C*x = scale*b
 433 *
 434                   CALL ZLATRS( 'Upper', 'No transpose', 'Nonunit',
 435      $                         NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
 436      $                         SCALE, RWORK, IERR )
 437                END IF
 438                NORMIN = 'Y'
 439                IF( SCALE.NE.ONE ) THEN
 440 *
 441 *                 Multiply by 1/SCALE if doing so will not cause
 442 *                 overflow.
 443 *
 444                   IX = IZAMAX( N-1, WORK, 1 )
 445                   XNORM = CABS1( WORK( IX, 1 ) )
 446                   IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
 447      $               GO TO 40
 448                   CALL ZDRSCL( N, SCALE, WORK, 1 )
 449                END IF
 450                GO TO 30
 451             END IF
 452 *
 453             SEP( KS ) = ONE / MAX( EST, SMLNUM )
 454          END IF
 455 *
 456    40    CONTINUE
 457          KS = KS + 1
 458    50 CONTINUE
 459       RETURN
 460 *
 461 *     End of ZTRSNA
 462 *
 463       END