LAPACK  3.10.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.
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 strsna.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  REAL S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
278  $ VR( LDVR, * ), WORK( LDWORK, * )
279 * ..
280 *
281 * =====================================================================
282 *
283 * .. Parameters ..
284  REAL ZERO, ONE, TWO
285  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+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  REAL 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  REAL DUMMY( 1 )
296 * ..
297 * .. External Functions ..
298  LOGICAL LSAME
299  REAL SDOT, SLAMCH, SLAPY2, SNRM2
300  EXTERNAL lsame, sdot, slamch, slapy2, snrm2
301 * ..
302 * .. External Subroutines ..
303  EXTERNAL slabad, slacn2, slacpy, slaqtr, strexc, 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( 'STRSNA', -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 = slamch( 'P' )
393  smlnum = slamch( 'S' ) / eps
394  bignum = one / smlnum
395  CALL slabad( 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 = sdot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
436  rnrm = snrm2( n, vr( 1, ks ), 1 )
437  lnrm = snrm2( n, vl( 1, ks ), 1 )
438  s( ks ) = abs( prod ) / ( rnrm*lnrm )
439  ELSE
440 *
441 * Complex eigenvalue.
442 *
443  prod1 = sdot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
444  prod1 = prod1 + sdot( n, vr( 1, ks+1 ), 1, vl( 1, ks+1 ),
445  $ 1 )
446  prod2 = sdot( n, vl( 1, ks ), 1, vr( 1, ks+1 ), 1 )
447  prod2 = prod2 - sdot( n, vl( 1, ks+1 ), 1, vr( 1, ks ),
448  $ 1 )
449  rnrm = slapy2( snrm2( n, vr( 1, ks ), 1 ),
450  $ snrm2( n, vr( 1, ks+1 ), 1 ) )
451  lnrm = slapy2( snrm2( n, vl( 1, ks ), 1 ),
452  $ snrm2( n, vl( 1, ks+1 ), 1 ) )
453  cond = slapy2( 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 slacpy( 'Full', n, n, t, ldt, work, ldwork )
468  ifst = k
469  ilst = 1
470  CALL strexc( '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 = slapy2( 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 slacn2( 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 slaqtr( .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 slaqtr( .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 slaqtr( .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 slaqtr( .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 STRSNA
599 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:63
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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:136
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:165
subroutine strexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
STREXC
Definition: strexc.f:148
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
real(wp) function snrm2(n, x, incx)
SNRM2
Definition: snrm2.f90:89
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: