LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ ztrsna()

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 2017
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.8.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 2017
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  $ dlabad
296 * ..
297 * .. Intrinsic Functions ..
298  INTRINSIC abs, dble, dimag, max
299 * ..
300 * .. Statement Functions ..
301  DOUBLE PRECISION cabs1
302 * ..
303 * .. Statement Function definitions ..
304  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( 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( 'ZTRSNA', -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 = dlamch( 'P' )
372  smlnum = dlamch( 'S' ) / eps
373  bignum = one / smlnum
374  CALL dlabad( 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 = zdotc( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
390  rnrm = dznrm2( n, vr( 1, ks ), 1 )
391  lnrm = dznrm2( 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 zlacpy( 'Full', n, n, t, ldt, work, ldwork )
405  CALL ztrexc( '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 zlacn2( 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 zlatrs( '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 zlatrs( '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 = izamax( 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 zdrscl( 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 ZTRSNA
463 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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
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
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
Definition: zdotc.f:85
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:73
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:77
subroutine zdrscl(N, SA, SX, INCX)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: zdrscl.f:86
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine ztrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
ZTREXC
Definition: ztrexc.f:128
Here is the call graph for this function:
Here is the caller graph for this function: