LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

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