LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dtgsna.f
Go to the documentation of this file.
1*> \brief \b DTGSNA
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DTGSNA + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsna.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsna.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsna.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
22* LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
23* IWORK, INFO )
24*
25* .. Scalar Arguments ..
26* CHARACTER HOWMNY, JOB
27* INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
28* ..
29* .. Array Arguments ..
30* LOGICAL SELECT( * )
31* INTEGER IWORK( * )
32* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
33* $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> DTGSNA estimates reciprocal condition numbers for specified
43*> eigenvalues and/or eigenvectors of a matrix pair (A, B) in
44*> generalized real Schur canonical form (or of any matrix pair
45*> (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where
46*> Z**T denotes the transpose of Z.
47*>
48*> (A, B) must be in generalized real Schur form (as returned by DGGES),
49*> i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
50*> blocks. B is upper triangular.
51*>
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] JOB
58*> \verbatim
59*> JOB is CHARACTER*1
60*> Specifies whether condition numbers are required for
61*> eigenvalues (S) or eigenvectors (DIF):
62*> = 'E': for eigenvalues only (S);
63*> = 'V': for eigenvectors only (DIF);
64*> = 'B': for both eigenvalues and eigenvectors (S and DIF).
65*> \endverbatim
66*>
67*> \param[in] HOWMNY
68*> \verbatim
69*> HOWMNY is CHARACTER*1
70*> = 'A': compute condition numbers for all eigenpairs;
71*> = 'S': compute condition numbers for selected eigenpairs
72*> specified by the array SELECT.
73*> \endverbatim
74*>
75*> \param[in] SELECT
76*> \verbatim
77*> SELECT is LOGICAL array, dimension (N)
78*> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
79*> condition numbers are required. To select condition numbers
80*> for the eigenpair corresponding to a real eigenvalue w(j),
81*> SELECT(j) must be set to .TRUE.. To select condition numbers
82*> corresponding to a complex conjugate pair of eigenvalues w(j)
83*> and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
84*> set to .TRUE..
85*> If HOWMNY = 'A', SELECT is not referenced.
86*> \endverbatim
87*>
88*> \param[in] N
89*> \verbatim
90*> N is INTEGER
91*> The order of the square matrix pair (A, B). N >= 0.
92*> \endverbatim
93*>
94*> \param[in] A
95*> \verbatim
96*> A is DOUBLE PRECISION array, dimension (LDA,N)
97*> The upper quasi-triangular matrix A in the pair (A,B).
98*> \endverbatim
99*>
100*> \param[in] LDA
101*> \verbatim
102*> LDA is INTEGER
103*> The leading dimension of the array A. LDA >= max(1,N).
104*> \endverbatim
105*>
106*> \param[in] B
107*> \verbatim
108*> B is DOUBLE PRECISION array, dimension (LDB,N)
109*> The upper triangular matrix B in the pair (A,B).
110*> \endverbatim
111*>
112*> \param[in] LDB
113*> \verbatim
114*> LDB is INTEGER
115*> The leading dimension of the array B. LDB >= max(1,N).
116*> \endverbatim
117*>
118*> \param[in] VL
119*> \verbatim
120*> VL is DOUBLE PRECISION array, dimension (LDVL,M)
121*> If JOB = 'E' or 'B', VL must contain left eigenvectors of
122*> (A, B), corresponding to the eigenpairs specified by HOWMNY
123*> and SELECT. The eigenvectors must be stored in consecutive
124*> columns of VL, as returned by DTGEVC.
125*> If JOB = 'V', VL is not referenced.
126*> \endverbatim
127*>
128*> \param[in] LDVL
129*> \verbatim
130*> LDVL is INTEGER
131*> The leading dimension of the array VL. LDVL >= 1.
132*> If JOB = 'E' or 'B', LDVL >= N.
133*> \endverbatim
134*>
135*> \param[in] VR
136*> \verbatim
137*> VR is DOUBLE PRECISION array, dimension (LDVR,M)
138*> If JOB = 'E' or 'B', VR must contain right eigenvectors of
139*> (A, B), corresponding to the eigenpairs specified by HOWMNY
140*> and SELECT. The eigenvectors must be stored in consecutive
141*> columns ov VR, as returned by DTGEVC.
142*> If JOB = 'V', VR is not referenced.
143*> \endverbatim
144*>
145*> \param[in] LDVR
146*> \verbatim
147*> LDVR is INTEGER
148*> The leading dimension of the array VR. LDVR >= 1.
149*> If JOB = 'E' or 'B', LDVR >= N.
150*> \endverbatim
151*>
152*> \param[out] S
153*> \verbatim
154*> S is DOUBLE PRECISION array, dimension (MM)
155*> If JOB = 'E' or 'B', the reciprocal condition numbers of the
156*> selected eigenvalues, stored in consecutive elements of the
157*> array. For a complex conjugate pair of eigenvalues two
158*> consecutive elements of S are set to the same value. Thus
159*> S(j), DIF(j), and the j-th columns of VL and VR all
160*> correspond to the same eigenpair (but not in general the
161*> j-th eigenpair, unless all eigenpairs are selected).
162*> If JOB = 'V', S is not referenced.
163*> \endverbatim
164*>
165*> \param[out] DIF
166*> \verbatim
167*> DIF is DOUBLE PRECISION array, dimension (MM)
168*> If JOB = 'V' or 'B', the estimated reciprocal condition
169*> numbers of the selected eigenvectors, stored in consecutive
170*> elements of the array. For a complex eigenvector two
171*> consecutive elements of DIF are set to the same value. If
172*> the eigenvalues cannot be reordered to compute DIF(j), DIF(j)
173*> is set to 0; this can only occur when the true value would be
174*> very small anyway.
175*> If JOB = 'E', DIF is not referenced.
176*> \endverbatim
177*>
178*> \param[in] MM
179*> \verbatim
180*> MM is INTEGER
181*> The number of elements in the arrays S and DIF. MM >= M.
182*> \endverbatim
183*>
184*> \param[out] M
185*> \verbatim
186*> M is INTEGER
187*> The number of elements of the arrays S and DIF used to store
188*> the specified condition numbers; for each selected real
189*> eigenvalue one element is used, and for each selected complex
190*> conjugate pair of eigenvalues, two elements are used.
191*> If HOWMNY = 'A', M is set to N.
192*> \endverbatim
193*>
194*> \param[out] WORK
195*> \verbatim
196*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
197*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
198*> \endverbatim
199*>
200*> \param[in] LWORK
201*> \verbatim
202*> LWORK is INTEGER
203*> The dimension of the array WORK. LWORK >= max(1,N).
204*> If JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16.
205*>
206*> If LWORK = -1, then a workspace query is assumed; the routine
207*> only calculates the optimal size of the WORK array, returns
208*> this value as the first entry of the WORK array, and no error
209*> message related to LWORK is issued by XERBLA.
210*> \endverbatim
211*>
212*> \param[out] IWORK
213*> \verbatim
214*> IWORK is INTEGER array, dimension (N + 6)
215*> If JOB = 'E', IWORK is not referenced.
216*> \endverbatim
217*>
218*> \param[out] INFO
219*> \verbatim
220*> INFO is INTEGER
221*> =0: Successful exit
222*> <0: If INFO = -i, the i-th argument had an illegal value
223*> \endverbatim
224*
225* Authors:
226* ========
227*
228*> \author Univ. of Tennessee
229*> \author Univ. of California Berkeley
230*> \author Univ. of Colorado Denver
231*> \author NAG Ltd.
232*
233*> \ingroup tgsna
234*
235*> \par Further Details:
236* =====================
237*>
238*> \verbatim
239*>
240*> The reciprocal of the condition number of a generalized eigenvalue
241*> w = (a, b) is defined as
242*>
243*> S(w) = (|u**TAv|**2 + |u**TBv|**2)**(1/2) / (norm(u)*norm(v))
244*>
245*> where u and v are the left and right eigenvectors of (A, B)
246*> corresponding to w; |z| denotes the absolute value of the complex
247*> number, and norm(u) denotes the 2-norm of the vector u.
248*> The pair (a, b) corresponds to an eigenvalue w = a/b (= u**TAv/u**TBv)
249*> of the matrix pair (A, B). If both a and b equal zero, then (A B) is
250*> singular and S(I) = -1 is returned.
251*>
252*> An approximate error bound on the chordal distance between the i-th
253*> computed generalized eigenvalue w and the corresponding exact
254*> eigenvalue lambda is
255*>
256*> chord(w, lambda) <= EPS * norm(A, B) / S(I)
257*>
258*> where EPS is the machine precision.
259*>
260*> The reciprocal of the condition number DIF(i) of right eigenvector u
261*> and left eigenvector v corresponding to the generalized eigenvalue w
262*> is defined as follows:
263*>
264*> a) If the i-th eigenvalue w = (a,b) is real
265*>
266*> Suppose U and V are orthogonal transformations such that
267*>
268*> U**T*(A, B)*V = (S, T) = ( a * ) ( b * ) 1
269*> ( 0 S22 ),( 0 T22 ) n-1
270*> 1 n-1 1 n-1
271*>
272*> Then the reciprocal condition number DIF(i) is
273*>
274*> Difl((a, b), (S22, T22)) = sigma-min( Zl ),
275*>
276*> where sigma-min(Zl) denotes the smallest singular value of the
277*> 2(n-1)-by-2(n-1) matrix
278*>
279*> Zl = [ kron(a, In-1) -kron(1, S22) ]
280*> [ kron(b, In-1) -kron(1, T22) ] .
281*>
282*> Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
283*> Kronecker product between the matrices X and Y.
284*>
285*> Note that if the default method for computing DIF(i) is wanted
286*> (see DLATDF), then the parameter DIFDRI (see below) should be
287*> changed from 3 to 4 (routine DLATDF(IJOB = 2 will be used)).
288*> See DTGSYL for more details.
289*>
290*> b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
291*>
292*> Suppose U and V are orthogonal transformations such that
293*>
294*> U**T*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2
295*> ( 0 S22 ),( 0 T22) n-2
296*> 2 n-2 2 n-2
297*>
298*> and (S11, T11) corresponds to the complex conjugate eigenvalue
299*> pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
300*> that
301*>
302*> U1**T*S11*V1 = ( s11 s12 ) and U1**T*T11*V1 = ( t11 t12 )
303*> ( 0 s22 ) ( 0 t22 )
304*>
305*> where the generalized eigenvalues w = s11/t11 and
306*> conjg(w) = s22/t22.
307*>
308*> Then the reciprocal condition number DIF(i) is bounded by
309*>
310*> min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
311*>
312*> where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
313*> Z1 is the complex 2-by-2 matrix
314*>
315*> Z1 = [ s11 -s22 ]
316*> [ t11 -t22 ],
317*>
318*> This is done by computing (using real arithmetic) the
319*> roots of the characteristical polynomial det(Z1**T * Z1 - lambda I),
320*> where Z1**T denotes the transpose of Z1 and det(X) denotes
321*> the determinant of X.
322*>
323*> and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
324*> upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
325*>
326*> Z2 = [ kron(S11**T, In-2) -kron(I2, S22) ]
327*> [ kron(T11**T, In-2) -kron(I2, T22) ]
328*>
329*> Note that if the default method for computing DIF is wanted (see
330*> DLATDF), then the parameter DIFDRI (see below) should be changed
331*> from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). See DTGSYL
332*> for more details.
333*>
334*> For each eigenvalue/vector specified by SELECT, DIF stores a
335*> Frobenius norm-based estimate of Difl.
336*>
337*> An approximate error bound for the i-th computed eigenvector VL(i) or
338*> VR(i) is given by
339*>
340*> EPS * norm(A, B) / DIF(i).
341*>
342*> See ref. [2-3] for more details and further references.
343*> \endverbatim
344*
345*> \par Contributors:
346* ==================
347*>
348*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
349*> Umea University, S-901 87 Umea, Sweden.
350*
351*> \par References:
352* ================
353*>
354*> \verbatim
355*>
356*> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
357*> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
358*> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
359*> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
360*>
361*> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
362*> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
363*> Estimation: Theory, Algorithms and Software,
364*> Report UMINF - 94.04, Department of Computing Science, Umea
365*> University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
366*> Note 87. To appear in Numerical Algorithms, 1996.
367*>
368*> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
369*> for Solving the Generalized Sylvester Equation and Estimating the
370*> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
371*> Department of Computing Science, Umea University, S-901 87 Umea,
372*> Sweden, December 1993, Revised April 1994, Also as LAPACK Working
373*> Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
374*> No 1, 1996.
375*> \endverbatim
376*>
377* =====================================================================
378 SUBROUTINE dtgsna( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
379 $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
380 $ IWORK, INFO )
381*
382* -- LAPACK computational routine --
383* -- LAPACK is a software package provided by Univ. of Tennessee, --
384* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
385*
386* .. Scalar Arguments ..
387 CHARACTER HOWMNY, JOB
388 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
389* ..
390* .. Array Arguments ..
391 LOGICAL SELECT( * )
392 INTEGER IWORK( * )
393 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
394 $ vl( ldvl, * ), vr( ldvr, * ), work( * )
395* ..
396*
397* =====================================================================
398*
399* .. Parameters ..
400 INTEGER DIFDRI
401 PARAMETER ( DIFDRI = 3 )
402 DOUBLE PRECISION ZERO, ONE, TWO, FOUR
403 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
404 $ four = 4.0d+0 )
405* ..
406* .. Local Scalars ..
407 LOGICAL LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS
408 INTEGER I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2
409 DOUBLE PRECISION ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND,
410 $ eps, lnrm, rnrm, root1, root2, scale, smlnum,
411 $ tmpii, tmpir, tmpri, tmprr, uhav, uhavi, uhbv,
412 $ uhbvi
413* ..
414* .. Local Arrays ..
415 DOUBLE PRECISION DUMMY( 1 ), DUMMY1( 1 )
416* ..
417* .. External Functions ..
418 LOGICAL LSAME
419 DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
420 EXTERNAL lsame, ddot, dlamch, dlapy2, dnrm2
421* ..
422* .. External Subroutines ..
423 EXTERNAL dgemv, dlacpy, dlag2, dtgexc, dtgsyl, xerbla
424* ..
425* .. Intrinsic Functions ..
426 INTRINSIC max, min, sqrt
427* ..
428* .. Executable Statements ..
429*
430* Decode and test the input parameters
431*
432 wantbh = lsame( job, 'B' )
433 wants = lsame( job, 'E' ) .OR. wantbh
434 wantdf = lsame( job, 'V' ) .OR. wantbh
435*
436 somcon = lsame( howmny, 'S' )
437*
438 info = 0
439 lquery = ( lwork.EQ.-1 )
440*
441 IF( .NOT.wants .AND. .NOT.wantdf ) THEN
442 info = -1
443 ELSE IF( .NOT.lsame( howmny, 'A' ) .AND. .NOT.somcon ) THEN
444 info = -2
445 ELSE IF( n.LT.0 ) THEN
446 info = -4
447 ELSE IF( lda.LT.max( 1, n ) ) THEN
448 info = -6
449 ELSE IF( ldb.LT.max( 1, n ) ) THEN
450 info = -8
451 ELSE IF( wants .AND. ldvl.LT.n ) THEN
452 info = -10
453 ELSE IF( wants .AND. ldvr.LT.n ) THEN
454 info = -12
455 ELSE
456*
457* Set M to the number of eigenpairs for which condition numbers
458* are required, and test MM.
459*
460 IF( somcon ) THEN
461 m = 0
462 pair = .false.
463 DO 10 k = 1, n
464 IF( pair ) THEN
465 pair = .false.
466 ELSE
467 IF( k.LT.n ) THEN
468 IF( a( k+1, k ).EQ.zero ) THEN
469 IF( SELECT( k ) )
470 $ m = m + 1
471 ELSE
472 pair = .true.
473 IF( SELECT( k ) .OR. SELECT( k+1 ) )
474 $ m = m + 2
475 END IF
476 ELSE
477 IF( SELECT( n ) )
478 $ m = m + 1
479 END IF
480 END IF
481 10 CONTINUE
482 ELSE
483 m = n
484 END IF
485*
486 IF( n.EQ.0 ) THEN
487 lwmin = 1
488 ELSE IF( lsame( job, 'V' ) .OR. lsame( job, 'B' ) ) THEN
489 lwmin = 2*n*( n + 2 ) + 16
490 ELSE
491 lwmin = n
492 END IF
493 work( 1 ) = lwmin
494*
495 IF( mm.LT.m ) THEN
496 info = -15
497 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
498 info = -18
499 END IF
500 END IF
501*
502 IF( info.NE.0 ) THEN
503 CALL xerbla( 'DTGSNA', -info )
504 RETURN
505 ELSE IF( lquery ) THEN
506 RETURN
507 END IF
508*
509* Quick return if possible
510*
511 IF( n.EQ.0 )
512 $ RETURN
513*
514* Get machine constants
515*
516 eps = dlamch( 'P' )
517 smlnum = dlamch( 'S' ) / eps
518 ks = 0
519 pair = .false.
520*
521 DO 20 k = 1, n
522*
523* Determine whether A(k,k) begins a 1-by-1 or 2-by-2 block.
524*
525 IF( pair ) THEN
526 pair = .false.
527 GO TO 20
528 ELSE
529 IF( k.LT.n )
530 $ pair = a( k+1, k ).NE.zero
531 END IF
532*
533* Determine whether condition numbers are required for the k-th
534* eigenpair.
535*
536 IF( somcon ) THEN
537 IF( pair ) THEN
538 IF( .NOT.SELECT( k ) .AND. .NOT.SELECT( k+1 ) )
539 $ GO TO 20
540 ELSE
541 IF( .NOT.SELECT( k ) )
542 $ GO TO 20
543 END IF
544 END IF
545*
546 ks = ks + 1
547*
548 IF( wants ) THEN
549*
550* Compute the reciprocal condition number of the k-th
551* eigenvalue.
552*
553 IF( pair ) THEN
554*
555* Complex eigenvalue pair.
556*
557 rnrm = dlapy2( dnrm2( n, vr( 1, ks ), 1 ),
558 $ dnrm2( n, vr( 1, ks+1 ), 1 ) )
559 lnrm = dlapy2( dnrm2( n, vl( 1, ks ), 1 ),
560 $ dnrm2( n, vl( 1, ks+1 ), 1 ) )
561 CALL dgemv( 'N', n, n, one, a, lda, vr( 1, ks ), 1, zero,
562 $ work, 1 )
563 tmprr = ddot( n, work, 1, vl( 1, ks ), 1 )
564 tmpri = ddot( n, work, 1, vl( 1, ks+1 ), 1 )
565 CALL dgemv( 'N', n, n, one, a, lda, vr( 1, ks+1 ), 1,
566 $ zero, work, 1 )
567 tmpii = ddot( n, work, 1, vl( 1, ks+1 ), 1 )
568 tmpir = ddot( n, work, 1, vl( 1, ks ), 1 )
569 uhav = tmprr + tmpii
570 uhavi = tmpir - tmpri
571 CALL dgemv( 'N', n, n, one, b, ldb, vr( 1, ks ), 1, zero,
572 $ work, 1 )
573 tmprr = ddot( n, work, 1, vl( 1, ks ), 1 )
574 tmpri = ddot( n, work, 1, vl( 1, ks+1 ), 1 )
575 CALL dgemv( 'N', n, n, one, b, ldb, vr( 1, ks+1 ), 1,
576 $ zero, work, 1 )
577 tmpii = ddot( n, work, 1, vl( 1, ks+1 ), 1 )
578 tmpir = ddot( n, work, 1, vl( 1, ks ), 1 )
579 uhbv = tmprr + tmpii
580 uhbvi = tmpir - tmpri
581 uhav = dlapy2( uhav, uhavi )
582 uhbv = dlapy2( uhbv, uhbvi )
583 cond = dlapy2( uhav, uhbv )
584 s( ks ) = cond / ( rnrm*lnrm )
585 s( ks+1 ) = s( ks )
586*
587 ELSE
588*
589* Real eigenvalue.
590*
591 rnrm = dnrm2( n, vr( 1, ks ), 1 )
592 lnrm = dnrm2( n, vl( 1, ks ), 1 )
593 CALL dgemv( 'N', n, n, one, a, lda, vr( 1, ks ), 1, zero,
594 $ work, 1 )
595 uhav = ddot( n, work, 1, vl( 1, ks ), 1 )
596 CALL dgemv( 'N', n, n, one, b, ldb, vr( 1, ks ), 1, zero,
597 $ work, 1 )
598 uhbv = ddot( n, work, 1, vl( 1, ks ), 1 )
599 cond = dlapy2( uhav, uhbv )
600 IF( cond.EQ.zero ) THEN
601 s( ks ) = -one
602 ELSE
603 s( ks ) = cond / ( rnrm*lnrm )
604 END IF
605 END IF
606 END IF
607*
608 IF( wantdf ) THEN
609 IF( n.EQ.1 ) THEN
610 dif( ks ) = dlapy2( a( 1, 1 ), b( 1, 1 ) )
611 GO TO 20
612 END IF
613*
614* Estimate the reciprocal condition number of the k-th
615* eigenvectors.
616 IF( pair ) THEN
617*
618* Copy the 2-by 2 pencil beginning at (A(k,k), B(k, k)).
619* Compute the eigenvalue(s) at position K.
620*
621 work( 1 ) = a( k, k )
622 work( 2 ) = a( k+1, k )
623 work( 3 ) = a( k, k+1 )
624 work( 4 ) = a( k+1, k+1 )
625 work( 5 ) = b( k, k )
626 work( 6 ) = b( k+1, k )
627 work( 7 ) = b( k, k+1 )
628 work( 8 ) = b( k+1, k+1 )
629 CALL dlag2( work, 2, work( 5 ), 2, smlnum*eps, beta,
630 $ dummy1( 1 ), alphar, dummy( 1 ), alphai )
631 alprqt = one
632 c1 = two*( alphar*alphar+alphai*alphai+beta*beta )
633 c2 = four*beta*beta*alphai*alphai
634 root1 = c1 + sqrt( c1*c1-4.0d0*c2 )
635 root1 = root1 / two
636 root2 = c2 / root1
637 cond = min( sqrt( root1 ), sqrt( root2 ) )
638 END IF
639*
640* Copy the matrix (A, B) to the array WORK and swap the
641* diagonal block beginning at A(k,k) to the (1,1) position.
642*
643 CALL dlacpy( 'Full', n, n, a, lda, work, n )
644 CALL dlacpy( 'Full', n, n, b, ldb, work( n*n+1 ), n )
645 ifst = k
646 ilst = 1
647*
648 CALL dtgexc( .false., .false., n, work, n, work( n*n+1 ), n,
649 $ dummy, 1, dummy1, 1, ifst, ilst,
650 $ work( n*n*2+1 ), lwork-2*n*n, ierr )
651*
652 IF( ierr.GT.0 ) THEN
653*
654* Ill-conditioned problem - swap rejected.
655*
656 dif( ks ) = zero
657 ELSE
658*
659* Reordering successful, solve generalized Sylvester
660* equation for R and L,
661* A22 * R - L * A11 = A12
662* B22 * R - L * B11 = B12,
663* and compute estimate of Difl((A11,B11), (A22, B22)).
664*
665 n1 = 1
666 IF( work( 2 ).NE.zero )
667 $ n1 = 2
668 n2 = n - n1
669 IF( n2.EQ.0 ) THEN
670 dif( ks ) = cond
671 ELSE
672 i = n*n + 1
673 iz = 2*n*n + 1
674 CALL dtgsyl( 'N', difdri, n2, n1, work( n*n1+n1+1 ),
675 $ n, work, n, work( n1+1 ), n,
676 $ work( n*n1+n1+i ), n, work( i ), n,
677 $ work( n1+i ), n, scale, dif( ks ),
678 $ work( iz+1 ), lwork-2*n*n, iwork, ierr )
679*
680 IF( pair )
681 $ dif( ks ) = min( max( one, alprqt )*dif( ks ),
682 $ cond )
683 END IF
684 END IF
685 IF( pair )
686 $ dif( ks+1 ) = dif( ks )
687 END IF
688 IF( pair )
689 $ ks = ks + 1
690*
691 20 CONTINUE
692 work( 1 ) = lwmin
693 RETURN
694*
695* End of DTGSNA
696*
697 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
subroutine dlag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition dlag2.f:156
subroutine dtgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, work, lwork, info)
DTGEXC
Definition dtgexc.f:220
subroutine dtgsna(job, howmny, select, n, a, lda, b, ldb, vl, ldvl, vr, ldvr, s, dif, mm, m, work, lwork, iwork, info)
DTGSNA
Definition dtgsna.f:381
subroutine dtgsyl(trans, ijob, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, dif, work, lwork, iwork, info)
DTGSYL
Definition dtgsyl.f:299