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

◆ ctrsna()

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

CTRSNA

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

Purpose:
 CTRSNA estimates reciprocal condition numbers for specified
 eigenvalues and/or right eigenvectors of a complex upper triangular
 matrix T (or of any matrix Q*T*Q**H with Q unitary).
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 j-th eigenpair, SELECT(j) 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 COMPLEX array, dimension (LDT,N)
          The upper triangular matrix T.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T. LDT >= max(1,N).
[in]VL
          VL is COMPLEX array, dimension (LDVL,M)
          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
          (or of any Q*T*Q**H with Q unitary), corresponding to the
          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
          must be stored in consecutive columns of VL, as returned by
          CHSEIN or CTREVC.
          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 COMPLEX array, dimension (LDVR,M)
          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
          (or of any Q*T*Q**H with Q unitary), corresponding to the
          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
          must be stored in consecutive columns of VR, as returned by
          CHSEIN or CTREVC.
          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. 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.
          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 COMPLEX 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]RWORK
          RWORK is REAL array, dimension (N)
          If JOB = 'E', RWORK 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**H*u| / (norm(u)*norm(v))

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

249*
250* -- LAPACK computational routine --
251* -- LAPACK is a software package provided by Univ. of Tennessee, --
252* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253*
254* .. Scalar Arguments ..
255 CHARACTER HOWMNY, JOB
256 INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
257* ..
258* .. Array Arguments ..
259 LOGICAL SELECT( * )
260 REAL RWORK( * ), S( * ), SEP( * )
261 COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
262 $ WORK( LDWORK, * )
263* ..
264*
265* =====================================================================
266*
267* .. Parameters ..
268 REAL ZERO, ONE
269 parameter( zero = 0.0e+0, one = 1.0+0 )
270* ..
271* .. Local Scalars ..
272 LOGICAL SOMCON, WANTBH, WANTS, WANTSP
273 CHARACTER NORMIN
274 INTEGER I, IERR, IX, J, K, KASE, KS
275 REAL BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
276 $ XNORM
277 COMPLEX CDUM, PROD
278* ..
279* .. Local Arrays ..
280 INTEGER ISAVE( 3 )
281 COMPLEX DUMMY( 1 )
282* ..
283* .. External Functions ..
284 LOGICAL LSAME
285 INTEGER ICAMAX
286 REAL SCNRM2, SLAMCH
287 COMPLEX CDOTC
288 EXTERNAL lsame, icamax, scnrm2, slamch, cdotc
289* ..
290* .. External Subroutines ..
291 EXTERNAL clacn2, clacpy, clatrs, csrscl, ctrexc, slabad,
292 $ xerbla
293* ..
294* .. Intrinsic Functions ..
295 INTRINSIC abs, aimag, max, real
296* ..
297* .. Statement Functions ..
298 REAL CABS1
299* ..
300* .. Statement Function definitions ..
301 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
302* ..
303* .. Executable Statements ..
304*
305* Decode and test the input parameters
306*
307 wantbh = lsame( job, 'B' )
308 wants = lsame( job, 'E' ) .OR. wantbh
309 wantsp = lsame( job, 'V' ) .OR. wantbh
310*
311 somcon = lsame( howmny, 'S' )
312*
313* Set M to the number of eigenpairs for which condition numbers are
314* to be computed.
315*
316 IF( somcon ) THEN
317 m = 0
318 DO 10 j = 1, n
319 IF( SELECT( j ) )
320 $ m = m + 1
321 10 CONTINUE
322 ELSE
323 m = n
324 END IF
325*
326 info = 0
327 IF( .NOT.wants .AND. .NOT.wantsp ) THEN
328 info = -1
329 ELSE IF( .NOT.lsame( howmny, 'A' ) .AND. .NOT.somcon ) THEN
330 info = -2
331 ELSE IF( n.LT.0 ) THEN
332 info = -4
333 ELSE IF( ldt.LT.max( 1, n ) ) THEN
334 info = -6
335 ELSE IF( ldvl.LT.1 .OR. ( wants .AND. ldvl.LT.n ) ) THEN
336 info = -8
337 ELSE IF( ldvr.LT.1 .OR. ( wants .AND. ldvr.LT.n ) ) THEN
338 info = -10
339 ELSE IF( mm.LT.m ) THEN
340 info = -13
341 ELSE IF( ldwork.LT.1 .OR. ( wantsp .AND. ldwork.LT.n ) ) THEN
342 info = -16
343 END IF
344 IF( info.NE.0 ) THEN
345 CALL xerbla( 'CTRSNA', -info )
346 RETURN
347 END IF
348*
349* Quick return if possible
350*
351 IF( n.EQ.0 )
352 $ RETURN
353*
354 IF( n.EQ.1 ) THEN
355 IF( somcon ) THEN
356 IF( .NOT.SELECT( 1 ) )
357 $ RETURN
358 END IF
359 IF( wants )
360 $ s( 1 ) = one
361 IF( wantsp )
362 $ sep( 1 ) = abs( t( 1, 1 ) )
363 RETURN
364 END IF
365*
366* Get machine constants
367*
368 eps = slamch( 'P' )
369 smlnum = slamch( 'S' ) / eps
370 bignum = one / smlnum
371 CALL slabad( smlnum, bignum )
372*
373 ks = 1
374 DO 50 k = 1, n
375*
376 IF( somcon ) THEN
377 IF( .NOT.SELECT( k ) )
378 $ GO TO 50
379 END IF
380*
381 IF( wants ) THEN
382*
383* Compute the reciprocal condition number of the k-th
384* eigenvalue.
385*
386 prod = cdotc( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
387 rnrm = scnrm2( n, vr( 1, ks ), 1 )
388 lnrm = scnrm2( n, vl( 1, ks ), 1 )
389 s( ks ) = abs( prod ) / ( rnrm*lnrm )
390*
391 END IF
392*
393 IF( wantsp ) THEN
394*
395* Estimate the reciprocal condition number of the k-th
396* eigenvector.
397*
398* Copy the matrix T to the array WORK and swap the k-th
399* diagonal element to the (1,1) position.
400*
401 CALL clacpy( 'Full', n, n, t, ldt, work, ldwork )
402 CALL ctrexc( 'No Q', n, work, ldwork, dummy, 1, k, 1, ierr )
403*
404* Form C = T22 - lambda*I in WORK(2:N,2:N).
405*
406 DO 20 i = 2, n
407 work( i, i ) = work( i, i ) - work( 1, 1 )
408 20 CONTINUE
409*
410* Estimate a lower bound for the 1-norm of inv(C**H). The 1st
411* and (N+1)th columns of WORK are used to store work vectors.
412*
413 sep( ks ) = zero
414 est = zero
415 kase = 0
416 normin = 'N'
417 30 CONTINUE
418 CALL clacn2( n-1, work( 1, n+1 ), work, est, kase, isave )
419*
420 IF( kase.NE.0 ) THEN
421 IF( kase.EQ.1 ) THEN
422*
423* Solve C**H*x = scale*b
424*
425 CALL clatrs( 'Upper', 'Conjugate transpose',
426 $ 'Nonunit', normin, n-1, work( 2, 2 ),
427 $ ldwork, work, scale, rwork, ierr )
428 ELSE
429*
430* Solve C*x = scale*b
431*
432 CALL clatrs( 'Upper', 'No transpose', 'Nonunit',
433 $ normin, n-1, work( 2, 2 ), ldwork, work,
434 $ scale, rwork, ierr )
435 END IF
436 normin = 'Y'
437 IF( scale.NE.one ) THEN
438*
439* Multiply by 1/SCALE if doing so will not cause
440* overflow.
441*
442 ix = icamax( n-1, work, 1 )
443 xnorm = cabs1( work( ix, 1 ) )
444 IF( scale.LT.xnorm*smlnum .OR. scale.EQ.zero )
445 $ GO TO 40
446 CALL csrscl( n, scale, work, 1 )
447 END IF
448 GO TO 30
449 END IF
450*
451 sep( ks ) = one / max( est, smlnum )
452 END IF
453*
454 40 CONTINUE
455 ks = ks + 1
456 50 CONTINUE
457 RETURN
458*
459* End of CTRSNA
460*
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:83
subroutine clatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
CLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition: clatrs.f:239
subroutine csrscl(N, SA, SX, INCX)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
Definition: csrscl.f:84
subroutine clacn2(N, V, X, EST, KASE, ISAVE)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: clacn2.f:133
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine ctrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
CTREXC
Definition: ctrexc.f:126
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition: scnrm2.f90:90
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: