LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ztrsna ( character  JOB,
character  HOWMNY,
logical, dimension( * )  SELECT,
integer  N,
complex*16, dimension( ldt, * )  T,
integer  LDT,
complex*16, dimension( ldvl, * )  VL,
integer  LDVL,
complex*16, dimension( ldvr, * )  VR,
integer  LDVR,
double precision, dimension( * )  S,
double precision, dimension( * )  SEP,
integer  MM,
integer  M,
complex*16, dimension( ldwork, * )  WORK,
integer  LDWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZTRSNA

Download ZTRSNA + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZTRSNA estimates reciprocal condition numbers for specified
 eigenvalues and/or right eigenvectors of a complex upper triangular
 matrix T (or of any matrix Q*T*Q**H with Q unitary).
Parameters
[in]JOB
          JOB is CHARACTER*1
          Specifies whether condition numbers are required for
          eigenvalues (S) or eigenvectors (SEP):
          = 'E': for eigenvalues only (S);
          = 'V': for eigenvectors only (SEP);
          = 'B': for both eigenvalues and eigenvectors (S and SEP).
[in]HOWMNY
          HOWMNY is CHARACTER*1
          = 'A': compute condition numbers for all eigenpairs;
          = 'S': compute condition numbers for selected eigenpairs
                 specified by the array SELECT.
[in]SELECT
          SELECT is LOGICAL array, dimension (N)
          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
          condition numbers are required. To select condition numbers
          for the j-th eigenpair, SELECT(j) must be set to .TRUE..
          If HOWMNY = 'A', SELECT is not referenced.
[in]N
          N is INTEGER
          The order of the matrix T. N >= 0.
[in]T
          T is COMPLEX*16 array, dimension (LDT,N)
          The upper triangular matrix T.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T. LDT >= max(1,N).
[in]VL
          VL is COMPLEX*16 array, dimension (LDVL,M)
          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
          (or of any Q*T*Q**H with Q unitary), corresponding to the
          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
          must be stored in consecutive columns of VL, as returned by
          ZHSEIN or ZTREVC.
          If JOB = 'V', VL is not referenced.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the array VL.
          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
[in]VR
          VR is COMPLEX*16 array, dimension (LDVR,M)
          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
          (or of any Q*T*Q**H with Q unitary), corresponding to the
          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
          must be stored in consecutive columns of VR, as returned by
          ZHSEIN or ZTREVC.
          If JOB = 'V', VR is not referenced.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the array VR.
          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
[out]S
          S is DOUBLE PRECISION array, dimension (MM)
          If JOB = 'E' or 'B', the reciprocal condition numbers of the
          selected eigenvalues, stored in consecutive elements of the
          array. Thus S(j), SEP(j), and the j-th columns of VL and VR
          all correspond to the same eigenpair (but not in general the
          j-th eigenpair, unless all eigenpairs are selected).
          If JOB = 'V', S is not referenced.
[out]SEP
          SEP is DOUBLE PRECISION array, dimension (MM)
          If JOB = 'V' or 'B', the estimated reciprocal condition
          numbers of the selected eigenvectors, stored in consecutive
          elements of the array.
          If JOB = 'E', SEP is not referenced.
[in]MM
          MM is INTEGER
          The number of elements in the arrays S (if JOB = 'E' or 'B')
           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
[out]M
          M is INTEGER
          The number of elements of the arrays S and/or SEP actually
          used to store the estimated condition numbers.
          If HOWMNY = 'A', M is set to N.
[out]WORK
          WORK is COMPLEX*16 array, dimension (LDWORK,N+6)
          If JOB = 'E', WORK is not referenced.
[in]LDWORK
          LDWORK is INTEGER
          The leading dimension of the array WORK.
          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
          If JOB = 'E', RWORK is not referenced.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
Further Details:
  The reciprocal of the condition number of an eigenvalue lambda is
  defined as

          S(lambda) = |v**H*u| / (norm(u)*norm(v))

  where u and v are the right and left eigenvectors of T corresponding
  to lambda; v**H denotes the conjugate transpose of v, and norm(u)
  denotes the Euclidean norm. These reciprocal condition numbers always
  lie between zero (very badly conditioned) and one (very well
  conditioned). If n = 1, S(lambda) is defined to be 1.

  An approximate error bound for a computed eigenvalue W(i) is given by

                      EPS * norm(T) / S(i)

  where EPS is the machine precision.

  The reciprocal of the condition number of the right eigenvector u
  corresponding to lambda is defined as follows. Suppose

              T = ( lambda  c  )
                  (   0    T22 )

  Then the reciprocal condition number is

          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )

  where sigma-min denotes the smallest singular value. We approximate
  the smallest singular value by the reciprocal of an estimate of the
  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
  defined to be abs(T(1,1)).

  An approximate error bound for a computed right eigenvector VR(i)
  is given by

                      EPS * norm(T) / SEP(i)

Definition at line 251 of file ztrsna.f.

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 *
subroutine zdrscl(N, SA, SX, INCX)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: zdrscl.f:86
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
Definition: zdotc.f:54
subroutine zlacn2(N, V, X, EST, KASE, ISAVE)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: zlacn2.f:135
double precision function dznrm2(N, X, INCX)
DZNRM2
Definition: dznrm2.f:56
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
subroutine ztrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
ZTREXC
Definition: ztrexc.f:126
subroutine zlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
ZLATRS solves a triangular system of equations with the scale factor set to prevent overflow...
Definition: zlatrs.f:241
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: