LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ strsna()

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

STRSNA

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

Purpose:
 STRSNA 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 SHSEQR), 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 REAL 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 REAL 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
          SHSEIN or STREVC.
          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 REAL 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
          SHSEIN or STREVC.
          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 REAL 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 REAL 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 REAL 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.
Date
December 2016
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 267 of file strsna.f.

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