LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 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*
396 ks = 0
397 pair = .false.
398 DO 60 k = 1, n
399*
400* Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
401*
402 IF( pair ) THEN
403 pair = .false.
404 GO TO 60
405 ELSE
406 IF( k.LT.n )
407 $ pair = t( k+1, k ).NE.zero
408 END IF
409*
410* Determine whether condition numbers are required for the k-th
411* eigenpair.
412*
413 IF( somcon ) THEN
414 IF( pair ) THEN
415 IF( .NOT.SELECT( k ) .AND. .NOT.SELECT( k+1 ) )
416 $ GO TO 60
417 ELSE
418 IF( .NOT.SELECT( k ) )
419 $ GO TO 60
420 END IF
421 END IF
422*
423 ks = ks + 1
424*
425 IF( wants ) THEN
426*
427* Compute the reciprocal condition number of the k-th
428* eigenvalue.
429*
430 IF( .NOT.pair ) THEN
431*
432* Real eigenvalue.
433*
434 prod = ddot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
435 rnrm = dnrm2( n, vr( 1, ks ), 1 )
436 lnrm = dnrm2( n, vl( 1, ks ), 1 )
437 s( ks ) = abs( prod ) / ( rnrm*lnrm )
438 ELSE
439*
440* Complex eigenvalue.
441*
442 prod1 = ddot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
443 prod1 = prod1 + ddot( n, vr( 1, ks+1 ), 1, vl( 1, ks+1 ),
444 $ 1 )
445 prod2 = ddot( n, vl( 1, ks ), 1, vr( 1, ks+1 ), 1 )
446 prod2 = prod2 - ddot( n, vl( 1, ks+1 ), 1, vr( 1, ks ),
447 $ 1 )
448 rnrm = dlapy2( dnrm2( n, vr( 1, ks ), 1 ),
449 $ dnrm2( n, vr( 1, ks+1 ), 1 ) )
450 lnrm = dlapy2( dnrm2( n, vl( 1, ks ), 1 ),
451 $ dnrm2( n, vl( 1, ks+1 ), 1 ) )
452 cond = dlapy2( prod1, prod2 ) / ( rnrm*lnrm )
453 s( ks ) = cond
454 s( ks+1 ) = cond
455 END IF
456 END IF
457*
458 IF( wantsp ) THEN
459*
460* Estimate the reciprocal condition number of the k-th
461* eigenvector.
462*
463* Copy the matrix T to the array WORK and swap the diagonal
464* block beginning at T(k,k) to the (1,1) position.
465*
466 CALL dlacpy( 'Full', n, n, t, ldt, work, ldwork )
467 ifst = k
468 ilst = 1
469 CALL dtrexc( 'No Q', n, work, ldwork, dummy, 1, ifst, ilst,
470 $ work( 1, n+1 ), ierr )
471*
472 IF( ierr.EQ.1 .OR. ierr.EQ.2 ) THEN
473*
474* Could not swap because blocks not well separated
475*
476 scale = one
477 est = bignum
478 ELSE
479*
480* Reordering successful
481*
482 IF( work( 2, 1 ).EQ.zero ) THEN
483*
484* Form C = T22 - lambda*I in WORK(2:N,2:N).
485*
486 DO 20 i = 2, n
487 work( i, i ) = work( i, i ) - work( 1, 1 )
488 20 CONTINUE
489 n2 = 1
490 nn = n - 1
491 ELSE
492*
493* Triangularize the 2 by 2 block by unitary
494* transformation U = [ cs i*ss ]
495* [ i*ss cs ].
496* such that the (1,1) position of WORK is complex
497* eigenvalue lambda with positive imaginary part. (2,2)
498* position of WORK is the complex eigenvalue lambda
499* with negative imaginary part.
500*
501 mu = sqrt( abs( work( 1, 2 ) ) )*
502 $ sqrt( abs( work( 2, 1 ) ) )
503 delta = dlapy2( mu, work( 2, 1 ) )
504 cs = mu / delta
505 sn = -work( 2, 1 ) / delta
506*
507* Form
508*
509* C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
510* [ mu ]
511* [ .. ]
512* [ .. ]
513* [ mu ]
514* where C**T is transpose of matrix C,
515* and RWORK is stored starting in the N+1-st column of
516* WORK.
517*
518 DO 30 j = 3, n
519 work( 2, j ) = cs*work( 2, j )
520 work( j, j ) = work( j, j ) - work( 1, 1 )
521 30 CONTINUE
522 work( 2, 2 ) = zero
523*
524 work( 1, n+1 ) = two*mu
525 DO 40 i = 2, n - 1
526 work( i, n+1 ) = sn*work( 1, i+1 )
527 40 CONTINUE
528 n2 = 2
529 nn = 2*( n-1 )
530 END IF
531*
532* Estimate norm(inv(C**T))
533*
534 est = zero
535 kase = 0
536 50 CONTINUE
537 CALL dlacn2( nn, work( 1, n+2 ), work( 1, n+4 ), iwork,
538 $ est, kase, isave )
539 IF( kase.NE.0 ) THEN
540 IF( kase.EQ.1 ) THEN
541 IF( n2.EQ.1 ) THEN
542*
543* Real eigenvalue: solve C**T*x = scale*c.
544*
545 CALL dlaqtr( .true., .true., n-1, work( 2, 2 ),
546 $ ldwork, dummy, dumm, scale,
547 $ work( 1, n+4 ), work( 1, n+6 ),
548 $ ierr )
549 ELSE
550*
551* Complex eigenvalue: solve
552* C**T*(p+iq) = scale*(c+id) in real arithmetic.
553*
554 CALL dlaqtr( .true., .false., n-1, work( 2, 2 ),
555 $ ldwork, work( 1, n+1 ), mu, scale,
556 $ work( 1, n+4 ), work( 1, n+6 ),
557 $ ierr )
558 END IF
559 ELSE
560 IF( n2.EQ.1 ) THEN
561*
562* Real eigenvalue: solve C*x = scale*c.
563*
564 CALL dlaqtr( .false., .true., n-1, work( 2, 2 ),
565 $ ldwork, dummy, dumm, scale,
566 $ work( 1, n+4 ), work( 1, n+6 ),
567 $ ierr )
568 ELSE
569*
570* Complex eigenvalue: solve
571* C*(p+iq) = scale*(c+id) in real arithmetic.
572*
573 CALL dlaqtr( .false., .false., n-1,
574 $ work( 2, 2 ), ldwork,
575 $ work( 1, n+1 ), mu, scale,
576 $ work( 1, n+4 ), work( 1, n+6 ),
577 $ ierr )
578*
579 END IF
580 END IF
581*
582 GO TO 50
583 END IF
584 END IF
585*
586 sep( ks ) = scale / max( est, smlnum )
587 IF( pair )
588 $ sep( ks+1 ) = sep( ks )
589 END IF
590*
591 IF( pair )
592 $ ks = ks + 1
593*
594 60 CONTINUE
595 RETURN
596*
597* End of DTRSNA
598*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
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 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 dlamch(cmach)
DLAMCH
Definition dlamch.f:69
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:63
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
subroutine dtrexc(compq, n, t, ldt, q, ldq, ifst, ilst, work, info)
DTREXC
Definition dtrexc.f:148
Here is the call graph for this function:
Here is the caller graph for this function: