LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dtrsna()

subroutine dtrsna ( character  JOB,
character  HOWMNY,
logical, dimension( * )  SELECT,
integer  N,
double precision, dimension( ldt, * )  T,
integer  LDT,
double precision, dimension( ldvl, * )  VL,
integer  LDVL,
double precision, dimension( ldvr, * )  VR,
integer  LDVR,
double precision, dimension( * )  S,
double precision, dimension( * )  SEP,
integer  MM,
integer  M,
double precision, dimension( ldwork, * )  WORK,
integer  LDWORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DTRSNA

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

Purpose:
 DTRSNA estimates reciprocal condition numbers for specified
 eigenvalues and/or right eigenvectors of a real upper
 quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
 orthogonal).

 T must be in Schur canonical form (as returned by DHSEQR), that is,
 block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
 2-by-2 diagonal block has its diagonal elements equal and its
 off-diagonal elements of opposite sign.
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 eigenpair corresponding to a real eigenvalue w(j),
          SELECT(j) must be set to .TRUE.. To select condition numbers
          corresponding to a complex conjugate pair of eigenvalues w(j)
          and w(j+1), either SELECT(j) or SELECT(j+1) or both, 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 DOUBLE PRECISION array, dimension (LDT,N)
          The upper quasi-triangular matrix T, in Schur canonical form.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T. LDT >= max(1,N).
[in]VL
          VL is DOUBLE PRECISION array, dimension (LDVL,M)
          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
          must be stored in consecutive columns of VL, as returned by
          DHSEIN or DTREVC.
          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 DOUBLE PRECISION array, dimension (LDVR,M)
          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
          must be stored in consecutive columns of VR, as returned by
          DHSEIN or DTREVC.
          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. For a complex conjugate pair of eigenvalues two
          consecutive elements of S are set to the same value. 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. For a complex eigenvector two
          consecutive elements of SEP are set to the same value. If
          the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
          is set to 0; this can only occur when the true value would be
          very small anyway.
          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 DOUBLE PRECISION 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]IWORK
          IWORK is INTEGER array, dimension (2*(N-1))
          If JOB = 'E', IWORK 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.
Further Details:
  The reciprocal of the condition number of an eigenvalue lambda is
  defined as

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

  where u and v are the right and left eigenvectors of T corresponding
  to lambda; v**T denotes the 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 262 of file dtrsna.f.

265 *
266 * -- LAPACK computational routine --
267 * -- LAPACK is a software package provided by Univ. of Tennessee, --
268 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
269 *
270 * .. Scalar Arguments ..
271  CHARACTER HOWMNY, JOB
272  INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
273 * ..
274 * .. Array Arguments ..
275  LOGICAL SELECT( * )
276  INTEGER IWORK( * )
277  DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
278  $ VR( LDVR, * ), WORK( LDWORK, * )
279 * ..
280 *
281 * =====================================================================
282 *
283 * .. Parameters ..
284  DOUBLE PRECISION ZERO, ONE, TWO
285  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
286 * ..
287 * .. Local Scalars ..
288  LOGICAL PAIR, SOMCON, WANTBH, WANTS, WANTSP
289  INTEGER I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
290  DOUBLE PRECISION BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
291  $ MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
292 * ..
293 * .. Local Arrays ..
294  INTEGER ISAVE( 3 )
295  DOUBLE PRECISION DUMMY( 1 )
296 * ..
297 * .. External Functions ..
298  LOGICAL LSAME
299  DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
300  EXTERNAL lsame, ddot, dlamch, dlapy2, dnrm2
301 * ..
302 * .. External Subroutines ..
303  EXTERNAL dlabad, dlacn2, dlacpy, dlaqtr, dtrexc, xerbla
304 * ..
305 * .. Intrinsic Functions ..
306  INTRINSIC abs, max, sqrt
307 * ..
308 * .. Executable Statements ..
309 *
310 * Decode and test the input parameters
311 *
312  wantbh = lsame( job, 'B' )
313  wants = lsame( job, 'E' ) .OR. wantbh
314  wantsp = lsame( job, 'V' ) .OR. wantbh
315 *
316  somcon = lsame( howmny, 'S' )
317 *
318  info = 0
319  IF( .NOT.wants .AND. .NOT.wantsp ) THEN
320  info = -1
321  ELSE IF( .NOT.lsame( howmny, 'A' ) .AND. .NOT.somcon ) THEN
322  info = -2
323  ELSE IF( n.LT.0 ) THEN
324  info = -4
325  ELSE IF( ldt.LT.max( 1, n ) ) THEN
326  info = -6
327  ELSE IF( ldvl.LT.1 .OR. ( wants .AND. ldvl.LT.n ) ) THEN
328  info = -8
329  ELSE IF( ldvr.LT.1 .OR. ( wants .AND. ldvr.LT.n ) ) THEN
330  info = -10
331  ELSE
332 *
333 * Set M to the number of eigenpairs for which condition numbers
334 * are required, and test MM.
335 *
336  IF( somcon ) THEN
337  m = 0
338  pair = .false.
339  DO 10 k = 1, n
340  IF( pair ) THEN
341  pair = .false.
342  ELSE
343  IF( k.LT.n ) THEN
344  IF( t( k+1, k ).EQ.zero ) THEN
345  IF( SELECT( k ) )
346  $ m = m + 1
347  ELSE
348  pair = .true.
349  IF( SELECT( k ) .OR. SELECT( k+1 ) )
350  $ m = m + 2
351  END IF
352  ELSE
353  IF( SELECT( n ) )
354  $ m = m + 1
355  END IF
356  END IF
357  10 CONTINUE
358  ELSE
359  m = n
360  END IF
361 *
362  IF( mm.LT.m ) THEN
363  info = -13
364  ELSE IF( ldwork.LT.1 .OR. ( wantsp .AND. ldwork.LT.n ) ) THEN
365  info = -16
366  END IF
367  END IF
368  IF( info.NE.0 ) THEN
369  CALL xerbla( 'DTRSNA', -info )
370  RETURN
371  END IF
372 *
373 * Quick return if possible
374 *
375  IF( n.EQ.0 )
376  $ RETURN
377 *
378  IF( n.EQ.1 ) THEN
379  IF( somcon ) THEN
380  IF( .NOT.SELECT( 1 ) )
381  $ RETURN
382  END IF
383  IF( wants )
384  $ s( 1 ) = one
385  IF( wantsp )
386  $ sep( 1 ) = abs( t( 1, 1 ) )
387  RETURN
388  END IF
389 *
390 * Get machine constants
391 *
392  eps = dlamch( 'P' )
393  smlnum = dlamch( 'S' ) / eps
394  bignum = one / smlnum
395  CALL dlabad( smlnum, bignum )
396 *
397  ks = 0
398  pair = .false.
399  DO 60 k = 1, n
400 *
401 * Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
402 *
403  IF( pair ) THEN
404  pair = .false.
405  GO TO 60
406  ELSE
407  IF( k.LT.n )
408  $ pair = t( k+1, k ).NE.zero
409  END IF
410 *
411 * Determine whether condition numbers are required for the k-th
412 * eigenpair.
413 *
414  IF( somcon ) THEN
415  IF( pair ) THEN
416  IF( .NOT.SELECT( k ) .AND. .NOT.SELECT( k+1 ) )
417  $ GO TO 60
418  ELSE
419  IF( .NOT.SELECT( k ) )
420  $ GO TO 60
421  END IF
422  END IF
423 *
424  ks = ks + 1
425 *
426  IF( wants ) THEN
427 *
428 * Compute the reciprocal condition number of the k-th
429 * eigenvalue.
430 *
431  IF( .NOT.pair ) THEN
432 *
433 * Real eigenvalue.
434 *
435  prod = ddot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
436  rnrm = dnrm2( n, vr( 1, ks ), 1 )
437  lnrm = dnrm2( n, vl( 1, ks ), 1 )
438  s( ks ) = abs( prod ) / ( rnrm*lnrm )
439  ELSE
440 *
441 * Complex eigenvalue.
442 *
443  prod1 = ddot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
444  prod1 = prod1 + ddot( n, vr( 1, ks+1 ), 1, vl( 1, ks+1 ),
445  $ 1 )
446  prod2 = ddot( n, vl( 1, ks ), 1, vr( 1, ks+1 ), 1 )
447  prod2 = prod2 - ddot( n, vl( 1, ks+1 ), 1, vr( 1, ks ),
448  $ 1 )
449  rnrm = dlapy2( dnrm2( n, vr( 1, ks ), 1 ),
450  $ dnrm2( n, vr( 1, ks+1 ), 1 ) )
451  lnrm = dlapy2( dnrm2( n, vl( 1, ks ), 1 ),
452  $ dnrm2( n, vl( 1, ks+1 ), 1 ) )
453  cond = dlapy2( prod1, prod2 ) / ( rnrm*lnrm )
454  s( ks ) = cond
455  s( ks+1 ) = cond
456  END IF
457  END IF
458 *
459  IF( wantsp ) THEN
460 *
461 * Estimate the reciprocal condition number of the k-th
462 * eigenvector.
463 *
464 * Copy the matrix T to the array WORK and swap the diagonal
465 * block beginning at T(k,k) to the (1,1) position.
466 *
467  CALL dlacpy( 'Full', n, n, t, ldt, work, ldwork )
468  ifst = k
469  ilst = 1
470  CALL dtrexc( 'No Q', n, work, ldwork, dummy, 1, ifst, ilst,
471  $ work( 1, n+1 ), ierr )
472 *
473  IF( ierr.EQ.1 .OR. ierr.EQ.2 ) THEN
474 *
475 * Could not swap because blocks not well separated
476 *
477  scale = one
478  est = bignum
479  ELSE
480 *
481 * Reordering successful
482 *
483  IF( work( 2, 1 ).EQ.zero ) THEN
484 *
485 * Form C = T22 - lambda*I in WORK(2:N,2:N).
486 *
487  DO 20 i = 2, n
488  work( i, i ) = work( i, i ) - work( 1, 1 )
489  20 CONTINUE
490  n2 = 1
491  nn = n - 1
492  ELSE
493 *
494 * Triangularize the 2 by 2 block by unitary
495 * transformation U = [ cs i*ss ]
496 * [ i*ss cs ].
497 * such that the (1,1) position of WORK is complex
498 * eigenvalue lambda with positive imaginary part. (2,2)
499 * position of WORK is the complex eigenvalue lambda
500 * with negative imaginary part.
501 *
502  mu = sqrt( abs( work( 1, 2 ) ) )*
503  $ sqrt( abs( work( 2, 1 ) ) )
504  delta = dlapy2( mu, work( 2, 1 ) )
505  cs = mu / delta
506  sn = -work( 2, 1 ) / delta
507 *
508 * Form
509 *
510 * C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
511 * [ mu ]
512 * [ .. ]
513 * [ .. ]
514 * [ mu ]
515 * where C**T is transpose of matrix C,
516 * and RWORK is stored starting in the N+1-st column of
517 * WORK.
518 *
519  DO 30 j = 3, n
520  work( 2, j ) = cs*work( 2, j )
521  work( j, j ) = work( j, j ) - work( 1, 1 )
522  30 CONTINUE
523  work( 2, 2 ) = zero
524 *
525  work( 1, n+1 ) = two*mu
526  DO 40 i = 2, n - 1
527  work( i, n+1 ) = sn*work( 1, i+1 )
528  40 CONTINUE
529  n2 = 2
530  nn = 2*( n-1 )
531  END IF
532 *
533 * Estimate norm(inv(C**T))
534 *
535  est = zero
536  kase = 0
537  50 CONTINUE
538  CALL dlacn2( nn, work( 1, n+2 ), work( 1, n+4 ), iwork,
539  $ est, kase, isave )
540  IF( kase.NE.0 ) THEN
541  IF( kase.EQ.1 ) THEN
542  IF( n2.EQ.1 ) THEN
543 *
544 * Real eigenvalue: solve C**T*x = scale*c.
545 *
546  CALL dlaqtr( .true., .true., n-1, work( 2, 2 ),
547  $ ldwork, dummy, dumm, scale,
548  $ work( 1, n+4 ), work( 1, n+6 ),
549  $ ierr )
550  ELSE
551 *
552 * Complex eigenvalue: solve
553 * C**T*(p+iq) = scale*(c+id) in real arithmetic.
554 *
555  CALL dlaqtr( .true., .false., n-1, work( 2, 2 ),
556  $ ldwork, work( 1, n+1 ), mu, scale,
557  $ work( 1, n+4 ), work( 1, n+6 ),
558  $ ierr )
559  END IF
560  ELSE
561  IF( n2.EQ.1 ) THEN
562 *
563 * Real eigenvalue: solve C*x = scale*c.
564 *
565  CALL dlaqtr( .false., .true., n-1, work( 2, 2 ),
566  $ ldwork, dummy, dumm, scale,
567  $ work( 1, n+4 ), work( 1, n+6 ),
568  $ ierr )
569  ELSE
570 *
571 * Complex eigenvalue: solve
572 * C*(p+iq) = scale*(c+id) in real arithmetic.
573 *
574  CALL dlaqtr( .false., .false., n-1,
575  $ work( 2, 2 ), ldwork,
576  $ work( 1, n+1 ), mu, scale,
577  $ work( 1, n+4 ), work( 1, n+6 ),
578  $ ierr )
579 *
580  END IF
581  END IF
582 *
583  GO TO 50
584  END IF
585  END IF
586 *
587  sep( ks ) = scale / max( est, smlnum )
588  IF( pair )
589  $ sep( ks+1 ) = sep( ks )
590  END IF
591 *
592  IF( pair )
593  $ ks = ks + 1
594 *
595  60 CONTINUE
596  RETURN
597 *
598 * End of DTRSNA
599 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:63
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:82
subroutine dlaqtr(LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO)
DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of sp...
Definition: dlaqtr.f:165
subroutine dlacn2(N, V, X, ISGN, EST, KASE, ISAVE)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: dlacn2.f:136
subroutine dtrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
DTREXC
Definition: dtrexc.f:148
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition: dnrm2.f90:89
Here is the call graph for this function:
Here is the caller graph for this function: