Annotation of rpl/lapack/lapack/ztgsen.f, revision 1.21
1.10 bertrand 1: *> \brief \b ZTGSEN
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.17 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.10 bertrand 7: *
8: *> \htmlonly
1.17 bertrand 9: *> Download ZTGSEN + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsen.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsen.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsen.f">
1.10 bertrand 15: *> [TXT]</a>
1.17 bertrand 16: *> \endhtmlonly
1.10 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
22: * ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
23: * WORK, LWORK, IWORK, LIWORK, INFO )
1.17 bertrand 24: *
1.10 bertrand 25: * .. Scalar Arguments ..
26: * LOGICAL WANTQ, WANTZ
27: * INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
28: * $ M, N
29: * DOUBLE PRECISION PL, PR
30: * ..
31: * .. Array Arguments ..
32: * LOGICAL SELECT( * )
33: * INTEGER IWORK( * )
34: * DOUBLE PRECISION DIF( * )
35: * COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
36: * $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * )
37: * ..
1.17 bertrand 38: *
1.10 bertrand 39: *
40: *> \par Purpose:
41: * =============
42: *>
43: *> \verbatim
44: *>
45: *> ZTGSEN reorders the generalized Schur decomposition of a complex
46: *> matrix pair (A, B) (in terms of an unitary equivalence trans-
47: *> formation Q**H * (A, B) * Z), so that a selected cluster of eigenvalues
48: *> appears in the leading diagonal blocks of the pair (A,B). The leading
49: *> columns of Q and Z form unitary bases of the corresponding left and
50: *> right eigenspaces (deflating subspaces). (A, B) must be in
51: *> generalized Schur canonical form, that is, A and B are both upper
52: *> triangular.
53: *>
54: *> ZTGSEN also computes the generalized eigenvalues
55: *>
56: *> w(j)= ALPHA(j) / BETA(j)
57: *>
58: *> of the reordered matrix pair (A, B).
59: *>
60: *> Optionally, the routine computes estimates of reciprocal condition
61: *> numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
62: *> (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
63: *> between the matrix pairs (A11, B11) and (A22,B22) that correspond to
64: *> the selected cluster and the eigenvalues outside the cluster, resp.,
65: *> and norms of "projections" onto left and right eigenspaces w.r.t.
66: *> the selected cluster in the (1,1)-block.
67: *>
68: *> \endverbatim
69: *
70: * Arguments:
71: * ==========
72: *
73: *> \param[in] IJOB
74: *> \verbatim
1.19 bertrand 75: *> IJOB is INTEGER
1.10 bertrand 76: *> Specifies whether condition numbers are required for the
77: *> cluster of eigenvalues (PL and PR) or the deflating subspaces
78: *> (Difu and Difl):
79: *> =0: Only reorder w.r.t. SELECT. No extras.
80: *> =1: Reciprocal of norms of "projections" onto left and right
81: *> eigenspaces w.r.t. the selected cluster (PL and PR).
82: *> =2: Upper bounds on Difu and Difl. F-norm-based estimate
83: *> (DIF(1:2)).
84: *> =3: Estimate of Difu and Difl. 1-norm-based estimate
85: *> (DIF(1:2)).
86: *> About 5 times as expensive as IJOB = 2.
87: *> =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic
88: *> version to get it all.
89: *> =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above)
90: *> \endverbatim
91: *>
92: *> \param[in] WANTQ
93: *> \verbatim
94: *> WANTQ is LOGICAL
95: *> .TRUE. : update the left transformation matrix Q;
96: *> .FALSE.: do not update Q.
97: *> \endverbatim
98: *>
99: *> \param[in] WANTZ
100: *> \verbatim
101: *> WANTZ is LOGICAL
102: *> .TRUE. : update the right transformation matrix Z;
103: *> .FALSE.: do not update Z.
104: *> \endverbatim
105: *>
106: *> \param[in] SELECT
107: *> \verbatim
108: *> SELECT is LOGICAL array, dimension (N)
109: *> SELECT specifies the eigenvalues in the selected cluster. To
110: *> select an eigenvalue w(j), SELECT(j) must be set to
111: *> .TRUE..
112: *> \endverbatim
113: *>
114: *> \param[in] N
115: *> \verbatim
116: *> N is INTEGER
117: *> The order of the matrices A and B. N >= 0.
118: *> \endverbatim
119: *>
120: *> \param[in,out] A
121: *> \verbatim
122: *> A is COMPLEX*16 array, dimension(LDA,N)
123: *> On entry, the upper triangular matrix A, in generalized
124: *> Schur canonical form.
125: *> On exit, A is overwritten by the reordered matrix A.
126: *> \endverbatim
127: *>
128: *> \param[in] LDA
129: *> \verbatim
130: *> LDA is INTEGER
131: *> The leading dimension of the array A. LDA >= max(1,N).
132: *> \endverbatim
133: *>
134: *> \param[in,out] B
135: *> \verbatim
136: *> B is COMPLEX*16 array, dimension(LDB,N)
137: *> On entry, the upper triangular matrix B, in generalized
138: *> Schur canonical form.
139: *> On exit, B is overwritten by the reordered matrix B.
140: *> \endverbatim
141: *>
142: *> \param[in] LDB
143: *> \verbatim
144: *> LDB is INTEGER
145: *> The leading dimension of the array B. LDB >= max(1,N).
146: *> \endverbatim
147: *>
148: *> \param[out] ALPHA
149: *> \verbatim
150: *> ALPHA is COMPLEX*16 array, dimension (N)
151: *> \endverbatim
152: *>
153: *> \param[out] BETA
154: *> \verbatim
155: *> BETA is COMPLEX*16 array, dimension (N)
156: *>
157: *> The diagonal elements of A and B, respectively,
158: *> when the pair (A,B) has been reduced to generalized Schur
159: *> form. ALPHA(i)/BETA(i) i=1,...,N are the generalized
160: *> eigenvalues.
161: *> \endverbatim
162: *>
163: *> \param[in,out] Q
164: *> \verbatim
165: *> Q is COMPLEX*16 array, dimension (LDQ,N)
166: *> On entry, if WANTQ = .TRUE., Q is an N-by-N matrix.
167: *> On exit, Q has been postmultiplied by the left unitary
168: *> transformation matrix which reorder (A, B); The leading M
169: *> columns of Q form orthonormal bases for the specified pair of
170: *> left eigenspaces (deflating subspaces).
171: *> If WANTQ = .FALSE., Q is not referenced.
172: *> \endverbatim
173: *>
174: *> \param[in] LDQ
175: *> \verbatim
176: *> LDQ is INTEGER
177: *> The leading dimension of the array Q. LDQ >= 1.
178: *> If WANTQ = .TRUE., LDQ >= N.
179: *> \endverbatim
180: *>
181: *> \param[in,out] Z
182: *> \verbatim
183: *> Z is COMPLEX*16 array, dimension (LDZ,N)
184: *> On entry, if WANTZ = .TRUE., Z is an N-by-N matrix.
185: *> On exit, Z has been postmultiplied by the left unitary
186: *> transformation matrix which reorder (A, B); The leading M
187: *> columns of Z form orthonormal bases for the specified pair of
188: *> left eigenspaces (deflating subspaces).
189: *> If WANTZ = .FALSE., Z is not referenced.
190: *> \endverbatim
191: *>
192: *> \param[in] LDZ
193: *> \verbatim
194: *> LDZ is INTEGER
195: *> The leading dimension of the array Z. LDZ >= 1.
196: *> If WANTZ = .TRUE., LDZ >= N.
197: *> \endverbatim
198: *>
199: *> \param[out] M
200: *> \verbatim
201: *> M is INTEGER
202: *> The dimension of the specified pair of left and right
203: *> eigenspaces, (deflating subspaces) 0 <= M <= N.
204: *> \endverbatim
205: *>
206: *> \param[out] PL
207: *> \verbatim
208: *> PL is DOUBLE PRECISION
209: *> \endverbatim
210: *>
211: *> \param[out] PR
212: *> \verbatim
213: *> PR is DOUBLE PRECISION
214: *>
215: *> If IJOB = 1, 4 or 5, PL, PR are lower bounds on the
216: *> reciprocal of the norm of "projections" onto left and right
217: *> eigenspace with respect to the selected cluster.
218: *> 0 < PL, PR <= 1.
219: *> If M = 0 or M = N, PL = PR = 1.
220: *> If IJOB = 0, 2 or 3 PL, PR are not referenced.
221: *> \endverbatim
222: *>
223: *> \param[out] DIF
224: *> \verbatim
225: *> DIF is DOUBLE PRECISION array, dimension (2).
226: *> If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
227: *> If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
228: *> Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based
229: *> estimates of Difu and Difl, computed using reversed
230: *> communication with ZLACN2.
231: *> If M = 0 or N, DIF(1:2) = F-norm([A, B]).
232: *> If IJOB = 0 or 1, DIF is not referenced.
233: *> \endverbatim
234: *>
235: *> \param[out] WORK
236: *> \verbatim
237: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
238: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
239: *> \endverbatim
240: *>
241: *> \param[in] LWORK
242: *> \verbatim
243: *> LWORK is INTEGER
244: *> The dimension of the array WORK. LWORK >= 1
245: *> If IJOB = 1, 2 or 4, LWORK >= 2*M*(N-M)
246: *> If IJOB = 3 or 5, LWORK >= 4*M*(N-M)
247: *>
248: *> If LWORK = -1, then a workspace query is assumed; the routine
249: *> only calculates the optimal size of the WORK array, returns
250: *> this value as the first entry of the WORK array, and no error
251: *> message related to LWORK is issued by XERBLA.
252: *> \endverbatim
253: *>
254: *> \param[out] IWORK
255: *> \verbatim
256: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
257: *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
258: *> \endverbatim
259: *>
260: *> \param[in] LIWORK
261: *> \verbatim
262: *> LIWORK is INTEGER
263: *> The dimension of the array IWORK. LIWORK >= 1.
264: *> If IJOB = 1, 2 or 4, LIWORK >= N+2;
265: *> If IJOB = 3 or 5, LIWORK >= MAX(N+2, 2*M*(N-M));
266: *>
267: *> If LIWORK = -1, then a workspace query is assumed; the
268: *> routine only calculates the optimal size of the IWORK array,
269: *> returns this value as the first entry of the IWORK array, and
270: *> no error message related to LIWORK is issued by XERBLA.
271: *> \endverbatim
272: *>
273: *> \param[out] INFO
274: *> \verbatim
275: *> INFO is INTEGER
276: *> =0: Successful exit.
277: *> <0: If INFO = -i, the i-th argument had an illegal value.
278: *> =1: Reordering of (A, B) failed because the transformed
279: *> matrix pair (A, B) would be too far from generalized
280: *> Schur form; the problem is very ill-conditioned.
281: *> (A, B) may have been partially reordered.
282: *> If requested, 0 is returned in DIF(*), PL and PR.
283: *> \endverbatim
284: *
285: * Authors:
286: * ========
287: *
1.17 bertrand 288: *> \author Univ. of Tennessee
289: *> \author Univ. of California Berkeley
290: *> \author Univ. of Colorado Denver
291: *> \author NAG Ltd.
1.10 bertrand 292: *
293: *> \ingroup complex16OTHERcomputational
294: *
295: *> \par Further Details:
296: * =====================
297: *>
298: *> \verbatim
299: *>
300: *> ZTGSEN first collects the selected eigenvalues by computing unitary
301: *> U and W that move them to the top left corner of (A, B). In other
302: *> words, the selected eigenvalues are the eigenvalues of (A11, B11) in
303: *>
304: *> U**H*(A, B)*W = (A11 A12) (B11 B12) n1
305: *> ( 0 A22),( 0 B22) n2
306: *> n1 n2 n1 n2
307: *>
308: *> where N = n1+n2 and U**H means the conjugate transpose of U. The first
309: *> n1 columns of U and W span the specified pair of left and right
310: *> eigenspaces (deflating subspaces) of (A, B).
311: *>
312: *> If (A, B) has been obtained from the generalized real Schur
313: *> decomposition of a matrix pair (C, D) = Q*(A, B)*Z**H, then the
314: *> reordered generalized Schur form of (C, D) is given by
315: *>
316: *> (C, D) = (Q*U)*(U**H *(A, B)*W)*(Z*W)**H,
317: *>
318: *> and the first n1 columns of Q*U and Z*W span the corresponding
319: *> deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
320: *>
321: *> Note that if the selected eigenvalue is sufficiently ill-conditioned,
322: *> then its value may differ significantly from its value before
323: *> reordering.
324: *>
325: *> The reciprocal condition numbers of the left and right eigenspaces
326: *> spanned by the first n1 columns of U and W (or Q*U and Z*W) may
327: *> be returned in DIF(1:2), corresponding to Difu and Difl, resp.
328: *>
329: *> The Difu and Difl are defined as:
330: *>
331: *> Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
332: *> and
333: *> Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
334: *>
335: *> where sigma-min(Zu) is the smallest singular value of the
336: *> (2*n1*n2)-by-(2*n1*n2) matrix
337: *>
338: *> Zu = [ kron(In2, A11) -kron(A22**H, In1) ]
339: *> [ kron(In2, B11) -kron(B22**H, In1) ].
340: *>
341: *> Here, Inx is the identity matrix of size nx and A22**H is the
342: *> conjugate transpose of A22. kron(X, Y) is the Kronecker product between
343: *> the matrices X and Y.
344: *>
345: *> When DIF(2) is small, small changes in (A, B) can cause large changes
346: *> in the deflating subspace. An approximate (asymptotic) bound on the
347: *> maximum angular error in the computed deflating subspaces is
348: *>
349: *> EPS * norm((A, B)) / DIF(2),
350: *>
351: *> where EPS is the machine precision.
352: *>
353: *> The reciprocal norm of the projectors on the left and right
354: *> eigenspaces associated with (A11, B11) may be returned in PL and PR.
355: *> They are computed as follows. First we compute L and R so that
356: *> P*(A, B)*Q is block diagonal, where
357: *>
358: *> P = ( I -L ) n1 Q = ( I R ) n1
359: *> ( 0 I ) n2 and ( 0 I ) n2
360: *> n1 n2 n1 n2
361: *>
362: *> and (L, R) is the solution to the generalized Sylvester equation
363: *>
364: *> A11*R - L*A22 = -A12
365: *> B11*R - L*B22 = -B12
366: *>
367: *> Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
368: *> An approximate (asymptotic) bound on the average absolute error of
369: *> the selected eigenvalues is
370: *>
371: *> EPS * norm((A, B)) / PL.
372: *>
373: *> There are also global error bounds which valid for perturbations up
374: *> to a certain restriction: A lower bound (x) on the smallest
375: *> F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
376: *> coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
377: *> (i.e. (A + E, B + F), is
378: *>
379: *> x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
380: *>
381: *> An approximate bound on x can be computed from DIF(1:2), PL and PR.
382: *>
383: *> If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
384: *> (L', R') and unperturbed (L, R) left and right deflating subspaces
385: *> associated with the selected cluster in the (1,1)-blocks can be
386: *> bounded as
387: *>
388: *> max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
389: *> max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
390: *>
391: *> See LAPACK User's Guide section 4.11 or the following references
392: *> for more information.
393: *>
394: *> Note that if the default method for computing the Frobenius-norm-
395: *> based estimate DIF is not wanted (see ZLATDF), then the parameter
396: *> IDIFJB (see below) should be changed from 3 to 4 (routine ZLATDF
397: *> (IJOB = 2 will be used)). See ZTGSYL for more details.
398: *> \endverbatim
399: *
400: *> \par Contributors:
401: * ==================
402: *>
403: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
404: *> Umea University, S-901 87 Umea, Sweden.
405: *
406: *> \par References:
407: * ================
408: *>
409: *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
410: *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
411: *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
412: *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
413: *> \n
414: *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
415: *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
416: *> Estimation: Theory, Algorithms and Software, Report
417: *> UMINF - 94.04, Department of Computing Science, Umea University,
418: *> S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
419: *> To appear in Numerical Algorithms, 1996.
420: *> \n
421: *> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
422: *> for Solving the Generalized Sylvester Equation and Estimating the
423: *> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
424: *> Department of Computing Science, Umea University, S-901 87 Umea,
425: *> Sweden, December 1993, Revised April 1994, Also as LAPACK working
426: *> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
427: *> 1996.
428: *>
429: * =====================================================================
1.1 bertrand 430: SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
431: $ ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
432: $ WORK, LWORK, IWORK, LIWORK, INFO )
433: *
1.21 ! bertrand 434: * -- LAPACK computational routine --
1.1 bertrand 435: * -- LAPACK is a software package provided by Univ. of Tennessee, --
436: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
437: *
438: * .. Scalar Arguments ..
439: LOGICAL WANTQ, WANTZ
440: INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
441: $ M, N
442: DOUBLE PRECISION PL, PR
443: * ..
444: * .. Array Arguments ..
445: LOGICAL SELECT( * )
446: INTEGER IWORK( * )
447: DOUBLE PRECISION DIF( * )
448: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
449: $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * )
450: * ..
451: *
452: * =====================================================================
453: *
454: * .. Parameters ..
455: INTEGER IDIFJB
456: PARAMETER ( IDIFJB = 3 )
457: DOUBLE PRECISION ZERO, ONE
458: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
459: * ..
460: * .. Local Scalars ..
461: LOGICAL LQUERY, SWAP, WANTD, WANTD1, WANTD2, WANTP
462: INTEGER I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2,
463: $ N1, N2
464: DOUBLE PRECISION DSCALE, DSUM, RDSCAL, SAFMIN
465: COMPLEX*16 TEMP1, TEMP2
466: * ..
467: * .. Local Arrays ..
468: INTEGER ISAVE( 3 )
469: * ..
470: * .. External Subroutines ..
471: EXTERNAL XERBLA, ZLACN2, ZLACPY, ZLASSQ, ZSCAL, ZTGEXC,
472: $ ZTGSYL
473: * ..
474: * .. Intrinsic Functions ..
475: INTRINSIC ABS, DCMPLX, DCONJG, MAX, SQRT
476: * ..
477: * .. External Functions ..
478: DOUBLE PRECISION DLAMCH
479: EXTERNAL DLAMCH
480: * ..
481: * .. Executable Statements ..
482: *
483: * Decode and test the input parameters
484: *
485: INFO = 0
486: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
487: *
488: IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN
489: INFO = -1
490: ELSE IF( N.LT.0 ) THEN
491: INFO = -5
492: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
493: INFO = -7
494: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
495: INFO = -9
496: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
497: INFO = -13
498: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
499: INFO = -15
500: END IF
501: *
502: IF( INFO.NE.0 ) THEN
503: CALL XERBLA( 'ZTGSEN', -INFO )
504: RETURN
505: END IF
506: *
507: IERR = 0
508: *
509: WANTP = IJOB.EQ.1 .OR. IJOB.GE.4
510: WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4
511: WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5
512: WANTD = WANTD1 .OR. WANTD2
513: *
514: * Set M to the dimension of the specified pair of deflating
515: * subspaces.
516: *
517: M = 0
1.15 bertrand 518: IF( .NOT.LQUERY .OR. IJOB.NE.0 ) THEN
1.1 bertrand 519: DO 10 K = 1, N
520: ALPHA( K ) = A( K, K )
521: BETA( K ) = B( K, K )
522: IF( K.LT.N ) THEN
523: IF( SELECT( K ) )
524: $ M = M + 1
525: ELSE
526: IF( SELECT( N ) )
527: $ M = M + 1
528: END IF
529: 10 CONTINUE
1.15 bertrand 530: END IF
1.1 bertrand 531: *
532: IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
533: LWMIN = MAX( 1, 2*M*( N-M ) )
534: LIWMIN = MAX( 1, N+2 )
535: ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
536: LWMIN = MAX( 1, 4*M*( N-M ) )
537: LIWMIN = MAX( 1, 2*M*( N-M ), N+2 )
538: ELSE
539: LWMIN = 1
540: LIWMIN = 1
541: END IF
542: *
543: WORK( 1 ) = LWMIN
544: IWORK( 1 ) = LIWMIN
545: *
546: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
547: INFO = -21
548: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
549: INFO = -23
550: END IF
551: *
552: IF( INFO.NE.0 ) THEN
553: CALL XERBLA( 'ZTGSEN', -INFO )
554: RETURN
555: ELSE IF( LQUERY ) THEN
556: RETURN
557: END IF
558: *
559: * Quick return if possible.
560: *
561: IF( M.EQ.N .OR. M.EQ.0 ) THEN
562: IF( WANTP ) THEN
563: PL = ONE
564: PR = ONE
565: END IF
566: IF( WANTD ) THEN
567: DSCALE = ZERO
568: DSUM = ONE
569: DO 20 I = 1, N
570: CALL ZLASSQ( N, A( 1, I ), 1, DSCALE, DSUM )
571: CALL ZLASSQ( N, B( 1, I ), 1, DSCALE, DSUM )
572: 20 CONTINUE
573: DIF( 1 ) = DSCALE*SQRT( DSUM )
574: DIF( 2 ) = DIF( 1 )
575: END IF
576: GO TO 70
577: END IF
578: *
579: * Get machine constant
580: *
581: SAFMIN = DLAMCH( 'S' )
582: *
583: * Collect the selected blocks at the top-left corner of (A, B).
584: *
585: KS = 0
586: DO 30 K = 1, N
587: SWAP = SELECT( K )
588: IF( SWAP ) THEN
589: KS = KS + 1
590: *
591: * Swap the K-th block to position KS. Compute unitary Q
592: * and Z that will swap adjacent diagonal blocks in (A, B).
593: *
594: IF( K.NE.KS )
595: $ CALL ZTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
596: $ LDZ, K, KS, IERR )
597: *
598: IF( IERR.GT.0 ) THEN
599: *
600: * Swap is rejected: exit.
601: *
602: INFO = 1
603: IF( WANTP ) THEN
604: PL = ZERO
605: PR = ZERO
606: END IF
607: IF( WANTD ) THEN
608: DIF( 1 ) = ZERO
609: DIF( 2 ) = ZERO
610: END IF
611: GO TO 70
612: END IF
613: END IF
614: 30 CONTINUE
615: IF( WANTP ) THEN
616: *
617: * Solve generalized Sylvester equation for R and L:
618: * A11 * R - L * A22 = A12
619: * B11 * R - L * B22 = B12
620: *
621: N1 = M
622: N2 = N - M
623: I = N1 + 1
624: CALL ZLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
625: CALL ZLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
626: $ N1 )
627: IJB = 0
628: CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
629: $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1,
630: $ DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
631: $ LWORK-2*N1*N2, IWORK, IERR )
632: *
633: * Estimate the reciprocal of norms of "projections" onto
634: * left and right eigenspaces
635: *
636: RDSCAL = ZERO
637: DSUM = ONE
638: CALL ZLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM )
639: PL = RDSCAL*SQRT( DSUM )
640: IF( PL.EQ.ZERO ) THEN
641: PL = ONE
642: ELSE
643: PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) )
644: END IF
645: RDSCAL = ZERO
646: DSUM = ONE
647: CALL ZLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM )
648: PR = RDSCAL*SQRT( DSUM )
649: IF( PR.EQ.ZERO ) THEN
650: PR = ONE
651: ELSE
652: PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) )
653: END IF
654: END IF
655: IF( WANTD ) THEN
656: *
657: * Compute estimates Difu and Difl.
658: *
659: IF( WANTD1 ) THEN
660: N1 = M
661: N2 = N - M
662: I = N1 + 1
663: IJB = IDIFJB
664: *
665: * Frobenius norm-based Difu estimate.
666: *
667: CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
668: $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
669: $ N1, DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
670: $ LWORK-2*N1*N2, IWORK, IERR )
671: *
672: * Frobenius norm-based Difl estimate.
673: *
674: CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
675: $ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
676: $ N2, DSCALE, DIF( 2 ), WORK( N1*N2*2+1 ),
677: $ LWORK-2*N1*N2, IWORK, IERR )
678: ELSE
679: *
680: * Compute 1-norm-based estimates of Difu and Difl using
681: * reversed communication with ZLACN2. In each step a
682: * generalized Sylvester equation or a transposed variant
683: * is solved.
684: *
685: KASE = 0
686: N1 = M
687: N2 = N - M
688: I = N1 + 1
689: IJB = 0
690: MN2 = 2*N1*N2
691: *
692: * 1-norm-based estimate of Difu.
693: *
694: 40 CONTINUE
695: CALL ZLACN2( MN2, WORK( MN2+1 ), WORK, DIF( 1 ), KASE,
696: $ ISAVE )
697: IF( KASE.NE.0 ) THEN
698: IF( KASE.EQ.1 ) THEN
699: *
700: * Solve generalized Sylvester equation
701: *
702: CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
703: $ WORK, N1, B, LDB, B( I, I ), LDB,
704: $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
705: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
706: $ IERR )
707: ELSE
708: *
709: * Solve the transposed variant.
710: *
711: CALL ZTGSYL( 'C', IJB, N1, N2, A, LDA, A( I, I ), LDA,
712: $ WORK, N1, B, LDB, B( I, I ), LDB,
713: $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
714: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
715: $ IERR )
716: END IF
717: GO TO 40
718: END IF
719: DIF( 1 ) = DSCALE / DIF( 1 )
720: *
721: * 1-norm-based estimate of Difl.
722: *
723: 50 CONTINUE
724: CALL ZLACN2( MN2, WORK( MN2+1 ), WORK, DIF( 2 ), KASE,
725: $ ISAVE )
726: IF( KASE.NE.0 ) THEN
727: IF( KASE.EQ.1 ) THEN
728: *
729: * Solve generalized Sylvester equation
730: *
731: CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
732: $ WORK, N2, B( I, I ), LDB, B, LDB,
733: $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
734: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
735: $ IERR )
736: ELSE
737: *
738: * Solve the transposed variant.
739: *
740: CALL ZTGSYL( 'C', IJB, N2, N1, A( I, I ), LDA, A, LDA,
741: $ WORK, N2, B, LDB, B( I, I ), LDB,
742: $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
743: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
744: $ IERR )
745: END IF
746: GO TO 50
747: END IF
748: DIF( 2 ) = DSCALE / DIF( 2 )
749: END IF
750: END IF
751: *
752: * If B(K,K) is complex, make it real and positive (normalization
753: * of the generalized Schur form) and Store the generalized
754: * eigenvalues of reordered pair (A, B)
755: *
756: DO 60 K = 1, N
757: DSCALE = ABS( B( K, K ) )
758: IF( DSCALE.GT.SAFMIN ) THEN
759: TEMP1 = DCONJG( B( K, K ) / DSCALE )
760: TEMP2 = B( K, K ) / DSCALE
761: B( K, K ) = DSCALE
762: CALL ZSCAL( N-K, TEMP1, B( K, K+1 ), LDB )
763: CALL ZSCAL( N-K+1, TEMP1, A( K, K ), LDA )
764: IF( WANTQ )
765: $ CALL ZSCAL( N, TEMP2, Q( 1, K ), 1 )
766: ELSE
767: B( K, K ) = DCMPLX( ZERO, ZERO )
768: END IF
769: *
770: ALPHA( K ) = A( K, K )
771: BETA( K ) = B( K, K )
772: *
773: 60 CONTINUE
774: *
775: 70 CONTINUE
776: *
777: WORK( 1 ) = LWMIN
778: IWORK( 1 ) = LIWMIN
779: *
780: RETURN
781: *
782: * End of ZTGSEN
783: *
784: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztgsen.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack