LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ztrsna.f
Go to the documentation of this file.
1*> \brief \b ZTRSNA
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZTRSNA + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsna.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsna.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsna.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
22* LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
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* DOUBLE PRECISION RWORK( * ), S( * ), SEP( * )
32* COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
33* $ WORK( LDWORK, * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> ZTRSNA estimates reciprocal condition numbers for specified
43*> eigenvalues and/or right eigenvectors of a complex upper triangular
44*> matrix T (or of any matrix Q*T*Q**H with Q unitary).
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] JOB
51*> \verbatim
52*> JOB is CHARACTER*1
53*> Specifies whether condition numbers are required for
54*> eigenvalues (S) or eigenvectors (SEP):
55*> = 'E': for eigenvalues only (S);
56*> = 'V': for eigenvectors only (SEP);
57*> = 'B': for both eigenvalues and eigenvectors (S and SEP).
58*> \endverbatim
59*>
60*> \param[in] HOWMNY
61*> \verbatim
62*> HOWMNY is CHARACTER*1
63*> = 'A': compute condition numbers for all eigenpairs;
64*> = 'S': compute condition numbers for selected eigenpairs
65*> specified by the array SELECT.
66*> \endverbatim
67*>
68*> \param[in] SELECT
69*> \verbatim
70*> SELECT is LOGICAL array, dimension (N)
71*> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
72*> condition numbers are required. To select condition numbers
73*> for the j-th eigenpair, SELECT(j) must be set to .TRUE..
74*> If HOWMNY = 'A', SELECT is not referenced.
75*> \endverbatim
76*>
77*> \param[in] N
78*> \verbatim
79*> N is INTEGER
80*> The order of the matrix T. N >= 0.
81*> \endverbatim
82*>
83*> \param[in] T
84*> \verbatim
85*> T is COMPLEX*16 array, dimension (LDT,N)
86*> The upper triangular matrix T.
87*> \endverbatim
88*>
89*> \param[in] LDT
90*> \verbatim
91*> LDT is INTEGER
92*> The leading dimension of the array T. LDT >= max(1,N).
93*> \endverbatim
94*>
95*> \param[in] VL
96*> \verbatim
97*> VL is COMPLEX*16 array, dimension (LDVL,M)
98*> If JOB = 'E' or 'B', VL must contain left eigenvectors of T
99*> (or of any Q*T*Q**H with Q unitary), corresponding to the
100*> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
101*> must be stored in consecutive columns of VL, as returned by
102*> ZHSEIN or ZTREVC.
103*> If JOB = 'V', VL is not referenced.
104*> \endverbatim
105*>
106*> \param[in] LDVL
107*> \verbatim
108*> LDVL is INTEGER
109*> The leading dimension of the array VL.
110*> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
111*> \endverbatim
112*>
113*> \param[in] VR
114*> \verbatim
115*> VR is COMPLEX*16 array, dimension (LDVR,M)
116*> If JOB = 'E' or 'B', VR must contain right eigenvectors of T
117*> (or of any Q*T*Q**H with Q unitary), corresponding to the
118*> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
119*> must be stored in consecutive columns of VR, as returned by
120*> ZHSEIN or ZTREVC.
121*> If JOB = 'V', VR is not referenced.
122*> \endverbatim
123*>
124*> \param[in] LDVR
125*> \verbatim
126*> LDVR is INTEGER
127*> The leading dimension of the array VR.
128*> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
129*> \endverbatim
130*>
131*> \param[out] S
132*> \verbatim
133*> S is DOUBLE PRECISION array, dimension (MM)
134*> If JOB = 'E' or 'B', the reciprocal condition numbers of the
135*> selected eigenvalues, stored in consecutive elements of the
136*> array. Thus S(j), SEP(j), and the j-th columns of VL and VR
137*> all correspond to the same eigenpair (but not in general the
138*> j-th eigenpair, unless all eigenpairs are selected).
139*> If JOB = 'V', S is not referenced.
140*> \endverbatim
141*>
142*> \param[out] SEP
143*> \verbatim
144*> SEP is DOUBLE PRECISION array, dimension (MM)
145*> If JOB = 'V' or 'B', the estimated reciprocal condition
146*> numbers of the selected eigenvectors, stored in consecutive
147*> elements of the array.
148*> If JOB = 'E', SEP is not referenced.
149*> \endverbatim
150*>
151*> \param[in] MM
152*> \verbatim
153*> MM is INTEGER
154*> The number of elements in the arrays S (if JOB = 'E' or 'B')
155*> and/or SEP (if JOB = 'V' or 'B'). MM >= M.
156*> \endverbatim
157*>
158*> \param[out] M
159*> \verbatim
160*> M is INTEGER
161*> The number of elements of the arrays S and/or SEP actually
162*> used to store the estimated condition numbers.
163*> If HOWMNY = 'A', M is set to N.
164*> \endverbatim
165*>
166*> \param[out] WORK
167*> \verbatim
168*> WORK is COMPLEX*16 array, dimension (LDWORK,N+6)
169*> If JOB = 'E', WORK is not referenced.
170*> \endverbatim
171*>
172*> \param[in] LDWORK
173*> \verbatim
174*> LDWORK is INTEGER
175*> The leading dimension of the array WORK.
176*> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
177*> \endverbatim
178*>
179*> \param[out] RWORK
180*> \verbatim
181*> RWORK is DOUBLE PRECISION array, dimension (N)
182*> If JOB = 'E', RWORK is not referenced.
183*> \endverbatim
184*>
185*> \param[out] INFO
186*> \verbatim
187*> INFO is INTEGER
188*> = 0: successful exit
189*> < 0: if INFO = -i, the i-th argument had an illegal value
190*> \endverbatim
191*
192* Authors:
193* ========
194*
195*> \author Univ. of Tennessee
196*> \author Univ. of California Berkeley
197*> \author Univ. of Colorado Denver
198*> \author NAG Ltd.
199*
200*> \ingroup trsna
201*
202*> \par Further Details:
203* =====================
204*>
205*> \verbatim
206*>
207*> The reciprocal of the condition number of an eigenvalue lambda is
208*> defined as
209*>
210*> S(lambda) = |v**H*u| / (norm(u)*norm(v))
211*>
212*> where u and v are the right and left eigenvectors of T corresponding
213*> to lambda; v**H denotes the conjugate transpose of v, and norm(u)
214*> denotes the Euclidean norm. These reciprocal condition numbers always
215*> lie between zero (very badly conditioned) and one (very well
216*> conditioned). If n = 1, S(lambda) is defined to be 1.
217*>
218*> An approximate error bound for a computed eigenvalue W(i) is given by
219*>
220*> EPS * norm(T) / S(i)
221*>
222*> where EPS is the machine precision.
223*>
224*> The reciprocal of the condition number of the right eigenvector u
225*> corresponding to lambda is defined as follows. Suppose
226*>
227*> T = ( lambda c )
228*> ( 0 T22 )
229*>
230*> Then the reciprocal condition number is
231*>
232*> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
233*>
234*> where sigma-min denotes the smallest singular value. We approximate
235*> the smallest singular value by the reciprocal of an estimate of the
236*> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
237*> defined to be abs(T(1,1)).
238*>
239*> An approximate error bound for a computed right eigenvector VR(i)
240*> is given by
241*>
242*> EPS * norm(T) / SEP(i)
243*> \endverbatim
244*>
245* =====================================================================
246 SUBROUTINE ztrsna( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
247 $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
248 $ INFO )
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 DOUBLE PRECISION RWORK( * ), S( * ), SEP( * )
261 COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
262 $ work( ldwork, * )
263* ..
264*
265* =====================================================================
266*
267* .. Parameters ..
268 DOUBLE PRECISION ZERO, ONE
269 PARAMETER ( ZERO = 0.0d+0, one = 1.0d0+0 )
270* ..
271* .. Local Scalars ..
272 LOGICAL SOMCON, WANTBH, WANTS, WANTSP
273 CHARACTER NORMIN
274 INTEGER I, IERR, IX, J, K, KASE, KS
275 DOUBLE PRECISION BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
276 $ xnorm
277 COMPLEX*16 CDUM, PROD
278* ..
279* .. Local Arrays ..
280 INTEGER ISAVE( 3 )
281 COMPLEX*16 DUMMY( 1 )
282* ..
283* .. External Functions ..
284 LOGICAL LSAME
285 INTEGER IZAMAX
286 DOUBLE PRECISION DLAMCH, DZNRM2
287 COMPLEX*16 ZDOTC
288 EXTERNAL lsame, izamax, dlamch, dznrm2, zdotc
289* ..
290* .. External Subroutines ..
291 EXTERNAL xerbla, zdrscl, zlacn2, zlacpy, zlatrs, ztrexc
292* ..
293* .. Intrinsic Functions ..
294 INTRINSIC abs, dble, dimag, max
295* ..
296* .. Statement Functions ..
297 DOUBLE PRECISION CABS1
298* ..
299* .. Statement Function definitions ..
300 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
301* ..
302* .. Executable Statements ..
303*
304* Decode and test the input parameters
305*
306 wantbh = lsame( job, 'B' )
307 wants = lsame( job, 'E' ) .OR. wantbh
308 wantsp = lsame( job, 'V' ) .OR. wantbh
309*
310 somcon = lsame( howmny, 'S' )
311*
312* Set M to the number of eigenpairs for which condition numbers are
313* to be computed.
314*
315 IF( somcon ) THEN
316 m = 0
317 DO 10 j = 1, n
318 IF( SELECT( j ) )
319 $ m = m + 1
320 10 CONTINUE
321 ELSE
322 m = n
323 END IF
324*
325 info = 0
326 IF( .NOT.wants .AND. .NOT.wantsp ) THEN
327 info = -1
328 ELSE IF( .NOT.lsame( howmny, 'A' ) .AND. .NOT.somcon ) THEN
329 info = -2
330 ELSE IF( n.LT.0 ) THEN
331 info = -4
332 ELSE IF( ldt.LT.max( 1, n ) ) THEN
333 info = -6
334 ELSE IF( ldvl.LT.1 .OR. ( wants .AND. ldvl.LT.n ) ) THEN
335 info = -8
336 ELSE IF( ldvr.LT.1 .OR. ( wants .AND. ldvr.LT.n ) ) THEN
337 info = -10
338 ELSE IF( mm.LT.m ) THEN
339 info = -13
340 ELSE IF( ldwork.LT.1 .OR. ( wantsp .AND. ldwork.LT.n ) ) THEN
341 info = -16
342 END IF
343 IF( info.NE.0 ) THEN
344 CALL xerbla( 'ZTRSNA', -info )
345 RETURN
346 END IF
347*
348* Quick return if possible
349*
350 IF( n.EQ.0 )
351 $ RETURN
352*
353 IF( n.EQ.1 ) THEN
354 IF( somcon ) THEN
355 IF( .NOT.SELECT( 1 ) )
356 $ RETURN
357 END IF
358 IF( wants )
359 $ s( 1 ) = one
360 IF( wantsp )
361 $ sep( 1 ) = abs( t( 1, 1 ) )
362 RETURN
363 END IF
364*
365* Get machine constants
366*
367 eps = dlamch( 'P' )
368 smlnum = dlamch( 'S' ) / eps
369 bignum = one / smlnum
370*
371 ks = 1
372 DO 50 k = 1, n
373*
374 IF( somcon ) THEN
375 IF( .NOT.SELECT( k ) )
376 $ GO TO 50
377 END IF
378*
379 IF( wants ) THEN
380*
381* Compute the reciprocal condition number of the k-th
382* eigenvalue.
383*
384 prod = zdotc( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
385 rnrm = dznrm2( n, vr( 1, ks ), 1 )
386 lnrm = dznrm2( n, vl( 1, ks ), 1 )
387 s( ks ) = abs( prod ) / ( rnrm*lnrm )
388*
389 END IF
390*
391 IF( wantsp ) THEN
392*
393* Estimate the reciprocal condition number of the k-th
394* eigenvector.
395*
396* Copy the matrix T to the array WORK and swap the k-th
397* diagonal element to the (1,1) position.
398*
399 CALL zlacpy( 'Full', n, n, t, ldt, work, ldwork )
400 CALL ztrexc( 'No Q', n, work, ldwork, dummy, 1, k, 1, ierr )
401*
402* Form C = T22 - lambda*I in WORK(2:N,2:N).
403*
404 DO 20 i = 2, n
405 work( i, i ) = work( i, i ) - work( 1, 1 )
406 20 CONTINUE
407*
408* Estimate a lower bound for the 1-norm of inv(C**H). The 1st
409* and (N+1)th columns of WORK are used to store work vectors.
410*
411 sep( ks ) = zero
412 est = zero
413 kase = 0
414 normin = 'N'
415 30 CONTINUE
416 CALL zlacn2( n-1, work( 1, n+1 ), work, est, kase, isave )
417*
418 IF( kase.NE.0 ) THEN
419 IF( kase.EQ.1 ) THEN
420*
421* Solve C**H*x = scale*b
422*
423 CALL zlatrs( 'Upper', 'Conjugate transpose',
424 $ 'Nonunit', normin, n-1, work( 2, 2 ),
425 $ ldwork, work, scale, rwork, ierr )
426 ELSE
427*
428* Solve C*x = scale*b
429*
430 CALL zlatrs( 'Upper', 'No transpose', 'Nonunit',
431 $ normin, n-1, work( 2, 2 ), ldwork, work,
432 $ scale, rwork, ierr )
433 END IF
434 normin = 'Y'
435 IF( scale.NE.one ) THEN
436*
437* Multiply by 1/SCALE if doing so will not cause
438* overflow.
439*
440 ix = izamax( n-1, work, 1 )
441 xnorm = cabs1( work( ix, 1 ) )
442 IF( scale.LT.xnorm*smlnum .OR. scale.EQ.zero )
443 $ GO TO 40
444 CALL zdrscl( n, scale, work, 1 )
445 END IF
446 GO TO 30
447 END IF
448*
449 sep( ks ) = one / max( est, smlnum )
450 END IF
451*
452 40 CONTINUE
453 ks = ks + 1
454 50 CONTINUE
455 RETURN
456*
457* End of ZTRSNA
458*
459 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zlacn2
subroutine zlacn2(n, v, x, est, kase, isave)
ZLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition zlacn2.f:133
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
zlatrs
subroutine zlatrs(uplo, trans, diag, normin, n, a, lda, x, scale, cnorm, info)
ZLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition zlatrs.f:239
zdrscl
subroutine zdrscl(n, sa, sx, incx)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition zdrscl.f:84
ztrexc
subroutine ztrexc(compq, n, t, ldt, q, ldq, ifst, ilst, info)
ZTREXC
Definition ztrexc.f:126
ztrsna
subroutine ztrsna(job, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, s, sep, mm, m, work, ldwork, rwork, info)
ZTRSNA
Definition ztrsna.f:249