LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
strsna.f
Go to the documentation of this file.
1 *> \brief \b STRSNA
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download STRSNA + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strsna.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strsna.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strsna.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
22 * LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
23 * INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER HOWMNY, JOB
27 * INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
28 * ..
29 * .. Array Arguments ..
30 * LOGICAL SELECT( * )
31 * INTEGER IWORK( * )
32 * REAL S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
33 * $ VR( LDVR, * ), WORK( LDWORK, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> STRSNA estimates reciprocal condition numbers for specified
43 *> eigenvalues and/or right eigenvectors of a real upper
44 *> quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
45 *> orthogonal).
46 *>
47 *> T must be in Schur canonical form (as returned by SHSEQR), that is,
48 *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
49 *> 2-by-2 diagonal block has its diagonal elements equal and its
50 *> off-diagonal elements of opposite sign.
51 *> \endverbatim
52 *
53 * Arguments:
54 * ==========
55 *
56 *> \param[in] JOB
57 *> \verbatim
58 *> JOB is CHARACTER*1
59 *> Specifies whether condition numbers are required for
60 *> eigenvalues (S) or eigenvectors (SEP):
61 *> = 'E': for eigenvalues only (S);
62 *> = 'V': for eigenvectors only (SEP);
63 *> = 'B': for both eigenvalues and eigenvectors (S and SEP).
64 *> \endverbatim
65 *>
66 *> \param[in] HOWMNY
67 *> \verbatim
68 *> HOWMNY is CHARACTER*1
69 *> = 'A': compute condition numbers for all eigenpairs;
70 *> = 'S': compute condition numbers for selected eigenpairs
71 *> specified by the array SELECT.
72 *> \endverbatim
73 *>
74 *> \param[in] SELECT
75 *> \verbatim
76 *> SELECT is LOGICAL array, dimension (N)
77 *> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
78 *> condition numbers are required. To select condition numbers
79 *> for the eigenpair corresponding to a real eigenvalue w(j),
80 *> SELECT(j) must be set to .TRUE.. To select condition numbers
81 *> corresponding to a complex conjugate pair of eigenvalues w(j)
82 *> and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
83 *> set to .TRUE..
84 *> If HOWMNY = 'A', SELECT is not referenced.
85 *> \endverbatim
86 *>
87 *> \param[in] N
88 *> \verbatim
89 *> N is INTEGER
90 *> The order of the matrix T. N >= 0.
91 *> \endverbatim
92 *>
93 *> \param[in] T
94 *> \verbatim
95 *> T is REAL array, dimension (LDT,N)
96 *> The upper quasi-triangular matrix T, in Schur canonical form.
97 *> \endverbatim
98 *>
99 *> \param[in] LDT
100 *> \verbatim
101 *> LDT is INTEGER
102 *> The leading dimension of the array T. LDT >= max(1,N).
103 *> \endverbatim
104 *>
105 *> \param[in] VL
106 *> \verbatim
107 *> VL is REAL array, dimension (LDVL,M)
108 *> If JOB = 'E' or 'B', VL must contain left eigenvectors of T
109 *> (or of any Q*T*Q**T with Q orthogonal), corresponding to the
110 *> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
111 *> must be stored in consecutive columns of VL, as returned by
112 *> SHSEIN or STREVC.
113 *> If JOB = 'V', VL is not referenced.
114 *> \endverbatim
115 *>
116 *> \param[in] LDVL
117 *> \verbatim
118 *> LDVL is INTEGER
119 *> The leading dimension of the array VL.
120 *> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
121 *> \endverbatim
122 *>
123 *> \param[in] VR
124 *> \verbatim
125 *> VR is REAL array, dimension (LDVR,M)
126 *> If JOB = 'E' or 'B', VR must contain right eigenvectors of T
127 *> (or of any Q*T*Q**T with Q orthogonal), corresponding to the
128 *> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
129 *> must be stored in consecutive columns of VR, as returned by
130 *> SHSEIN or STREVC.
131 *> If JOB = 'V', VR is not referenced.
132 *> \endverbatim
133 *>
134 *> \param[in] LDVR
135 *> \verbatim
136 *> LDVR is INTEGER
137 *> The leading dimension of the array VR.
138 *> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
139 *> \endverbatim
140 *>
141 *> \param[out] S
142 *> \verbatim
143 *> S is REAL array, dimension (MM)
144 *> If JOB = 'E' or 'B', the reciprocal condition numbers of the
145 *> selected eigenvalues, stored in consecutive elements of the
146 *> array. For a complex conjugate pair of eigenvalues two
147 *> consecutive elements of S are set to the same value. Thus
148 *> S(j), SEP(j), and the j-th columns of VL and VR all
149 *> correspond to the same eigenpair (but not in general the
150 *> j-th eigenpair, unless all eigenpairs are selected).
151 *> If JOB = 'V', S is not referenced.
152 *> \endverbatim
153 *>
154 *> \param[out] SEP
155 *> \verbatim
156 *> SEP is REAL array, dimension (MM)
157 *> If JOB = 'V' or 'B', the estimated reciprocal condition
158 *> numbers of the selected eigenvectors, stored in consecutive
159 *> elements of the array. For a complex eigenvector two
160 *> consecutive elements of SEP are set to the same value. If
161 *> the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
162 *> is set to 0; this can only occur when the true value would be
163 *> very small anyway.
164 *> If JOB = 'E', SEP is not referenced.
165 *> \endverbatim
166 *>
167 *> \param[in] MM
168 *> \verbatim
169 *> MM is INTEGER
170 *> The number of elements in the arrays S (if JOB = 'E' or 'B')
171 *> and/or SEP (if JOB = 'V' or 'B'). MM >= M.
172 *> \endverbatim
173 *>
174 *> \param[out] M
175 *> \verbatim
176 *> M is INTEGER
177 *> The number of elements of the arrays S and/or SEP actually
178 *> used to store the estimated condition numbers.
179 *> If HOWMNY = 'A', M is set to N.
180 *> \endverbatim
181 *>
182 *> \param[out] WORK
183 *> \verbatim
184 *> WORK is REAL array, dimension (LDWORK,N+6)
185 *> If JOB = 'E', WORK is not referenced.
186 *> \endverbatim
187 *>
188 *> \param[in] LDWORK
189 *> \verbatim
190 *> LDWORK is INTEGER
191 *> The leading dimension of the array WORK.
192 *> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
193 *> \endverbatim
194 *>
195 *> \param[out] IWORK
196 *> \verbatim
197 *> IWORK is INTEGER array, dimension (2*(N-1))
198 *> If JOB = 'E', IWORK is not referenced.
199 *> \endverbatim
200 *>
201 *> \param[out] INFO
202 *> \verbatim
203 *> INFO is INTEGER
204 *> = 0: successful exit
205 *> < 0: if INFO = -i, the i-th argument had an illegal value
206 *> \endverbatim
207 *
208 * Authors:
209 * ========
210 *
211 *> \author Univ. of Tennessee
212 *> \author Univ. of California Berkeley
213 *> \author Univ. of Colorado Denver
214 *> \author NAG Ltd.
215 *
216 *> \ingroup realOTHERcomputational
217 *
218 *> \par Further Details:
219 * =====================
220 *>
221 *> \verbatim
222 *>
223 *> The reciprocal of the condition number of an eigenvalue lambda is
224 *> defined as
225 *>
226 *> S(lambda) = |v**T*u| / (norm(u)*norm(v))
227 *>
228 *> where u and v are the right and left eigenvectors of T corresponding
229 *> to lambda; v**T denotes the transpose of v, and norm(u)
230 *> denotes the Euclidean norm. These reciprocal condition numbers always
231 *> lie between zero (very badly conditioned) and one (very well
232 *> conditioned). If n = 1, S(lambda) is defined to be 1.
233 *>
234 *> An approximate error bound for a computed eigenvalue W(i) is given by
235 *>
236 *> EPS * norm(T) / S(i)
237 *>
238 *> where EPS is the machine precision.
239 *>
240 *> The reciprocal of the condition number of the right eigenvector u
241 *> corresponding to lambda is defined as follows. Suppose
242 *>
243 *> T = ( lambda c )
244 *> ( 0 T22 )
245 *>
246 *> Then the reciprocal condition number is
247 *>
248 *> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
249 *>
250 *> where sigma-min denotes the smallest singular value. We approximate
251 *> the smallest singular value by the reciprocal of an estimate of the
252 *> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
253 *> defined to be abs(T(1,1)).
254 *>
255 *> An approximate error bound for a computed right eigenvector VR(i)
256 *> is given by
257 *>
258 *> EPS * norm(T) / SEP(i)
259 *> \endverbatim
260 *>
261 * =====================================================================
262  SUBROUTINE strsna( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
263  $ LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
264  $ INFO )
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 *
600  END
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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
subroutine strsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, INFO)
STRSNA
Definition: strsna.f:265