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

◆ 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: