Annotation of rpl/lapack/lapack/ztgsna.f, revision 1.11
1.9 bertrand 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: *> \date November 2011
217: *
218: *> \ingroup complex16OTHERcomputational
219: *
220: *> \par Further Details:
221: * =====================
222: *>
223: *> \verbatim
224: *>
225: *> The reciprocal of the condition number of the i-th generalized
226: *> eigenvalue w = (a, b) is defined as
227: *>
228: *> S(I) = (|v**HAu|**2 + |v**HBu|**2)**(1/2) / (norm(u)*norm(v))
229: *>
230: *> where u and v are the right and left eigenvectors of (A, B)
231: *> corresponding to w; |z| denotes the absolute value of the complex
232: *> number, and norm(u) denotes the 2-norm of the vector u. The pair
233: *> (a, b) corresponds to an eigenvalue w = a/b (= v**HAu/v**HBu) of the
234: *> matrix pair (A, B). If both a and b equal zero, then (A,B) is
235: *> singular and S(I) = -1 is returned.
236: *>
237: *> An approximate error bound on the chordal distance between the i-th
238: *> computed generalized eigenvalue w and the corresponding exact
239: *> eigenvalue lambda is
240: *>
241: *> chord(w, lambda) <= EPS * norm(A, B) / S(I),
242: *>
243: *> where EPS is the machine precision.
244: *>
245: *> The reciprocal of the condition number of the right eigenvector u
246: *> and left eigenvector v corresponding to the generalized eigenvalue w
247: *> is defined as follows. Suppose
248: *>
249: *> (A, B) = ( a * ) ( b * ) 1
250: *> ( 0 A22 ),( 0 B22 ) n-1
251: *> 1 n-1 1 n-1
252: *>
253: *> Then the reciprocal condition number DIF(I) is
254: *>
255: *> Difl[(a, b), (A22, B22)] = sigma-min( Zl )
256: *>
257: *> where sigma-min(Zl) denotes the smallest singular value of
258: *>
259: *> Zl = [ kron(a, In-1) -kron(1, A22) ]
260: *> [ kron(b, In-1) -kron(1, B22) ].
261: *>
262: *> Here In-1 is the identity matrix of size n-1 and X**H is the conjugate
263: *> transpose of X. kron(X, Y) is the Kronecker product between the
264: *> matrices X and Y.
265: *>
266: *> We approximate the smallest singular value of Zl with an upper
267: *> bound. This is done by ZLATDF.
268: *>
269: *> An approximate error bound for a computed eigenvector VL(i) or
270: *> VR(i) is given by
271: *>
272: *> EPS * norm(A, B) / DIF(i).
273: *>
274: *> See ref. [2-3] for more details and further references.
275: *> \endverbatim
276: *
277: *> \par Contributors:
278: * ==================
279: *>
280: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
281: *> Umea University, S-901 87 Umea, Sweden.
282: *
283: *> \par References:
284: * ================
285: *>
286: *> \verbatim
287: *>
288: *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
289: *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
290: *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
291: *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
292: *>
293: *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
294: *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
295: *> Estimation: Theory, Algorithms and Software, Report
296: *> UMINF - 94.04, Department of Computing Science, Umea University,
297: *> S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
298: *> To appear in Numerical Algorithms, 1996.
299: *>
300: *> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
301: *> for Solving the Generalized Sylvester Equation and Estimating the
302: *> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
303: *> Department of Computing Science, Umea University, S-901 87 Umea,
304: *> Sweden, December 1993, Revised April 1994, Also as LAPACK Working
305: *> Note 75.
306: *> To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
307: *> \endverbatim
308: *>
309: * =====================================================================
1.1 bertrand 310: SUBROUTINE ZTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
311: $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
312: $ IWORK, INFO )
313: *
1.9 bertrand 314: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 315: * -- LAPACK is a software package provided by Univ. of Tennessee, --
316: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 bertrand 317: * November 2011
1.1 bertrand 318: *
319: * .. Scalar Arguments ..
320: CHARACTER HOWMNY, JOB
321: INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
322: * ..
323: * .. Array Arguments ..
324: LOGICAL SELECT( * )
325: INTEGER IWORK( * )
326: DOUBLE PRECISION DIF( * ), S( * )
327: COMPLEX*16 A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
328: $ VR( LDVR, * ), WORK( * )
329: * ..
330: *
331: * =====================================================================
332: *
333: * .. Parameters ..
334: DOUBLE PRECISION ZERO, ONE
335: INTEGER IDIFJB
336: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, IDIFJB = 3 )
337: * ..
338: * .. Local Scalars ..
339: LOGICAL LQUERY, SOMCON, WANTBH, WANTDF, WANTS
340: INTEGER I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2
341: DOUBLE PRECISION BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM
342: COMPLEX*16 YHAX, YHBX
343: * ..
344: * .. Local Arrays ..
345: COMPLEX*16 DUMMY( 1 ), DUMMY1( 1 )
346: * ..
347: * .. External Functions ..
348: LOGICAL LSAME
349: DOUBLE PRECISION DLAMCH, DLAPY2, DZNRM2
350: COMPLEX*16 ZDOTC
351: EXTERNAL LSAME, DLAMCH, DLAPY2, DZNRM2, ZDOTC
352: * ..
353: * .. External Subroutines ..
354: EXTERNAL DLABAD, XERBLA, ZGEMV, ZLACPY, ZTGEXC, ZTGSYL
355: * ..
356: * .. Intrinsic Functions ..
357: INTRINSIC ABS, DCMPLX, MAX
358: * ..
359: * .. Executable Statements ..
360: *
361: * Decode and test the input parameters
362: *
363: WANTBH = LSAME( JOB, 'B' )
364: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
365: WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
366: *
367: SOMCON = LSAME( HOWMNY, 'S' )
368: *
369: INFO = 0
370: LQUERY = ( LWORK.EQ.-1 )
371: *
372: IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
373: INFO = -1
374: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
375: INFO = -2
376: ELSE IF( N.LT.0 ) THEN
377: INFO = -4
378: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
379: INFO = -6
380: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
381: INFO = -8
382: ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
383: INFO = -10
384: ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
385: INFO = -12
386: ELSE
387: *
388: * Set M to the number of eigenpairs for which condition numbers
389: * are required, and test MM.
390: *
391: IF( SOMCON ) THEN
392: M = 0
393: DO 10 K = 1, N
394: IF( SELECT( K ) )
395: $ M = M + 1
396: 10 CONTINUE
397: ELSE
398: M = N
399: END IF
400: *
401: IF( N.EQ.0 ) THEN
402: LWMIN = 1
403: ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
404: LWMIN = 2*N*N
405: ELSE
406: LWMIN = N
407: END IF
408: WORK( 1 ) = LWMIN
409: *
410: IF( MM.LT.M ) THEN
411: INFO = -15
412: ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
413: INFO = -18
414: END IF
415: END IF
416: *
417: IF( INFO.NE.0 ) THEN
418: CALL XERBLA( 'ZTGSNA', -INFO )
419: RETURN
420: ELSE IF( LQUERY ) THEN
421: RETURN
422: END IF
423: *
424: * Quick return if possible
425: *
426: IF( N.EQ.0 )
427: $ RETURN
428: *
429: * Get machine constants
430: *
431: EPS = DLAMCH( 'P' )
432: SMLNUM = DLAMCH( 'S' ) / EPS
433: BIGNUM = ONE / SMLNUM
434: CALL DLABAD( SMLNUM, BIGNUM )
435: KS = 0
436: DO 20 K = 1, N
437: *
438: * Determine whether condition numbers are required for the k-th
439: * eigenpair.
440: *
441: IF( SOMCON ) THEN
442: IF( .NOT.SELECT( K ) )
443: $ GO TO 20
444: END IF
445: *
446: KS = KS + 1
447: *
448: IF( WANTS ) THEN
449: *
450: * Compute the reciprocal condition number of the k-th
451: * eigenvalue.
452: *
453: RNRM = DZNRM2( N, VR( 1, KS ), 1 )
454: LNRM = DZNRM2( N, VL( 1, KS ), 1 )
455: CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), A, LDA,
456: $ VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
457: YHAX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
458: CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), B, LDB,
459: $ VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
460: YHBX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
461: COND = DLAPY2( ABS( YHAX ), ABS( YHBX ) )
462: IF( COND.EQ.ZERO ) THEN
463: S( KS ) = -ONE
464: ELSE
465: S( KS ) = COND / ( RNRM*LNRM )
466: END IF
467: END IF
468: *
469: IF( WANTDF ) THEN
470: IF( N.EQ.1 ) THEN
471: DIF( KS ) = DLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) )
472: ELSE
473: *
474: * Estimate the reciprocal condition number of the k-th
475: * eigenvectors.
476: *
477: * Copy the matrix (A, B) to the array WORK and move the
478: * (k,k)th pair to the (1,1) position.
479: *
480: CALL ZLACPY( 'Full', N, N, A, LDA, WORK, N )
481: CALL ZLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
482: IFST = K
483: ILST = 1
484: *
485: CALL ZTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ),
486: $ N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR )
487: *
488: IF( IERR.GT.0 ) THEN
489: *
490: * Ill-conditioned problem - swap rejected.
491: *
492: DIF( KS ) = ZERO
493: ELSE
494: *
495: * Reordering successful, solve generalized Sylvester
496: * equation for R and L,
497: * A22 * R - L * A11 = A12
498: * B22 * R - L * B11 = B12,
499: * and compute estimate of Difl[(A11,B11), (A22, B22)].
500: *
501: N1 = 1
502: N2 = N - N1
503: I = N*N + 1
504: CALL ZTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ),
505: $ N, WORK, N, WORK( N1+1 ), N,
506: $ WORK( N*N1+N1+I ), N, WORK( I ), N,
507: $ WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY,
508: $ 1, IWORK, IERR )
509: END IF
510: END IF
511: END IF
512: *
513: 20 CONTINUE
514: WORK( 1 ) = LWMIN
515: RETURN
516: *
517: * End of ZTGSNA
518: *
519: 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