LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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
9 *> Download DTRSNA + dependencies
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 *> \date November 2011
217 *
218 *> \ingroup doubleOTHERcomputational
219 *
220 *> \par Further Details:
221 * =====================
222 *>
223 *> \verbatim
224 *>
225 *> The reciprocal of the condition number of an eigenvalue lambda is
226 *> defined as
227 *>
228 *> S(lambda) = |v**T*u| / (norm(u)*norm(v))
229 *>
230 *> where u and v are the right and left eigenvectors of T corresponding
231 *> to lambda; v**T denotes the transpose of v, and norm(u)
232 *> denotes the Euclidean norm. These reciprocal condition numbers always
233 *> lie between zero (very badly conditioned) and one (very well
234 *> conditioned). If n = 1, S(lambda) is defined to be 1.
235 *>
236 *> An approximate error bound for a computed eigenvalue W(i) is given by
237 *>
238 *> EPS * norm(T) / S(i)
239 *>
240 *> where EPS is the machine precision.
241 *>
242 *> The reciprocal of the condition number of the right eigenvector u
243 *> corresponding to lambda is defined as follows. Suppose
244 *>
245 *> T = ( lambda c )
246 *> ( 0 T22 )
247 *>
248 *> Then the reciprocal condition number is
249 *>
250 *> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
251 *>
252 *> where sigma-min denotes the smallest singular value. We approximate
253 *> the smallest singular value by the reciprocal of an estimate of the
254 *> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
255 *> defined to be abs(T(1,1)).
256 *>
257 *> An approximate error bound for a computed right eigenvector VR(i)
258 *> is given by
259 *>
260 *> EPS * norm(T) / SEP(i)
261 *> \endverbatim
262 *>
263 * =====================================================================
264  SUBROUTINE dtrsna( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
265  $ ldvr, s, sep, mm, m, work, ldwork, iwork,
266  $ info )
267 *
268 * -- LAPACK computational routine (version 3.4.0) --
269 * -- LAPACK is a software package provided by Univ. of Tennessee, --
270 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
271 * November 2011
272 *
273 * .. Scalar Arguments ..
274  CHARACTER howmny, job
275  INTEGER info, ldt, ldvl, ldvr, ldwork, m, mm, n
276 * ..
277 * .. Array Arguments ..
278  LOGICAL select( * )
279  INTEGER iwork( * )
280  DOUBLE PRECISION s( * ), sep( * ), t( ldt, * ), vl( ldvl, * ),
281  $ vr( ldvr, * ), work( ldwork, * )
282 * ..
283 *
284 * =====================================================================
285 *
286 * .. Parameters ..
287  DOUBLE PRECISION zero, one, two
288  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
289 * ..
290 * .. Local Scalars ..
291  LOGICAL pair, somcon, wantbh, wants, wantsp
292  INTEGER i, ierr, ifst, ilst, j, k, kase, ks, n2, nn
293  DOUBLE PRECISION bignum, cond, cs, delta, dumm, eps, est, lnrm,
294  $ mu, prod, prod1, prod2, rnrm, scale, smlnum, sn
295 * ..
296 * .. Local Arrays ..
297  INTEGER isave( 3 )
298  DOUBLE PRECISION dummy( 1 )
299 * ..
300 * .. External Functions ..
301  LOGICAL lsame
302  DOUBLE PRECISION ddot, dlamch, dlapy2, dnrm2
303  EXTERNAL lsame, ddot, dlamch, dlapy2, dnrm2
304 * ..
305 * .. External Subroutines ..
306  EXTERNAL dlacn2, dlacpy, dlaqtr, dtrexc, xerbla
307 * ..
308 * .. Intrinsic Functions ..
309  INTRINSIC abs, max, sqrt
310 * ..
311 * .. Executable Statements ..
312 *
313 * Decode and test the input parameters
314 *
315  wantbh = lsame( job, 'B' )
316  wants = lsame( job, 'E' ) .OR. wantbh
317  wantsp = lsame( job, 'V' ) .OR. wantbh
318 *
319  somcon = lsame( howmny, 'S' )
320 *
321  info = 0
322  IF( .NOT.wants .AND. .NOT.wantsp ) THEN
323  info = -1
324  ELSE IF( .NOT.lsame( howmny, 'A' ) .AND. .NOT.somcon ) THEN
325  info = -2
326  ELSE IF( n.LT.0 ) THEN
327  info = -4
328  ELSE IF( ldt.LT.max( 1, n ) ) THEN
329  info = -6
330  ELSE IF( ldvl.LT.1 .OR. ( wants .AND. ldvl.LT.n ) ) THEN
331  info = -8
332  ELSE IF( ldvr.LT.1 .OR. ( wants .AND. ldvr.LT.n ) ) THEN
333  info = -10
334  ELSE
335 *
336 * Set M to the number of eigenpairs for which condition numbers
337 * are required, and test MM.
338 *
339  IF( somcon ) THEN
340  m = 0
341  pair = .false.
342  DO 10 k = 1, n
343  IF( pair ) THEN
344  pair = .false.
345  ELSE
346  IF( k.LT.n ) THEN
347  IF( t( k+1, k ).EQ.zero ) THEN
348  IF( SELECT( k ) )
349  $ m = m + 1
350  ELSE
351  pair = .true.
352  IF( SELECT( k ) .OR. SELECT( k+1 ) )
353  $ m = m + 2
354  END IF
355  ELSE
356  IF( SELECT( n ) )
357  $ m = m + 1
358  END IF
359  END IF
360  10 continue
361  ELSE
362  m = n
363  END IF
364 *
365  IF( mm.LT.m ) THEN
366  info = -13
367  ELSE IF( ldwork.LT.1 .OR. ( wantsp .AND. ldwork.LT.n ) ) THEN
368  info = -16
369  END IF
370  END IF
371  IF( info.NE.0 ) THEN
372  CALL xerbla( 'DTRSNA', -info )
373  return
374  END IF
375 *
376 * Quick return if possible
377 *
378  IF( n.EQ.0 )
379  $ return
380 *
381  IF( n.EQ.1 ) THEN
382  IF( somcon ) THEN
383  IF( .NOT.SELECT( 1 ) )
384  $ return
385  END IF
386  IF( wants )
387  $ s( 1 ) = one
388  IF( wantsp )
389  $ sep( 1 ) = abs( t( 1, 1 ) )
390  return
391  END IF
392 *
393 * Get machine constants
394 *
395  eps = dlamch( 'P' )
396  smlnum = dlamch( 'S' ) / eps
397  bignum = one / smlnum
398  CALL dlabad( smlnum, bignum )
399 *
400  ks = 0
401  pair = .false.
402  DO 60 k = 1, n
403 *
404 * Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
405 *
406  IF( pair ) THEN
407  pair = .false.
408  go to 60
409  ELSE
410  IF( k.LT.n )
411  $ pair = t( k+1, k ).NE.zero
412  END IF
413 *
414 * Determine whether condition numbers are required for the k-th
415 * eigenpair.
416 *
417  IF( somcon ) THEN
418  IF( pair ) THEN
419  IF( .NOT.SELECT( k ) .AND. .NOT.SELECT( k+1 ) )
420  $ go to 60
421  ELSE
422  IF( .NOT.SELECT( k ) )
423  $ go to 60
424  END IF
425  END IF
426 *
427  ks = ks + 1
428 *
429  IF( wants ) THEN
430 *
431 * Compute the reciprocal condition number of the k-th
432 * eigenvalue.
433 *
434  IF( .NOT.pair ) THEN
435 *
436 * Real eigenvalue.
437 *
438  prod = ddot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
439  rnrm = dnrm2( n, vr( 1, ks ), 1 )
440  lnrm = dnrm2( n, vl( 1, ks ), 1 )
441  s( ks ) = abs( prod ) / ( rnrm*lnrm )
442  ELSE
443 *
444 * Complex eigenvalue.
445 *
446  prod1 = ddot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
447  prod1 = prod1 + ddot( n, vr( 1, ks+1 ), 1, vl( 1, ks+1 ),
448  $ 1 )
449  prod2 = ddot( n, vl( 1, ks ), 1, vr( 1, ks+1 ), 1 )
450  prod2 = prod2 - ddot( n, vl( 1, ks+1 ), 1, vr( 1, ks ),
451  $ 1 )
452  rnrm = dlapy2( dnrm2( n, vr( 1, ks ), 1 ),
453  $ dnrm2( n, vr( 1, ks+1 ), 1 ) )
454  lnrm = dlapy2( dnrm2( n, vl( 1, ks ), 1 ),
455  $ dnrm2( n, vl( 1, ks+1 ), 1 ) )
456  cond = dlapy2( prod1, prod2 ) / ( rnrm*lnrm )
457  s( ks ) = cond
458  s( ks+1 ) = cond
459  END IF
460  END IF
461 *
462  IF( wantsp ) THEN
463 *
464 * Estimate the reciprocal condition number of the k-th
465 * eigenvector.
466 *
467 * Copy the matrix T to the array WORK and swap the diagonal
468 * block beginning at T(k,k) to the (1,1) position.
469 *
470  CALL dlacpy( 'Full', n, n, t, ldt, work, ldwork )
471  ifst = k
472  ilst = 1
473  CALL dtrexc( 'No Q', n, work, ldwork, dummy, 1, ifst, ilst,
474  $ work( 1, n+1 ), ierr )
475 *
476  IF( ierr.EQ.1 .OR. ierr.EQ.2 ) THEN
477 *
478 * Could not swap because blocks not well separated
479 *
480  scale = one
481  est = bignum
482  ELSE
483 *
484 * Reordering successful
485 *
486  IF( work( 2, 1 ).EQ.zero ) THEN
487 *
488 * Form C = T22 - lambda*I in WORK(2:N,2:N).
489 *
490  DO 20 i = 2, n
491  work( i, i ) = work( i, i ) - work( 1, 1 )
492  20 continue
493  n2 = 1
494  nn = n - 1
495  ELSE
496 *
497 * Triangularize the 2 by 2 block by unitary
498 * transformation U = [ cs i*ss ]
499 * [ i*ss cs ].
500 * such that the (1,1) position of WORK is complex
501 * eigenvalue lambda with positive imaginary part. (2,2)
502 * position of WORK is the complex eigenvalue lambda
503 * with negative imaginary part.
504 *
505  mu = sqrt( abs( work( 1, 2 ) ) )*
506  $ sqrt( abs( work( 2, 1 ) ) )
507  delta = dlapy2( mu, work( 2, 1 ) )
508  cs = mu / delta
509  sn = -work( 2, 1 ) / delta
510 *
511 * Form
512 *
513 * C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
514 * [ mu ]
515 * [ .. ]
516 * [ .. ]
517 * [ mu ]
518 * where C**T is transpose of matrix C,
519 * and RWORK is stored starting in the N+1-st column of
520 * WORK.
521 *
522  DO 30 j = 3, n
523  work( 2, j ) = cs*work( 2, j )
524  work( j, j ) = work( j, j ) - work( 1, 1 )
525  30 continue
526  work( 2, 2 ) = zero
527 *
528  work( 1, n+1 ) = two*mu
529  DO 40 i = 2, n - 1
530  work( i, n+1 ) = sn*work( 1, i+1 )
531  40 continue
532  n2 = 2
533  nn = 2*( n-1 )
534  END IF
535 *
536 * Estimate norm(inv(C**T))
537 *
538  est = zero
539  kase = 0
540  50 continue
541  CALL dlacn2( nn, work( 1, n+2 ), work( 1, n+4 ), iwork,
542  $ est, kase, isave )
543  IF( kase.NE.0 ) THEN
544  IF( kase.EQ.1 ) THEN
545  IF( n2.EQ.1 ) THEN
546 *
547 * Real eigenvalue: solve C**T*x = scale*c.
548 *
549  CALL dlaqtr( .true., .true., n-1, work( 2, 2 ),
550  $ ldwork, dummy, dumm, scale,
551  $ work( 1, n+4 ), work( 1, n+6 ),
552  $ ierr )
553  ELSE
554 *
555 * Complex eigenvalue: solve
556 * C**T*(p+iq) = scale*(c+id) in real arithmetic.
557 *
558  CALL dlaqtr( .true., .false., n-1, work( 2, 2 ),
559  $ ldwork, work( 1, n+1 ), mu, scale,
560  $ work( 1, n+4 ), work( 1, n+6 ),
561  $ ierr )
562  END IF
563  ELSE
564  IF( n2.EQ.1 ) THEN
565 *
566 * Real eigenvalue: solve C*x = scale*c.
567 *
568  CALL dlaqtr( .false., .true., n-1, work( 2, 2 ),
569  $ ldwork, dummy, dumm, scale,
570  $ work( 1, n+4 ), work( 1, n+6 ),
571  $ ierr )
572  ELSE
573 *
574 * Complex eigenvalue: solve
575 * C*(p+iq) = scale*(c+id) in real arithmetic.
576 *
577  CALL dlaqtr( .false., .false., n-1,
578  $ work( 2, 2 ), ldwork,
579  $ work( 1, n+1 ), mu, scale,
580  $ work( 1, n+4 ), work( 1, n+6 ),
581  $ ierr )
582 *
583  END IF
584  END IF
585 *
586  go to 50
587  END IF
588  END IF
589 *
590  sep( ks ) = scale / max( est, smlnum )
591  IF( pair )
592  $ sep( ks+1 ) = sep( ks )
593  END IF
594 *
595  IF( pair )
596  $ ks = ks + 1
597 *
598  60 continue
599  return
600 *
601 * End of DTRSNA
602 *
603  END