LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
dtrsna.f
Go to the documentation of this file.
1*> \brief \b DTRSNA
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsna.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsna.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsna.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DTRSNA( 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* DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
33* \$ VR( LDVR, * ), WORK( LDWORK, * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> DTRSNA 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 DHSEQR), 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*> DHSEIN or DTREVC.
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 DOUBLE PRECISION 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*> DHSEIN or DTREVC.
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleOTHERcomputational
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 dtrsna( 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 DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
278 \$ vr( ldvr, * ), work( ldwork, * )
279* ..
280*
281* =====================================================================
282*
283* .. Parameters ..
284 DOUBLE PRECISION ZERO, ONE, TWO
285 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0, two = 2.0d+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 DOUBLE PRECISION 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 DOUBLE PRECISION DUMMY( 1 )
296* ..
297* .. External Functions ..
298 LOGICAL LSAME
299 DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
300 EXTERNAL lsame, ddot, dlamch, dlapy2, dnrm2
301* ..
302* .. External Subroutines ..
303 EXTERNAL dlabad, dlacn2, dlacpy, dlaqtr, dtrexc, 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( 'DTRSNA', -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 = dlamch( 'P' )
393 smlnum = dlamch( 'S' ) / eps
394 bignum = one / smlnum
395 CALL dlabad( 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 = ddot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
436 rnrm = dnrm2( n, vr( 1, ks ), 1 )
437 lnrm = dnrm2( n, vl( 1, ks ), 1 )
438 s( ks ) = abs( prod ) / ( rnrm*lnrm )
439 ELSE
440*
441* Complex eigenvalue.
442*
443 prod1 = ddot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
444 prod1 = prod1 + ddot( n, vr( 1, ks+1 ), 1, vl( 1, ks+1 ),
445 \$ 1 )
446 prod2 = ddot( n, vl( 1, ks ), 1, vr( 1, ks+1 ), 1 )
447 prod2 = prod2 - ddot( n, vl( 1, ks+1 ), 1, vr( 1, ks ),
448 \$ 1 )
449 rnrm = dlapy2( dnrm2( n, vr( 1, ks ), 1 ),
450 \$ dnrm2( n, vr( 1, ks+1 ), 1 ) )
451 lnrm = dlapy2( dnrm2( n, vl( 1, ks ), 1 ),
452 \$ dnrm2( n, vl( 1, ks+1 ), 1 ) )
453 cond = dlapy2( 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 dlacpy( 'Full', n, n, t, ldt, work, ldwork )
468 ifst = k
469 ilst = 1
470 CALL dtrexc( '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 = dlapy2( 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 dlacn2( 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 dlaqtr( .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 dlaqtr( .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 dlaqtr( .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 dlaqtr( .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 DTRSNA
599*
600 END
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlaqtr(LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO)
DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of sp...
Definition: dlaqtr.f:165
subroutine dlacn2(N, V, X, ISGN, EST, KASE, ISAVE)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: dlacn2.f:136
subroutine dtrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
DTREXC
Definition: dtrexc.f:148
subroutine dtrsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, INFO)
DTRSNA
Definition: dtrsna.f:265