LAPACK 3.12.1
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*> Download ZTRSNA + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsna.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsna.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsna.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
20* LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
21* INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER HOWMNY, JOB
25* INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
26* ..
27* .. Array Arguments ..
28* LOGICAL SELECT( * )
29* DOUBLE PRECISION RWORK( * ), S( * ), SEP( * )
30* COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
31* $ WORK( LDWORK, * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> ZTRSNA estimates reciprocal condition numbers for specified
41*> eigenvalues and/or right eigenvectors of a complex upper triangular
42*> matrix T (or of any matrix Q*T*Q**H with Q unitary).
43*> \endverbatim
44*
45* Arguments:
46* ==========
47*
48*> \param[in] JOB
49*> \verbatim
50*> JOB is CHARACTER*1
51*> Specifies whether condition numbers are required for
52*> eigenvalues (S) or eigenvectors (SEP):
53*> = 'E': for eigenvalues only (S);
54*> = 'V': for eigenvectors only (SEP);
55*> = 'B': for both eigenvalues and eigenvectors (S and SEP).
56*> \endverbatim
57*>
58*> \param[in] HOWMNY
59*> \verbatim
60*> HOWMNY is CHARACTER*1
61*> = 'A': compute condition numbers for all eigenpairs;
62*> = 'S': compute condition numbers for selected eigenpairs
63*> specified by the array SELECT.
64*> \endverbatim
65*>
66*> \param[in] SELECT
67*> \verbatim
68*> SELECT is LOGICAL array, dimension (N)
69*> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
70*> condition numbers are required. To select condition numbers
71*> for the j-th eigenpair, SELECT(j) must be set to .TRUE..
72*> If HOWMNY = 'A', SELECT is not referenced.
73*> \endverbatim
74*>
75*> \param[in] N
76*> \verbatim
77*> N is INTEGER
78*> The order of the matrix T. N >= 0.
79*> \endverbatim
80*>
81*> \param[in] T
82*> \verbatim
83*> T is COMPLEX*16 array, dimension (LDT,N)
84*> The upper triangular matrix T.
85*> \endverbatim
86*>
87*> \param[in] LDT
88*> \verbatim
89*> LDT is INTEGER
90*> The leading dimension of the array T. LDT >= max(1,N).
91*> \endverbatim
92*>
93*> \param[in] VL
94*> \verbatim
95*> VL is COMPLEX*16 array, dimension (LDVL,M)
96*> If JOB = 'E' or 'B', VL must contain left eigenvectors of T
97*> (or of any Q*T*Q**H with Q unitary), corresponding to the
98*> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
99*> must be stored in consecutive columns of VL, as returned by
100*> ZHSEIN or ZTREVC.
101*> If JOB = 'V', VL is not referenced.
102*> \endverbatim
103*>
104*> \param[in] LDVL
105*> \verbatim
106*> LDVL is INTEGER
107*> The leading dimension of the array VL.
108*> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
109*> \endverbatim
110*>
111*> \param[in] VR
112*> \verbatim
113*> VR is COMPLEX*16 array, dimension (LDVR,M)
114*> If JOB = 'E' or 'B', VR must contain right eigenvectors of T
115*> (or of any Q*T*Q**H with Q unitary), corresponding to the
116*> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
117*> must be stored in consecutive columns of VR, as returned by
118*> ZHSEIN or ZTREVC.
119*> If JOB = 'V', VR is not referenced.
120*> \endverbatim
121*>
122*> \param[in] LDVR
123*> \verbatim
124*> LDVR is INTEGER
125*> The leading dimension of the array VR.
126*> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
127*> \endverbatim
128*>
129*> \param[out] S
130*> \verbatim
131*> S is DOUBLE PRECISION array, dimension (MM)
132*> If JOB = 'E' or 'B', the reciprocal condition numbers of the
133*> selected eigenvalues, stored in consecutive elements of the
134*> array. Thus S(j), SEP(j), and the j-th columns of VL and VR
135*> all correspond to the same eigenpair (but not in general the
136*> j-th eigenpair, unless all eigenpairs are selected).
137*> If JOB = 'V', S is not referenced.
138*> \endverbatim
139*>
140*> \param[out] SEP
141*> \verbatim
142*> SEP is DOUBLE PRECISION array, dimension (MM)
143*> If JOB = 'V' or 'B', the estimated reciprocal condition
144*> numbers of the selected eigenvectors, stored in consecutive
145*> elements of the array.
146*> If JOB = 'E', SEP is not referenced.
147*> \endverbatim
148*>
149*> \param[in] MM
150*> \verbatim
151*> MM is INTEGER
152*> The number of elements in the arrays S (if JOB = 'E' or 'B')
153*> and/or SEP (if JOB = 'V' or 'B'). MM >= M.
154*> \endverbatim
155*>
156*> \param[out] M
157*> \verbatim
158*> M is INTEGER
159*> The number of elements of the arrays S and/or SEP actually
160*> used to store the estimated condition numbers.
161*> If HOWMNY = 'A', M is set to N.
162*> \endverbatim
163*>
164*> \param[out] WORK
165*> \verbatim
166*> WORK is COMPLEX*16 array, dimension (LDWORK,N+6)
167*> If JOB = 'E', WORK is not referenced.
168*> \endverbatim
169*>
170*> \param[in] LDWORK
171*> \verbatim
172*> LDWORK is INTEGER
173*> The leading dimension of the array WORK.
174*> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
175*> \endverbatim
176*>
177*> \param[out] RWORK
178*> \verbatim
179*> RWORK is DOUBLE PRECISION array, dimension (N)
180*> If JOB = 'E', RWORK is not referenced.
181*> \endverbatim
182*>
183*> \param[out] INFO
184*> \verbatim
185*> INFO is INTEGER
186*> = 0: successful exit
187*> < 0: if INFO = -i, the i-th argument had an illegal value
188*> \endverbatim
189*
190* Authors:
191* ========
192*
193*> \author Univ. of Tennessee
194*> \author Univ. of California Berkeley
195*> \author Univ. of Colorado Denver
196*> \author NAG Ltd.
197*
198*> \ingroup trsna
199*
200*> \par Further Details:
201* =====================
202*>
203*> \verbatim
204*>
205*> The reciprocal of the condition number of an eigenvalue lambda is
206*> defined as
207*>
208*> S(lambda) = |v**H*u| / (norm(u)*norm(v))
209*>
210*> where u and v are the right and left eigenvectors of T corresponding
211*> to lambda; v**H denotes the conjugate transpose of v, and norm(u)
212*> denotes the Euclidean norm. These reciprocal condition numbers always
213*> lie between zero (very badly conditioned) and one (very well
214*> conditioned). If n = 1, S(lambda) is defined to be 1.
215*>
216*> An approximate error bound for a computed eigenvalue W(i) is given by
217*>
218*> EPS * norm(T) / S(i)
219*>
220*> where EPS is the machine precision.
221*>
222*> The reciprocal of the condition number of the right eigenvector u
223*> corresponding to lambda is defined as follows. Suppose
224*>
225*> T = ( lambda c )
226*> ( 0 T22 )
227*>
228*> Then the reciprocal condition number is
229*>
230*> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
231*>
232*> where sigma-min denotes the smallest singular value. We approximate
233*> the smallest singular value by the reciprocal of an estimate of the
234*> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
235*> defined to be abs(T(1,1)).
236*>
237*> An approximate error bound for a computed right eigenvector VR(i)
238*> is given by
239*>
240*> EPS * norm(T) / SEP(i)
241*> \endverbatim
242*>
243* =====================================================================
244 SUBROUTINE ztrsna( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
245 $ VR,
246 $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
247 $ INFO )
248*
249* -- LAPACK computational routine --
250* -- LAPACK is a software package provided by Univ. of Tennessee, --
251* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
252*
253* .. Scalar Arguments ..
254 CHARACTER HOWMNY, JOB
255 INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
256* ..
257* .. Array Arguments ..
258 LOGICAL SELECT( * )
259 DOUBLE PRECISION RWORK( * ), S( * ), SEP( * )
260 COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
261 $ work( ldwork, * )
262* ..
263*
264* =====================================================================
265*
266* .. Parameters ..
267 DOUBLE PRECISION ZERO, ONE
268 PARAMETER ( ZERO = 0.0d+0, one = 1.0d0+0 )
269* ..
270* .. Local Scalars ..
271 LOGICAL SOMCON, WANTBH, WANTS, WANTSP
272 CHARACTER NORMIN
273 INTEGER I, IERR, IX, J, K, KASE, KS
274 DOUBLE PRECISION BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
275 $ xnorm
276 COMPLEX*16 CDUM, PROD
277* ..
278* .. Local Arrays ..
279 INTEGER ISAVE( 3 )
280 COMPLEX*16 DUMMY( 1 )
281* ..
282* .. External Functions ..
283 LOGICAL LSAME
284 INTEGER IZAMAX
285 DOUBLE PRECISION DLAMCH, DZNRM2
286 COMPLEX*16 ZDOTC
287 EXTERNAL lsame, izamax, dlamch, dznrm2,
288 $ zdotc
289* ..
290* .. External Subroutines ..
291 EXTERNAL xerbla, zdrscl, zlacn2, zlacpy, zlatrs,
292 $ ztrexc
293* ..
294* .. Intrinsic Functions ..
295 INTRINSIC abs, dble, dimag, max
296* ..
297* .. Statement Functions ..
298 DOUBLE PRECISION CABS1
299* ..
300* .. Statement Function definitions ..
301 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( 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( 'ZTRSNA', -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 = dlamch( 'P' )
369 smlnum = dlamch( 'S' ) / eps
370 bignum = one / smlnum
371*
372 ks = 1
373 DO 50 k = 1, n
374*
375 IF( somcon ) THEN
376 IF( .NOT.SELECT( k ) )
377 $ GO TO 50
378 END IF
379*
380 IF( wants ) THEN
381*
382* Compute the reciprocal condition number of the k-th
383* eigenvalue.
384*
385 prod = zdotc( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
386 rnrm = dznrm2( n, vr( 1, ks ), 1 )
387 lnrm = dznrm2( n, vl( 1, ks ), 1 )
388 s( ks ) = abs( prod ) / ( rnrm*lnrm )
389*
390 END IF
391*
392 IF( wantsp ) THEN
393*
394* Estimate the reciprocal condition number of the k-th
395* eigenvector.
396*
397* Copy the matrix T to the array WORK and swap the k-th
398* diagonal element to the (1,1) position.
399*
400 CALL zlacpy( 'Full', n, n, t, ldt, work, ldwork )
401 CALL ztrexc( 'No Q', n, work, ldwork, dummy, 1, k, 1,
402 $ 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 zlacn2( n-1, work( 1, n+1 ), work, est, kase,
419 $ isave )
420*
421 IF( kase.NE.0 ) THEN
422 IF( kase.EQ.1 ) THEN
423*
424* Solve C**H*x = scale*b
425*
426 CALL zlatrs( 'Upper', 'Conjugate transpose',
427 $ 'Nonunit', normin, n-1, work( 2, 2 ),
428 $ ldwork, work, scale, rwork, ierr )
429 ELSE
430*
431* Solve C*x = scale*b
432*
433 CALL zlatrs( 'Upper', 'No transpose', 'Nonunit',
434 $ normin, n-1, work( 2, 2 ), ldwork, work,
435 $ scale, rwork, ierr )
436 END IF
437 normin = 'Y'
438 IF( scale.NE.one ) THEN
439*
440* Multiply by 1/SCALE if doing so will not cause
441* overflow.
442*
443 ix = izamax( n-1, work, 1 )
444 xnorm = cabs1( work( ix, 1 ) )
445 IF( scale.LT.xnorm*smlnum .OR. scale.EQ.zero )
446 $ GO TO 40
447 CALL zdrscl( n, scale, work, 1 )
448 END IF
449 GO TO 30
450 END IF
451*
452 sep( ks ) = one / max( est, smlnum )
453 END IF
454*
455 40 CONTINUE
456 ks = ks + 1
457 50 CONTINUE
458 RETURN
459*
460* End of ZTRSNA
461*
462 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:131
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:101
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:238
zdrscl
subroutine zdrscl(n, sa, sx, incx)
ZDRSCL multiplies a vector by the reciprocal of a real scalar.
Definition zdrscl.f:82
ztrexc
subroutine ztrexc(compq, n, t, ldt, q, ldq, ifst, ilst, info)
ZTREXC
Definition ztrexc.f:124
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:248