SRC/ctrsna.f

   1 *> \brief \b CTRSNA
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CTRSNA + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrsna.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrsna.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrsna.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CTRSNA( 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 *       REAL               RWORK( * ), S( * ), SEP( * )
  32 *       COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
  33 *      $                   WORK( LDWORK, * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> CTRSNA 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 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 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 *>          CHSEIN or CTREVC.
 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 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 *>          CHSEIN or CTREVC.
 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 REAL 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 REAL 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 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 REAL 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 complexOTHERcomputational
 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 CTRSNA( 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       REAL               RWORK( * ), S( * ), SEP( * )
 264       COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
 265      $                   WORK( LDWORK, * )
 266 *     ..
 267 *
 268 *  =====================================================================
 269 *
 270 *     .. Parameters ..
 271       REAL               ZERO, ONE
 272       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0+0 )
 273 *     ..
 274 *     .. Local Scalars ..
 275       LOGICAL            SOMCON, WANTBH, WANTS, WANTSP
 276       CHARACTER          NORMIN
 277       INTEGER            I, IERR, IX, J, K, KASE, KS
 278       REAL               BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
 279      $                   XNORM
 280       COMPLEX            CDUM, PROD
 281 *     ..
 282 *     .. Local Arrays ..
 283       INTEGER            ISAVE( 3 )
 284       COMPLEX            DUMMY( 1 )
 285 *     ..
 286 *     .. External Functions ..
 287       LOGICAL            LSAME
 288       INTEGER            ICAMAX
 289       REAL               SCNRM2, SLAMCH
 290       COMPLEX            CDOTC
 291       EXTERNAL           LSAME, ICAMAX, SCNRM2, SLAMCH, CDOTC
 292 *     ..
 293 *     .. External Subroutines ..
 294       EXTERNAL           CLACN2, CLACPY, CLATRS, CSRSCL, CTREXC, SLABAD,
 295      $                   XERBLA
 296 *     ..
 297 *     .. Intrinsic Functions ..
 298       INTRINSIC          ABS, AIMAG, MAX, REAL
 299 *     ..
 300 *     .. Statement Functions ..
 301       REAL               CABS1
 302 *     ..
 303 *     .. Statement Function definitions ..
 304       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
 305 *     ..
 306 *     .. Executable Statements ..
 307 *
 308 *     Decode and test the input parameters
 309 *
 310       WANTBH = LSAME( JOB, 'B' )
 311       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
 312       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
 313 *
 314       SOMCON = LSAME( HOWMNY, 'S' )
 315 *
 316 *     Set M to the number of eigenpairs for which condition numbers are
 317 *     to be computed.
 318 *
 319       IF( SOMCON ) THEN
 320          M = 0
 321          DO 10 J = 1, N
 322             IF( SELECT( J ) )
 323      $         M = M + 1
 324    10    CONTINUE
 325       ELSE
 326          M = N
 327       END IF
 328 *
 329       INFO = 0
 330       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
 331          INFO = -1
 332       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
 333          INFO = -2
 334       ELSE IF( N.LT.0 ) THEN
 335          INFO = -4
 336       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
 337          INFO = -6
 338       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
 339          INFO = -8
 340       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
 341          INFO = -10
 342       ELSE IF( MM.LT.M ) THEN
 343          INFO = -13
 344       ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
 345          INFO = -16
 346       END IF
 347       IF( INFO.NE.0 ) THEN
 348          CALL XERBLA( 'CTRSNA', -INFO )
 349          RETURN
 350       END IF
 351 *
 352 *     Quick return if possible
 353 *
 354       IF( N.EQ.0 )
 355      $   RETURN
 356 *
 357       IF( N.EQ.1 ) THEN
 358          IF( SOMCON ) THEN
 359             IF( .NOT.SELECT( 1 ) )
 360      $         RETURN
 361          END IF
 362          IF( WANTS )
 363      $      S( 1 ) = ONE
 364          IF( WANTSP )
 365      $      SEP( 1 ) = ABS( T( 1, 1 ) )
 366          RETURN
 367       END IF
 368 *
 369 *     Get machine constants
 370 *
 371       EPS = SLAMCH( 'P' )
 372       SMLNUM = SLAMCH( 'S' ) / EPS
 373       BIGNUM = ONE / SMLNUM
 374       CALL SLABAD( SMLNUM, BIGNUM )
 375 *
 376       KS = 1
 377       DO 50 K = 1, N
 378 *
 379          IF( SOMCON ) THEN
 380             IF( .NOT.SELECT( K ) )
 381      $         GO TO 50
 382          END IF
 383 *
 384          IF( WANTS ) THEN
 385 *
 386 *           Compute the reciprocal condition number of the k-th
 387 *           eigenvalue.
 388 *
 389             PROD = CDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
 390             RNRM = SCNRM2( N, VR( 1, KS ), 1 )
 391             LNRM = SCNRM2( N, VL( 1, KS ), 1 )
 392             S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
 393 *
 394          END IF
 395 *
 396          IF( WANTSP ) THEN
 397 *
 398 *           Estimate the reciprocal condition number of the k-th
 399 *           eigenvector.
 400 *
 401 *           Copy the matrix T to the array WORK and swap the k-th
 402 *           diagonal element to the (1,1) position.
 403 *
 404             CALL CLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
 405             CALL CTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR )
 406 *
 407 *           Form  C = T22 - lambda*I in WORK(2:N,2:N).
 408 *
 409             DO 20 I = 2, N
 410                WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
 411    20       CONTINUE
 412 *
 413 *           Estimate a lower bound for the 1-norm of inv(C**H). The 1st
 414 *           and (N+1)th columns of WORK are used to store work vectors.
 415 *
 416             SEP( KS ) = ZERO
 417             EST = ZERO
 418             KASE = 0
 419             NORMIN = 'N'
 420    30       CONTINUE
 421             CALL CLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE )
 422 *
 423             IF( KASE.NE.0 ) THEN
 424                IF( KASE.EQ.1 ) THEN
 425 *
 426 *                 Solve C**H*x = scale*b
 427 *
 428                   CALL CLATRS( 'Upper', 'Conjugate transpose',
 429      $                         'Nonunit', NORMIN, N-1, WORK( 2, 2 ),
 430      $                         LDWORK, WORK, SCALE, RWORK, IERR )
 431                ELSE
 432 *
 433 *                 Solve C*x = scale*b
 434 *
 435                   CALL CLATRS( 'Upper', 'No transpose', 'Nonunit',
 436      $                         NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK,
 437      $                         SCALE, RWORK, IERR )
 438                END IF
 439                NORMIN = 'Y'
 440                IF( SCALE.NE.ONE ) THEN
 441 *
 442 *                 Multiply by 1/SCALE if doing so will not cause
 443 *                 overflow.
 444 *
 445                   IX = ICAMAX( N-1, WORK, 1 )
 446                   XNORM = CABS1( WORK( IX, 1 ) )
 447                   IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
 448      $               GO TO 40
 449                   CALL CSRSCL( N, SCALE, WORK, 1 )
 450                END IF
 451                GO TO 30
 452             END IF
 453 *
 454             SEP( KS ) = ONE / MAX( EST, SMLNUM )
 455          END IF
 456 *
 457    40    CONTINUE
 458          KS = KS + 1
 459    50 CONTINUE
 460       RETURN
 461 *
 462 *     End of CTRSNA
 463 *
 464       END