LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
ctrsna.f
Go to the documentation of this file.
1 *> \brief \b CTRSNA
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CTRSNA + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrsna.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrsna.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrsna.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CTRSNA( 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 * REAL RWORK( * ), S( * ), SEP( * )
32 * COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
33 * $ WORK( LDWORK, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> CTRSNA 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 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 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 *> CHSEIN or CTREVC.
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 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 *> CHSEIN or CTREVC.
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 REAL 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 REAL 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 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 REAL 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 complexOTHERcomputational
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 ctrsna( 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  REAL rwork( * ), s( * ), sep( * )
264  COMPLEX t( ldt, * ), vl( ldvl, * ), vr( ldvr, * ),
265  $ work( ldwork, * )
266 * ..
267 *
268 * =====================================================================
269 *
270 * .. Parameters ..
271  REAL zero, one
272  parameter( zero = 0.0e+0, one = 1.0+0 )
273 * ..
274 * .. Local Scalars ..
275  LOGICAL somcon, wantbh, wants, wantsp
276  CHARACTER normin
277  INTEGER i, ierr, ix, j, k, kase, ks
278  REAL bignum, eps, est, lnrm, rnrm, scale, smlnum,
279  $ xnorm
280  COMPLEX cdum, prod
281 * ..
282 * .. Local Arrays ..
283  INTEGER isave( 3 )
284  COMPLEX dummy( 1 )
285 * ..
286 * .. External Functions ..
287  LOGICAL lsame
288  INTEGER icamax
289  REAL scnrm2, slamch
290  COMPLEX cdotc
291  EXTERNAL lsame, icamax, scnrm2, slamch, cdotc
292 * ..
293 * .. External Subroutines ..
294  EXTERNAL clacn2, clacpy, clatrs, csrscl, ctrexc, slabad,
295  $ xerbla
296 * ..
297 * .. Intrinsic Functions ..
298  INTRINSIC abs, aimag, max, real
299 * ..
300 * .. Statement Functions ..
301  REAL cabs1
302 * ..
303 * .. Statement Function definitions ..
304  cabs1( cdum ) = abs( REAL( CDUM ) ) + abs( aimag( cdum ) )
305 * ..
306 * .. Executable Statements ..
307 *
308 * Decode and test the input parameters
309 *
310  wantbh = lsame( job, 'B' )
311  wants = lsame( job, 'E' ) .OR. wantbh
312  wantsp = lsame( job, 'V' ) .OR. wantbh
313 *
314  somcon = lsame( howmny, 'S' )
315 *
316 * Set M to the number of eigenpairs for which condition numbers are
317 * to be computed.
318 *
319  IF( somcon ) THEN
320  m = 0
321  DO 10 j = 1, n
322  IF( SELECT( j ) )
323  $ m = m + 1
324  10 continue
325  ELSE
326  m = n
327  END IF
328 *
329  info = 0
330  IF( .NOT.wants .AND. .NOT.wantsp ) THEN
331  info = -1
332  ELSE IF( .NOT.lsame( howmny, 'A' ) .AND. .NOT.somcon ) THEN
333  info = -2
334  ELSE IF( n.LT.0 ) THEN
335  info = -4
336  ELSE IF( ldt.LT.max( 1, n ) ) THEN
337  info = -6
338  ELSE IF( ldvl.LT.1 .OR. ( wants .AND. ldvl.LT.n ) ) THEN
339  info = -8
340  ELSE IF( ldvr.LT.1 .OR. ( wants .AND. ldvr.LT.n ) ) THEN
341  info = -10
342  ELSE IF( mm.LT.m ) THEN
343  info = -13
344  ELSE IF( ldwork.LT.1 .OR. ( wantsp .AND. ldwork.LT.n ) ) THEN
345  info = -16
346  END IF
347  IF( info.NE.0 ) THEN
348  CALL xerbla( 'CTRSNA', -info )
349  return
350  END IF
351 *
352 * Quick return if possible
353 *
354  IF( n.EQ.0 )
355  $ return
356 *
357  IF( n.EQ.1 ) THEN
358  IF( somcon ) THEN
359  IF( .NOT.SELECT( 1 ) )
360  $ return
361  END IF
362  IF( wants )
363  $ s( 1 ) = one
364  IF( wantsp )
365  $ sep( 1 ) = abs( t( 1, 1 ) )
366  return
367  END IF
368 *
369 * Get machine constants
370 *
371  eps = slamch( 'P' )
372  smlnum = slamch( 'S' ) / eps
373  bignum = one / smlnum
374  CALL slabad( smlnum, bignum )
375 *
376  ks = 1
377  DO 50 k = 1, n
378 *
379  IF( somcon ) THEN
380  IF( .NOT.SELECT( k ) )
381  $ go to 50
382  END IF
383 *
384  IF( wants ) THEN
385 *
386 * Compute the reciprocal condition number of the k-th
387 * eigenvalue.
388 *
389  prod = cdotc( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
390  rnrm = scnrm2( n, vr( 1, ks ), 1 )
391  lnrm = scnrm2( n, vl( 1, ks ), 1 )
392  s( ks ) = abs( prod ) / ( rnrm*lnrm )
393 *
394  END IF
395 *
396  IF( wantsp ) THEN
397 *
398 * Estimate the reciprocal condition number of the k-th
399 * eigenvector.
400 *
401 * Copy the matrix T to the array WORK and swap the k-th
402 * diagonal element to the (1,1) position.
403 *
404  CALL clacpy( 'Full', n, n, t, ldt, work, ldwork )
405  CALL ctrexc( 'No Q', n, work, ldwork, dummy, 1, k, 1, ierr )
406 *
407 * Form C = T22 - lambda*I in WORK(2:N,2:N).
408 *
409  DO 20 i = 2, n
410  work( i, i ) = work( i, i ) - work( 1, 1 )
411  20 continue
412 *
413 * Estimate a lower bound for the 1-norm of inv(C**H). The 1st
414 * and (N+1)th columns of WORK are used to store work vectors.
415 *
416  sep( ks ) = zero
417  est = zero
418  kase = 0
419  normin = 'N'
420  30 continue
421  CALL clacn2( n-1, work( 1, n+1 ), work, est, kase, isave )
422 *
423  IF( kase.NE.0 ) THEN
424  IF( kase.EQ.1 ) THEN
425 *
426 * Solve C**H*x = scale*b
427 *
428  CALL clatrs( 'Upper', 'Conjugate transpose',
429  $ 'Nonunit', normin, n-1, work( 2, 2 ),
430  $ ldwork, work, scale, rwork, ierr )
431  ELSE
432 *
433 * Solve C*x = scale*b
434 *
435  CALL clatrs( 'Upper', 'No transpose', 'Nonunit',
436  $ normin, n-1, work( 2, 2 ), ldwork, work,
437  $ scale, rwork, ierr )
438  END IF
439  normin = 'Y'
440  IF( scale.NE.one ) THEN
441 *
442 * Multiply by 1/SCALE if doing so will not cause
443 * overflow.
444 *
445  ix = icamax( n-1, work, 1 )
446  xnorm = cabs1( work( ix, 1 ) )
447  IF( scale.LT.xnorm*smlnum .OR. scale.EQ.zero )
448  $ go to 40
449  CALL csrscl( n, scale, work, 1 )
450  END IF
451  go to 30
452  END IF
453 *
454  sep( ks ) = one / max( est, smlnum )
455  END IF
456 *
457  40 continue
458  ks = ks + 1
459  50 continue
460  return
461 *
462 * End of CTRSNA
463 *
464  END