1: *> \brief \b ZTGSNA
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZTGSNA + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsna.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsna.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsna.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTGSNA( 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 DIF( * ), S( * )
33: * COMPLEX*16 A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
34: * $ VR( LDVR, * ), WORK( * )
35: * ..
36: *
37: *
38: *> \par Purpose:
39: * =============
40: *>
41: *> \verbatim
42: *>
43: *> ZTGSNA estimates reciprocal condition numbers for specified
44: *> eigenvalues and/or eigenvectors of a matrix pair (A, B).
45: *>
46: *> (A, B) must be in generalized Schur canonical form, that is, A and
47: *> B are both upper triangular.
48: *> \endverbatim
49: *
50: * Arguments:
51: * ==========
52: *
53: *> \param[in] JOB
54: *> \verbatim
55: *> JOB is CHARACTER*1
56: *> Specifies whether condition numbers are required for
57: *> eigenvalues (S) or eigenvectors (DIF):
58: *> = 'E': for eigenvalues only (S);
59: *> = 'V': for eigenvectors only (DIF);
60: *> = 'B': for both eigenvalues and eigenvectors (S and DIF).
61: *> \endverbatim
62: *>
63: *> \param[in] HOWMNY
64: *> \verbatim
65: *> HOWMNY is CHARACTER*1
66: *> = 'A': compute condition numbers for all eigenpairs;
67: *> = 'S': compute condition numbers for selected eigenpairs
68: *> specified by the array SELECT.
69: *> \endverbatim
70: *>
71: *> \param[in] SELECT
72: *> \verbatim
73: *> SELECT is LOGICAL array, dimension (N)
74: *> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
75: *> condition numbers are required. To select condition numbers
76: *> for the corresponding j-th eigenvalue and/or eigenvector,
77: *> SELECT(j) must be set to .TRUE..
78: *> If HOWMNY = 'A', SELECT is not referenced.
79: *> \endverbatim
80: *>
81: *> \param[in] N
82: *> \verbatim
83: *> N is INTEGER
84: *> The order of the square matrix pair (A, B). N >= 0.
85: *> \endverbatim
86: *>
87: *> \param[in] A
88: *> \verbatim
89: *> A is COMPLEX*16 array, dimension (LDA,N)
90: *> The upper triangular matrix A in the pair (A,B).
91: *> \endverbatim
92: *>
93: *> \param[in] LDA
94: *> \verbatim
95: *> LDA is INTEGER
96: *> The leading dimension of the array A. LDA >= max(1,N).
97: *> \endverbatim
98: *>
99: *> \param[in] B
100: *> \verbatim
101: *> B is COMPLEX*16 array, dimension (LDB,N)
102: *> The upper triangular matrix B in the pair (A, B).
103: *> \endverbatim
104: *>
105: *> \param[in] LDB
106: *> \verbatim
107: *> LDB is INTEGER
108: *> The leading dimension of the array B. LDB >= max(1,N).
109: *> \endverbatim
110: *>
111: *> \param[in] VL
112: *> \verbatim
113: *> VL is COMPLEX*16 array, dimension (LDVL,M)
114: *> IF JOB = 'E' or 'B', VL must contain left eigenvectors of
115: *> (A, B), corresponding to the eigenpairs specified by HOWMNY
116: *> and SELECT. The eigenvectors must be stored in consecutive
117: *> columns of VL, as returned by ZTGEVC.
118: *> If JOB = 'V', VL is not referenced.
119: *> \endverbatim
120: *>
121: *> \param[in] LDVL
122: *> \verbatim
123: *> LDVL is INTEGER
124: *> The leading dimension of the array VL. LDVL >= 1; and
125: *> If JOB = 'E' or 'B', LDVL >= N.
126: *> \endverbatim
127: *>
128: *> \param[in] VR
129: *> \verbatim
130: *> VR is COMPLEX*16 array, dimension (LDVR,M)
131: *> IF JOB = 'E' or 'B', VR must contain right eigenvectors of
132: *> (A, B), corresponding to the eigenpairs specified by HOWMNY
133: *> and SELECT. The eigenvectors must be stored in consecutive
134: *> columns of VR, as returned by ZTGEVC.
135: *> If JOB = 'V', VR is not referenced.
136: *> \endverbatim
137: *>
138: *> \param[in] LDVR
139: *> \verbatim
140: *> LDVR is INTEGER
141: *> The leading dimension of the array VR. LDVR >= 1;
142: *> If JOB = 'E' or 'B', LDVR >= N.
143: *> \endverbatim
144: *>
145: *> \param[out] S
146: *> \verbatim
147: *> S is DOUBLE PRECISION array, dimension (MM)
148: *> If JOB = 'E' or 'B', the reciprocal condition numbers of the
149: *> selected eigenvalues, stored in consecutive elements of the
150: *> array.
151: *> If JOB = 'V', S is not referenced.
152: *> \endverbatim
153: *>
154: *> \param[out] DIF
155: *> \verbatim
156: *> DIF 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.
160: *> If the eigenvalues cannot be reordered to compute DIF(j),
161: *> DIF(j) is set to 0; this can only occur when the true value
162: *> would be very small anyway.
163: *> For each eigenvalue/vector specified by SELECT, DIF stores
164: *> a Frobenius norm-based estimate of Difl.
165: *> If JOB = 'E', DIF is not referenced.
166: *> \endverbatim
167: *>
168: *> \param[in] MM
169: *> \verbatim
170: *> MM is INTEGER
171: *> The number of elements in the arrays S and DIF. 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 DIF used to store
178: *> the specified condition numbers; for each selected eigenvalue
179: *> one element is used. If HOWMNY = 'A', M is set to N.
180: *> \endverbatim
181: *>
182: *> \param[out] WORK
183: *> \verbatim
184: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
185: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
186: *> \endverbatim
187: *>
188: *> \param[in] LWORK
189: *> \verbatim
190: *> LWORK is INTEGER
191: *> The dimension of the array WORK. LWORK >= max(1,N).
192: *> If JOB = 'V' or 'B', LWORK >= max(1,2*N*N).
193: *> \endverbatim
194: *>
195: *> \param[out] IWORK
196: *> \verbatim
197: *> IWORK is INTEGER array, dimension (N+2)
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: *> \ingroup complex16OTHERcomputational
217: *
218: *> \par Further Details:
219: * =====================
220: *>
221: *> \verbatim
222: *>
223: *> The reciprocal of the condition number of the i-th generalized
224: *> eigenvalue w = (a, b) is defined as
225: *>
226: *> S(I) = (|v**HAu|**2 + |v**HBu|**2)**(1/2) / (norm(u)*norm(v))
227: *>
228: *> where u and v are the right and left eigenvectors of (A, B)
229: *> corresponding to w; |z| denotes the absolute value of the complex
230: *> number, and norm(u) denotes the 2-norm of the vector u. The pair
231: *> (a, b) corresponds to an eigenvalue w = a/b (= v**HAu/v**HBu) of the
232: *> matrix pair (A, B). If both a and b equal zero, then (A,B) is
233: *> singular and S(I) = -1 is returned.
234: *>
235: *> An approximate error bound on the chordal distance between the i-th
236: *> computed generalized eigenvalue w and the corresponding exact
237: *> eigenvalue lambda is
238: *>
239: *> chord(w, lambda) <= EPS * norm(A, B) / S(I),
240: *>
241: *> where EPS is the machine precision.
242: *>
243: *> The reciprocal of the condition number of the right eigenvector u
244: *> and left eigenvector v corresponding to the generalized eigenvalue w
245: *> is defined as follows. Suppose
246: *>
247: *> (A, B) = ( a * ) ( b * ) 1
248: *> ( 0 A22 ),( 0 B22 ) n-1
249: *> 1 n-1 1 n-1
250: *>
251: *> Then the reciprocal condition number DIF(I) is
252: *>
253: *> Difl[(a, b), (A22, B22)] = sigma-min( Zl )
254: *>
255: *> where sigma-min(Zl) denotes the smallest singular value of
256: *>
257: *> Zl = [ kron(a, In-1) -kron(1, A22) ]
258: *> [ kron(b, In-1) -kron(1, B22) ].
259: *>
260: *> Here In-1 is the identity matrix of size n-1 and X**H is the conjugate
261: *> transpose of X. kron(X, Y) is the Kronecker product between the
262: *> matrices X and Y.
263: *>
264: *> We approximate the smallest singular value of Zl with an upper
265: *> bound. This is done by ZLATDF.
266: *>
267: *> An approximate error bound for a computed eigenvector VL(i) or
268: *> VR(i) is given by
269: *>
270: *> EPS * norm(A, B) / DIF(i).
271: *>
272: *> See ref. [2-3] for more details and further references.
273: *> \endverbatim
274: *
275: *> \par Contributors:
276: * ==================
277: *>
278: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
279: *> Umea University, S-901 87 Umea, Sweden.
280: *
281: *> \par References:
282: * ================
283: *>
284: *> \verbatim
285: *>
286: *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
287: *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
288: *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
289: *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
290: *>
291: *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
292: *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
293: *> Estimation: Theory, Algorithms and Software, Report
294: *> UMINF - 94.04, Department of Computing Science, Umea University,
295: *> S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
296: *> To appear in Numerical Algorithms, 1996.
297: *>
298: *> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
299: *> for Solving the Generalized Sylvester Equation and Estimating the
300: *> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
301: *> Department of Computing Science, Umea University, S-901 87 Umea,
302: *> Sweden, December 1993, Revised April 1994, Also as LAPACK Working
303: *> Note 75.
304: *> To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
305: *> \endverbatim
306: *>
307: * =====================================================================
308: SUBROUTINE ZTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
309: $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
310: $ IWORK, INFO )
311: *
312: * -- LAPACK computational routine --
313: * -- LAPACK is a software package provided by Univ. of Tennessee, --
314: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
315: *
316: * .. Scalar Arguments ..
317: CHARACTER HOWMNY, JOB
318: INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
319: * ..
320: * .. Array Arguments ..
321: LOGICAL SELECT( * )
322: INTEGER IWORK( * )
323: DOUBLE PRECISION DIF( * ), S( * )
324: COMPLEX*16 A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
325: $ VR( LDVR, * ), WORK( * )
326: * ..
327: *
328: * =====================================================================
329: *
330: * .. Parameters ..
331: DOUBLE PRECISION ZERO, ONE
332: INTEGER IDIFJB
333: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, IDIFJB = 3 )
334: * ..
335: * .. Local Scalars ..
336: LOGICAL LQUERY, SOMCON, WANTBH, WANTDF, WANTS
337: INTEGER I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2
338: DOUBLE PRECISION BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM
339: COMPLEX*16 YHAX, YHBX
340: * ..
341: * .. Local Arrays ..
342: COMPLEX*16 DUMMY( 1 ), DUMMY1( 1 )
343: * ..
344: * .. External Functions ..
345: LOGICAL LSAME
346: DOUBLE PRECISION DLAMCH, DLAPY2, DZNRM2
347: COMPLEX*16 ZDOTC
348: EXTERNAL LSAME, DLAMCH, DLAPY2, DZNRM2, ZDOTC
349: * ..
350: * .. External Subroutines ..
351: EXTERNAL DLABAD, XERBLA, ZGEMV, ZLACPY, ZTGEXC, ZTGSYL
352: * ..
353: * .. Intrinsic Functions ..
354: INTRINSIC ABS, DCMPLX, MAX
355: * ..
356: * .. Executable Statements ..
357: *
358: * Decode and test the input parameters
359: *
360: WANTBH = LSAME( JOB, 'B' )
361: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
362: WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
363: *
364: SOMCON = LSAME( HOWMNY, 'S' )
365: *
366: INFO = 0
367: LQUERY = ( LWORK.EQ.-1 )
368: *
369: IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
370: INFO = -1
371: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
372: INFO = -2
373: ELSE IF( N.LT.0 ) THEN
374: INFO = -4
375: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
376: INFO = -6
377: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
378: INFO = -8
379: ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
380: INFO = -10
381: ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
382: INFO = -12
383: ELSE
384: *
385: * Set M to the number of eigenpairs for which condition numbers
386: * are required, and test MM.
387: *
388: IF( SOMCON ) THEN
389: M = 0
390: DO 10 K = 1, N
391: IF( SELECT( K ) )
392: $ M = M + 1
393: 10 CONTINUE
394: ELSE
395: M = N
396: END IF
397: *
398: IF( N.EQ.0 ) THEN
399: LWMIN = 1
400: ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
401: LWMIN = 2*N*N
402: ELSE
403: LWMIN = N
404: END IF
405: WORK( 1 ) = LWMIN
406: *
407: IF( MM.LT.M ) THEN
408: INFO = -15
409: ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
410: INFO = -18
411: END IF
412: END IF
413: *
414: IF( INFO.NE.0 ) THEN
415: CALL XERBLA( 'ZTGSNA', -INFO )
416: RETURN
417: ELSE IF( LQUERY ) THEN
418: RETURN
419: END IF
420: *
421: * Quick return if possible
422: *
423: IF( N.EQ.0 )
424: $ RETURN
425: *
426: * Get machine constants
427: *
428: EPS = DLAMCH( 'P' )
429: SMLNUM = DLAMCH( 'S' ) / EPS
430: BIGNUM = ONE / SMLNUM
431: CALL DLABAD( SMLNUM, BIGNUM )
432: KS = 0
433: DO 20 K = 1, N
434: *
435: * Determine whether condition numbers are required for the k-th
436: * eigenpair.
437: *
438: IF( SOMCON ) THEN
439: IF( .NOT.SELECT( K ) )
440: $ GO TO 20
441: END IF
442: *
443: KS = KS + 1
444: *
445: IF( WANTS ) THEN
446: *
447: * Compute the reciprocal condition number of the k-th
448: * eigenvalue.
449: *
450: RNRM = DZNRM2( N, VR( 1, KS ), 1 )
451: LNRM = DZNRM2( N, VL( 1, KS ), 1 )
452: CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), A, LDA,
453: $ VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
454: YHAX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
455: CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), B, LDB,
456: $ VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
457: YHBX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
458: COND = DLAPY2( ABS( YHAX ), ABS( YHBX ) )
459: IF( COND.EQ.ZERO ) THEN
460: S( KS ) = -ONE
461: ELSE
462: S( KS ) = COND / ( RNRM*LNRM )
463: END IF
464: END IF
465: *
466: IF( WANTDF ) THEN
467: IF( N.EQ.1 ) THEN
468: DIF( KS ) = DLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) )
469: ELSE
470: *
471: * Estimate the reciprocal condition number of the k-th
472: * eigenvectors.
473: *
474: * Copy the matrix (A, B) to the array WORK and move the
475: * (k,k)th pair to the (1,1) position.
476: *
477: CALL ZLACPY( 'Full', N, N, A, LDA, WORK, N )
478: CALL ZLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
479: IFST = K
480: ILST = 1
481: *
482: CALL ZTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ),
483: $ N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR )
484: *
485: IF( IERR.GT.0 ) THEN
486: *
487: * Ill-conditioned problem - swap rejected.
488: *
489: DIF( KS ) = ZERO
490: ELSE
491: *
492: * Reordering successful, solve generalized Sylvester
493: * equation for R and L,
494: * A22 * R - L * A11 = A12
495: * B22 * R - L * B11 = B12,
496: * and compute estimate of Difl[(A11,B11), (A22, B22)].
497: *
498: N1 = 1
499: N2 = N - N1
500: I = N*N + 1
501: CALL ZTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ),
502: $ N, WORK, N, WORK( N1+1 ), N,
503: $ WORK( N*N1+N1+I ), N, WORK( I ), N,
504: $ WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY,
505: $ 1, IWORK, IERR )
506: END IF
507: END IF
508: END IF
509: *
510: 20 CONTINUE
511: WORK( 1 ) = LWMIN
512: RETURN
513: *
514: * End of ZTGSNA
515: *
516: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztgsna.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack