LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
201 *
202 *> \ingroup complex16OTHERcomputational
203 *
204 *> \par Further Details:
205 * =====================
206 *>
207 *> \verbatim
208 *>
209 *> The reciprocal of the condition number of an eigenvalue lambda is
210 *> defined as
211 *>
212 *> S(lambda) = |v**H*u| / (norm(u)*norm(v))
213 *>
214 *> where u and v are the right and left eigenvectors of T corresponding
215 *> to lambda; v**H denotes the conjugate transpose of v, and norm(u)
216 *> denotes the Euclidean norm. These reciprocal condition numbers always
217 *> lie between zero (very badly conditioned) and one (very well
218 *> conditioned). If n = 1, S(lambda) is defined to be 1.
219 *>
220 *> An approximate error bound for a computed eigenvalue W(i) is given by
221 *>
222 *> EPS * norm(T) / S(i)
223 *>
224 *> where EPS is the machine precision.
225 *>
226 *> The reciprocal of the condition number of the right eigenvector u
227 *> corresponding to lambda is defined as follows. Suppose
228 *>
229 *> T = ( lambda c )
230 *> ( 0 T22 )
231 *>
232 *> Then the reciprocal condition number is
233 *>
234 *> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
235 *>
236 *> where sigma-min denotes the smallest singular value. We approximate
237 *> the smallest singular value by the reciprocal of an estimate of the
238 *> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
239 *> defined to be abs(T(1,1)).
240 *>
241 *> An approximate error bound for a computed right eigenvector VR(i)
242 *> is given by
243 *>
244 *> EPS * norm(T) / SEP(i)
245 *> \endverbatim
246 *>
247 * =====================================================================
248  SUBROUTINE ztrsna( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
249  $ ldvr, s, sep, mm, m, work, ldwork, rwork,
250  $ info )
251 *
252 * -- LAPACK computational routine (version 3.4.0) --
253 * -- LAPACK is a software package provided by Univ. of Tennessee, --
254 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
255 * November 2011
256 *
257 * .. Scalar Arguments ..
258  CHARACTER howmny, job
259  INTEGER info, ldt, ldvl, ldvr, ldwork, m, mm, n
260 * ..
261 * .. Array Arguments ..
262  LOGICAL select( * )
263  DOUBLE PRECISION rwork( * ), s( * ), sep( * )
264  COMPLEX*16 t( ldt, * ), vl( ldvl, * ), vr( ldvr, * ),
265  $ work( ldwork, * )
266 * ..
267 *
268 * =====================================================================
269 *
270 * .. Parameters ..
271  DOUBLE PRECISION zero, one
272  parameter( zero = 0.0d+0, one = 1.0d0+0 )
273 * ..
274 * .. Local Scalars ..
275  LOGICAL somcon, wantbh, wants, wantsp
276  CHARACTER normin
277  INTEGER i, ierr, ix, j, k, kase, ks
278  DOUBLE PRECISION bignum, eps, est, lnrm, rnrm, scale, smlnum,
279  $ xnorm
280  COMPLEX*16 cdum, prod
281 * ..
282 * .. Local Arrays ..
283  INTEGER isave( 3 )
284  COMPLEX*16 dummy( 1 )
285 * ..
286 * .. External Functions ..
287  LOGICAL lsame
288  INTEGER izamax
289  DOUBLE PRECISION dlamch, dznrm2
290  COMPLEX*16 zdotc
291  EXTERNAL lsame, izamax, dlamch, dznrm2, zdotc
292 * ..
293 * .. External Subroutines ..
294  EXTERNAL xerbla, zdrscl, zlacn2, zlacpy, zlatrs, ztrexc
295 * ..
296 * .. Intrinsic Functions ..
297  INTRINSIC abs, dble, dimag, max
298 * ..
299 * .. Statement Functions ..
300  DOUBLE PRECISION cabs1
301 * ..
302 * .. Statement Function definitions ..
303  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
304 * ..
305 * .. Executable Statements ..
306 *
307 * Decode and test the input parameters
308 *
309  wantbh = lsame( job, 'B' )
310  wants = lsame( job, 'E' ) .OR. wantbh
311  wantsp = lsame( job, 'V' ) .OR. wantbh
312 *
313  somcon = lsame( howmny, 'S' )
314 *
315 * Set M to the number of eigenpairs for which condition numbers are
316 * to be computed.
317 *
318  IF( somcon ) THEN
319  m = 0
320  DO 10 j = 1, n
321  IF( SELECT( j ) )
322  $ m = m + 1
323  10 continue
324  ELSE
325  m = n
326  END IF
327 *
328  info = 0
329  IF( .NOT.wants .AND. .NOT.wantsp ) THEN
330  info = -1
331  ELSE IF( .NOT.lsame( howmny, 'A' ) .AND. .NOT.somcon ) THEN
332  info = -2
333  ELSE IF( n.LT.0 ) THEN
334  info = -4
335  ELSE IF( ldt.LT.max( 1, n ) ) THEN
336  info = -6
337  ELSE IF( ldvl.LT.1 .OR. ( wants .AND. ldvl.LT.n ) ) THEN
338  info = -8
339  ELSE IF( ldvr.LT.1 .OR. ( wants .AND. ldvr.LT.n ) ) THEN
340  info = -10
341  ELSE IF( mm.LT.m ) THEN
342  info = -13
343  ELSE IF( ldwork.LT.1 .OR. ( wantsp .AND. ldwork.LT.n ) ) THEN
344  info = -16
345  END IF
346  IF( info.NE.0 ) THEN
347  CALL xerbla( 'ZTRSNA', -info )
348  return
349  END IF
350 *
351 * Quick return if possible
352 *
353  IF( n.EQ.0 )
354  $ return
355 *
356  IF( n.EQ.1 ) THEN
357  IF( somcon ) THEN
358  IF( .NOT.SELECT( 1 ) )
359  $ return
360  END IF
361  IF( wants )
362  $ s( 1 ) = one
363  IF( wantsp )
364  $ sep( 1 ) = abs( t( 1, 1 ) )
365  return
366  END IF
367 *
368 * Get machine constants
369 *
370  eps = dlamch( 'P' )
371  smlnum = dlamch( 'S' ) / eps
372  bignum = one / smlnum
373  CALL dlabad( smlnum, bignum )
374 *
375  ks = 1
376  DO 50 k = 1, n
377 *
378  IF( somcon ) THEN
379  IF( .NOT.SELECT( k ) )
380  $ go to 50
381  END IF
382 *
383  IF( wants ) THEN
384 *
385 * Compute the reciprocal condition number of the k-th
386 * eigenvalue.
387 *
388  prod = zdotc( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
389  rnrm = dznrm2( n, vr( 1, ks ), 1 )
390  lnrm = dznrm2( n, vl( 1, ks ), 1 )
391  s( ks ) = abs( prod ) / ( rnrm*lnrm )
392 *
393  END IF
394 *
395  IF( wantsp ) THEN
396 *
397 * Estimate the reciprocal condition number of the k-th
398 * eigenvector.
399 *
400 * Copy the matrix T to the array WORK and swap the k-th
401 * diagonal element to the (1,1) position.
402 *
403  CALL zlacpy( 'Full', n, n, t, ldt, work, ldwork )
404  CALL ztrexc( 'No Q', n, work, ldwork, dummy, 1, k, 1, ierr )
405 *
406 * Form C = T22 - lambda*I in WORK(2:N,2:N).
407 *
408  DO 20 i = 2, n
409  work( i, i ) = work( i, i ) - work( 1, 1 )
410  20 continue
411 *
412 * Estimate a lower bound for the 1-norm of inv(C**H). The 1st
413 * and (N+1)th columns of WORK are used to store work vectors.
414 *
415  sep( ks ) = zero
416  est = zero
417  kase = 0
418  normin = 'N'
419  30 continue
420  CALL zlacn2( n-1, work( 1, n+1 ), work, est, kase, isave )
421 *
422  IF( kase.NE.0 ) THEN
423  IF( kase.EQ.1 ) THEN
424 *
425 * Solve C**H*x = scale*b
426 *
427  CALL zlatrs( 'Upper', 'Conjugate transpose',
428  $ 'Nonunit', normin, n-1, work( 2, 2 ),
429  $ ldwork, work, scale, rwork, ierr )
430  ELSE
431 *
432 * Solve C*x = scale*b
433 *
434  CALL zlatrs( 'Upper', 'No transpose', 'Nonunit',
435  $ normin, n-1, work( 2, 2 ), ldwork, work,
436  $ scale, rwork, ierr )
437  END IF
438  normin = 'Y'
439  IF( scale.NE.one ) THEN
440 *
441 * Multiply by 1/SCALE if doing so will not cause
442 * overflow.
443 *
444  ix = izamax( n-1, work, 1 )
445  xnorm = cabs1( work( ix, 1 ) )
446  IF( scale.LT.xnorm*smlnum .OR. scale.EQ.zero )
447  $ go to 40
448  CALL zdrscl( n, scale, work, 1 )
449  END IF
450  go to 30
451  END IF
452 *
453  sep( ks ) = one / max( est, smlnum )
454  END IF
455 *
456  40 continue
457  ks = ks + 1
458  50 continue
459  return
460 *
461 * End of ZTRSNA
462 *
463  END