LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 trsna
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 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*
396 ks = 0
397 pair = .false.
398 DO 60 k = 1, n
399*
400* Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
401*
402 IF( pair ) THEN
403 pair = .false.
404 GO TO 60
405 ELSE
406 IF( k.LT.n )
407 $ pair = t( k+1, k ).NE.zero
408 END IF
409*
410* Determine whether condition numbers are required for the k-th
411* eigenpair.
412*
413 IF( somcon ) THEN
414 IF( pair ) THEN
415 IF( .NOT.SELECT( k ) .AND. .NOT.SELECT( k+1 ) )
416 $ GO TO 60
417 ELSE
418 IF( .NOT.SELECT( k ) )
419 $ GO TO 60
420 END IF
421 END IF
422*
423 ks = ks + 1
424*
425 IF( wants ) THEN
426*
427* Compute the reciprocal condition number of the k-th
428* eigenvalue.
429*
430 IF( .NOT.pair ) THEN
431*
432* Real eigenvalue.
433*
434 prod = sdot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
435 rnrm = snrm2( n, vr( 1, ks ), 1 )
436 lnrm = snrm2( n, vl( 1, ks ), 1 )
437 s( ks ) = abs( prod ) / ( rnrm*lnrm )
438 ELSE
439*
440* Complex eigenvalue.
441*
442 prod1 = sdot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
443 prod1 = prod1 + sdot( n, vr( 1, ks+1 ), 1, vl( 1, ks+1 ),
444 $ 1 )
445 prod2 = sdot( n, vl( 1, ks ), 1, vr( 1, ks+1 ), 1 )
446 prod2 = prod2 - sdot( n, vl( 1, ks+1 ), 1, vr( 1, ks ),
447 $ 1 )
448 rnrm = slapy2( snrm2( n, vr( 1, ks ), 1 ),
449 $ snrm2( n, vr( 1, ks+1 ), 1 ) )
450 lnrm = slapy2( snrm2( n, vl( 1, ks ), 1 ),
451 $ snrm2( n, vl( 1, ks+1 ), 1 ) )
452 cond = slapy2( prod1, prod2 ) / ( rnrm*lnrm )
453 s( ks ) = cond
454 s( ks+1 ) = cond
455 END IF
456 END IF
457*
458 IF( wantsp ) THEN
459*
460* Estimate the reciprocal condition number of the k-th
461* eigenvector.
462*
463* Copy the matrix T to the array WORK and swap the diagonal
464* block beginning at T(k,k) to the (1,1) position.
465*
466 CALL slacpy( 'Full', n, n, t, ldt, work, ldwork )
467 ifst = k
468 ilst = 1
469 CALL strexc( 'No Q', n, work, ldwork, dummy, 1, ifst, ilst,
470 $ work( 1, n+1 ), ierr )
471*
472 IF( ierr.EQ.1 .OR. ierr.EQ.2 ) THEN
473*
474* Could not swap because blocks not well separated
475*
476 scale = one
477 est = bignum
478 ELSE
479*
480* Reordering successful
481*
482 IF( work( 2, 1 ).EQ.zero ) THEN
483*
484* Form C = T22 - lambda*I in WORK(2:N,2:N).
485*
486 DO 20 i = 2, n
487 work( i, i ) = work( i, i ) - work( 1, 1 )
488 20 CONTINUE
489 n2 = 1
490 nn = n - 1
491 ELSE
492*
493* Triangularize the 2 by 2 block by unitary
494* transformation U = [ cs i*ss ]
495* [ i*ss cs ].
496* such that the (1,1) position of WORK is complex
497* eigenvalue lambda with positive imaginary part. (2,2)
498* position of WORK is the complex eigenvalue lambda
499* with negative imaginary part.
500*
501 mu = sqrt( abs( work( 1, 2 ) ) )*
502 $ sqrt( abs( work( 2, 1 ) ) )
503 delta = slapy2( mu, work( 2, 1 ) )
504 cs = mu / delta
505 sn = -work( 2, 1 ) / delta
506*
507* Form
508*
509* C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
510* [ mu ]
511* [ .. ]
512* [ .. ]
513* [ mu ]
514* where C**T is transpose of matrix C,
515* and RWORK is stored starting in the N+1-st column of
516* WORK.
517*
518 DO 30 j = 3, n
519 work( 2, j ) = cs*work( 2, j )
520 work( j, j ) = work( j, j ) - work( 1, 1 )
521 30 CONTINUE
522 work( 2, 2 ) = zero
523*
524 work( 1, n+1 ) = two*mu
525 DO 40 i = 2, n - 1
526 work( i, n+1 ) = sn*work( 1, i+1 )
527 40 CONTINUE
528 n2 = 2
529 nn = 2*( n-1 )
530 END IF
531*
532* Estimate norm(inv(C**T))
533*
534 est = zero
535 kase = 0
536 50 CONTINUE
537 CALL slacn2( nn, work( 1, n+2 ), work( 1, n+4 ), iwork,
538 $ est, kase, isave )
539 IF( kase.NE.0 ) THEN
540 IF( kase.EQ.1 ) THEN
541 IF( n2.EQ.1 ) THEN
542*
543* Real eigenvalue: solve C**T*x = scale*c.
544*
545 CALL slaqtr( .true., .true., n-1, work( 2, 2 ),
546 $ ldwork, dummy, dumm, scale,
547 $ work( 1, n+4 ), work( 1, n+6 ),
548 $ ierr )
549 ELSE
550*
551* Complex eigenvalue: solve
552* C**T*(p+iq) = scale*(c+id) in real arithmetic.
553*
554 CALL slaqtr( .true., .false., n-1, work( 2, 2 ),
555 $ ldwork, work( 1, n+1 ), mu, scale,
556 $ work( 1, n+4 ), work( 1, n+6 ),
557 $ ierr )
558 END IF
559 ELSE
560 IF( n2.EQ.1 ) THEN
561*
562* Real eigenvalue: solve C*x = scale*c.
563*
564 CALL slaqtr( .false., .true., n-1, work( 2, 2 ),
565 $ ldwork, dummy, dumm, scale,
566 $ work( 1, n+4 ), work( 1, n+6 ),
567 $ ierr )
568 ELSE
569*
570* Complex eigenvalue: solve
571* C*(p+iq) = scale*(c+id) in real arithmetic.
572*
573 CALL slaqtr( .false., .false., n-1,
574 $ work( 2, 2 ), ldwork,
575 $ work( 1, n+1 ), mu, scale,
576 $ work( 1, n+4 ), work( 1, n+6 ),
577 $ ierr )
578*
579 END IF
580 END IF
581*
582 GO TO 50
583 END IF
584 END IF
585*
586 sep( ks ) = scale / max( est, smlnum )
587 IF( pair )
588 $ sep( ks+1 ) = sep( ks )
589 END IF
590*
591 IF( pair )
592 $ ks = ks + 1
593*
594 60 CONTINUE
595 RETURN
596*
597* End of STRSNA
598*
599 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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 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 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