1: *> \brief \b ZTGSEN
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
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">
15: *> [TXT]</a>
16: *> \endhtmlonly
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 )
24: *
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: * ..
38: *
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
75: *> IJOB is integer
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: *
288: *> \author Univ. of Tennessee
289: *> \author Univ. of California Berkeley
290: *> \author Univ. of Colorado Denver
291: *> \author NAG Ltd.
292: *
293: *> \date June 2016
294: *
295: *> \ingroup complex16OTHERcomputational
296: *
297: *> \par Further Details:
298: * =====================
299: *>
300: *> \verbatim
301: *>
302: *> ZTGSEN first collects the selected eigenvalues by computing unitary
303: *> U and W that move them to the top left corner of (A, B). In other
304: *> words, the selected eigenvalues are the eigenvalues of (A11, B11) in
305: *>
306: *> U**H*(A, B)*W = (A11 A12) (B11 B12) n1
307: *> ( 0 A22),( 0 B22) n2
308: *> n1 n2 n1 n2
309: *>
310: *> where N = n1+n2 and U**H means the conjugate transpose of U. The first
311: *> n1 columns of U and W span the specified pair of left and right
312: *> eigenspaces (deflating subspaces) of (A, B).
313: *>
314: *> If (A, B) has been obtained from the generalized real Schur
315: *> decomposition of a matrix pair (C, D) = Q*(A, B)*Z**H, then the
316: *> reordered generalized Schur form of (C, D) is given by
317: *>
318: *> (C, D) = (Q*U)*(U**H *(A, B)*W)*(Z*W)**H,
319: *>
320: *> and the first n1 columns of Q*U and Z*W span the corresponding
321: *> deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
322: *>
323: *> Note that if the selected eigenvalue is sufficiently ill-conditioned,
324: *> then its value may differ significantly from its value before
325: *> reordering.
326: *>
327: *> The reciprocal condition numbers of the left and right eigenspaces
328: *> spanned by the first n1 columns of U and W (or Q*U and Z*W) may
329: *> be returned in DIF(1:2), corresponding to Difu and Difl, resp.
330: *>
331: *> The Difu and Difl are defined as:
332: *>
333: *> Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
334: *> and
335: *> Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
336: *>
337: *> where sigma-min(Zu) is the smallest singular value of the
338: *> (2*n1*n2)-by-(2*n1*n2) matrix
339: *>
340: *> Zu = [ kron(In2, A11) -kron(A22**H, In1) ]
341: *> [ kron(In2, B11) -kron(B22**H, In1) ].
342: *>
343: *> Here, Inx is the identity matrix of size nx and A22**H is the
344: *> conjugate transpose of A22. kron(X, Y) is the Kronecker product between
345: *> the matrices X and Y.
346: *>
347: *> When DIF(2) is small, small changes in (A, B) can cause large changes
348: *> in the deflating subspace. An approximate (asymptotic) bound on the
349: *> maximum angular error in the computed deflating subspaces is
350: *>
351: *> EPS * norm((A, B)) / DIF(2),
352: *>
353: *> where EPS is the machine precision.
354: *>
355: *> The reciprocal norm of the projectors on the left and right
356: *> eigenspaces associated with (A11, B11) may be returned in PL and PR.
357: *> They are computed as follows. First we compute L and R so that
358: *> P*(A, B)*Q is block diagonal, where
359: *>
360: *> P = ( I -L ) n1 Q = ( I R ) n1
361: *> ( 0 I ) n2 and ( 0 I ) n2
362: *> n1 n2 n1 n2
363: *>
364: *> and (L, R) is the solution to the generalized Sylvester equation
365: *>
366: *> A11*R - L*A22 = -A12
367: *> B11*R - L*B22 = -B12
368: *>
369: *> Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
370: *> An approximate (asymptotic) bound on the average absolute error of
371: *> the selected eigenvalues is
372: *>
373: *> EPS * norm((A, B)) / PL.
374: *>
375: *> There are also global error bounds which valid for perturbations up
376: *> to a certain restriction: A lower bound (x) on the smallest
377: *> F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
378: *> coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
379: *> (i.e. (A + E, B + F), is
380: *>
381: *> x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
382: *>
383: *> An approximate bound on x can be computed from DIF(1:2), PL and PR.
384: *>
385: *> If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
386: *> (L', R') and unperturbed (L, R) left and right deflating subspaces
387: *> associated with the selected cluster in the (1,1)-blocks can be
388: *> bounded as
389: *>
390: *> max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
391: *> max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
392: *>
393: *> See LAPACK User's Guide section 4.11 or the following references
394: *> for more information.
395: *>
396: *> Note that if the default method for computing the Frobenius-norm-
397: *> based estimate DIF is not wanted (see ZLATDF), then the parameter
398: *> IDIFJB (see below) should be changed from 3 to 4 (routine ZLATDF
399: *> (IJOB = 2 will be used)). See ZTGSYL for more details.
400: *> \endverbatim
401: *
402: *> \par Contributors:
403: * ==================
404: *>
405: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
406: *> Umea University, S-901 87 Umea, Sweden.
407: *
408: *> \par References:
409: * ================
410: *>
411: *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
412: *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
413: *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
414: *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
415: *> \n
416: *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
417: *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
418: *> Estimation: Theory, Algorithms and Software, Report
419: *> UMINF - 94.04, Department of Computing Science, Umea University,
420: *> S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
421: *> To appear in Numerical Algorithms, 1996.
422: *> \n
423: *> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
424: *> for Solving the Generalized Sylvester Equation and Estimating the
425: *> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
426: *> Department of Computing Science, Umea University, S-901 87 Umea,
427: *> Sweden, December 1993, Revised April 1994, Also as LAPACK working
428: *> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
429: *> 1996.
430: *>
431: * =====================================================================
432: SUBROUTINE ZTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
433: $ ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF,
434: $ WORK, LWORK, IWORK, LIWORK, INFO )
435: *
436: * -- LAPACK computational routine (version 3.6.1) --
437: * -- LAPACK is a software package provided by Univ. of Tennessee, --
438: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
439: * June 2016
440: *
441: * .. Scalar Arguments ..
442: LOGICAL WANTQ, WANTZ
443: INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
444: $ M, N
445: DOUBLE PRECISION PL, PR
446: * ..
447: * .. Array Arguments ..
448: LOGICAL SELECT( * )
449: INTEGER IWORK( * )
450: DOUBLE PRECISION DIF( * )
451: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
452: $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * )
453: * ..
454: *
455: * =====================================================================
456: *
457: * .. Parameters ..
458: INTEGER IDIFJB
459: PARAMETER ( IDIFJB = 3 )
460: DOUBLE PRECISION ZERO, ONE
461: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
462: * ..
463: * .. Local Scalars ..
464: LOGICAL LQUERY, SWAP, WANTD, WANTD1, WANTD2, WANTP
465: INTEGER I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2,
466: $ N1, N2
467: DOUBLE PRECISION DSCALE, DSUM, RDSCAL, SAFMIN
468: COMPLEX*16 TEMP1, TEMP2
469: * ..
470: * .. Local Arrays ..
471: INTEGER ISAVE( 3 )
472: * ..
473: * .. External Subroutines ..
474: EXTERNAL XERBLA, ZLACN2, ZLACPY, ZLASSQ, ZSCAL, ZTGEXC,
475: $ ZTGSYL
476: * ..
477: * .. Intrinsic Functions ..
478: INTRINSIC ABS, DCMPLX, DCONJG, MAX, SQRT
479: * ..
480: * .. External Functions ..
481: DOUBLE PRECISION DLAMCH
482: EXTERNAL DLAMCH
483: * ..
484: * .. Executable Statements ..
485: *
486: * Decode and test the input parameters
487: *
488: INFO = 0
489: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
490: *
491: IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN
492: INFO = -1
493: ELSE IF( N.LT.0 ) THEN
494: INFO = -5
495: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
496: INFO = -7
497: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
498: INFO = -9
499: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
500: INFO = -13
501: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
502: INFO = -15
503: END IF
504: *
505: IF( INFO.NE.0 ) THEN
506: CALL XERBLA( 'ZTGSEN', -INFO )
507: RETURN
508: END IF
509: *
510: IERR = 0
511: *
512: WANTP = IJOB.EQ.1 .OR. IJOB.GE.4
513: WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4
514: WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5
515: WANTD = WANTD1 .OR. WANTD2
516: *
517: * Set M to the dimension of the specified pair of deflating
518: * subspaces.
519: *
520: M = 0
521: IF( .NOT.LQUERY .OR. IJOB.NE.0 ) THEN
522: DO 10 K = 1, N
523: ALPHA( K ) = A( K, K )
524: BETA( K ) = B( K, K )
525: IF( K.LT.N ) THEN
526: IF( SELECT( K ) )
527: $ M = M + 1
528: ELSE
529: IF( SELECT( N ) )
530: $ M = M + 1
531: END IF
532: 10 CONTINUE
533: END IF
534: *
535: IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
536: LWMIN = MAX( 1, 2*M*( N-M ) )
537: LIWMIN = MAX( 1, N+2 )
538: ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
539: LWMIN = MAX( 1, 4*M*( N-M ) )
540: LIWMIN = MAX( 1, 2*M*( N-M ), N+2 )
541: ELSE
542: LWMIN = 1
543: LIWMIN = 1
544: END IF
545: *
546: WORK( 1 ) = LWMIN
547: IWORK( 1 ) = LIWMIN
548: *
549: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
550: INFO = -21
551: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
552: INFO = -23
553: END IF
554: *
555: IF( INFO.NE.0 ) THEN
556: CALL XERBLA( 'ZTGSEN', -INFO )
557: RETURN
558: ELSE IF( LQUERY ) THEN
559: RETURN
560: END IF
561: *
562: * Quick return if possible.
563: *
564: IF( M.EQ.N .OR. M.EQ.0 ) THEN
565: IF( WANTP ) THEN
566: PL = ONE
567: PR = ONE
568: END IF
569: IF( WANTD ) THEN
570: DSCALE = ZERO
571: DSUM = ONE
572: DO 20 I = 1, N
573: CALL ZLASSQ( N, A( 1, I ), 1, DSCALE, DSUM )
574: CALL ZLASSQ( N, B( 1, I ), 1, DSCALE, DSUM )
575: 20 CONTINUE
576: DIF( 1 ) = DSCALE*SQRT( DSUM )
577: DIF( 2 ) = DIF( 1 )
578: END IF
579: GO TO 70
580: END IF
581: *
582: * Get machine constant
583: *
584: SAFMIN = DLAMCH( 'S' )
585: *
586: * Collect the selected blocks at the top-left corner of (A, B).
587: *
588: KS = 0
589: DO 30 K = 1, N
590: SWAP = SELECT( K )
591: IF( SWAP ) THEN
592: KS = KS + 1
593: *
594: * Swap the K-th block to position KS. Compute unitary Q
595: * and Z that will swap adjacent diagonal blocks in (A, B).
596: *
597: IF( K.NE.KS )
598: $ CALL ZTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
599: $ LDZ, K, KS, IERR )
600: *
601: IF( IERR.GT.0 ) THEN
602: *
603: * Swap is rejected: exit.
604: *
605: INFO = 1
606: IF( WANTP ) THEN
607: PL = ZERO
608: PR = ZERO
609: END IF
610: IF( WANTD ) THEN
611: DIF( 1 ) = ZERO
612: DIF( 2 ) = ZERO
613: END IF
614: GO TO 70
615: END IF
616: END IF
617: 30 CONTINUE
618: IF( WANTP ) THEN
619: *
620: * Solve generalized Sylvester equation for R and L:
621: * A11 * R - L * A22 = A12
622: * B11 * R - L * B22 = B12
623: *
624: N1 = M
625: N2 = N - M
626: I = N1 + 1
627: CALL ZLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
628: CALL ZLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
629: $ N1 )
630: IJB = 0
631: CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
632: $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1,
633: $ DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
634: $ LWORK-2*N1*N2, IWORK, IERR )
635: *
636: * Estimate the reciprocal of norms of "projections" onto
637: * left and right eigenspaces
638: *
639: RDSCAL = ZERO
640: DSUM = ONE
641: CALL ZLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM )
642: PL = RDSCAL*SQRT( DSUM )
643: IF( PL.EQ.ZERO ) THEN
644: PL = ONE
645: ELSE
646: PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) )
647: END IF
648: RDSCAL = ZERO
649: DSUM = ONE
650: CALL ZLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM )
651: PR = RDSCAL*SQRT( DSUM )
652: IF( PR.EQ.ZERO ) THEN
653: PR = ONE
654: ELSE
655: PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) )
656: END IF
657: END IF
658: IF( WANTD ) THEN
659: *
660: * Compute estimates Difu and Difl.
661: *
662: IF( WANTD1 ) THEN
663: N1 = M
664: N2 = N - M
665: I = N1 + 1
666: IJB = IDIFJB
667: *
668: * Frobenius norm-based Difu estimate.
669: *
670: CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
671: $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
672: $ N1, DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
673: $ LWORK-2*N1*N2, IWORK, IERR )
674: *
675: * Frobenius norm-based Difl estimate.
676: *
677: CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
678: $ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
679: $ N2, DSCALE, DIF( 2 ), WORK( N1*N2*2+1 ),
680: $ LWORK-2*N1*N2, IWORK, IERR )
681: ELSE
682: *
683: * Compute 1-norm-based estimates of Difu and Difl using
684: * reversed communication with ZLACN2. In each step a
685: * generalized Sylvester equation or a transposed variant
686: * is solved.
687: *
688: KASE = 0
689: N1 = M
690: N2 = N - M
691: I = N1 + 1
692: IJB = 0
693: MN2 = 2*N1*N2
694: *
695: * 1-norm-based estimate of Difu.
696: *
697: 40 CONTINUE
698: CALL ZLACN2( MN2, WORK( MN2+1 ), WORK, DIF( 1 ), KASE,
699: $ ISAVE )
700: IF( KASE.NE.0 ) THEN
701: IF( KASE.EQ.1 ) THEN
702: *
703: * Solve generalized Sylvester equation
704: *
705: CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
706: $ WORK, N1, B, LDB, B( I, I ), LDB,
707: $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
708: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
709: $ IERR )
710: ELSE
711: *
712: * Solve the transposed variant.
713: *
714: CALL ZTGSYL( 'C', IJB, N1, N2, A, LDA, A( I, I ), LDA,
715: $ WORK, N1, B, LDB, B( I, I ), LDB,
716: $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
717: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
718: $ IERR )
719: END IF
720: GO TO 40
721: END IF
722: DIF( 1 ) = DSCALE / DIF( 1 )
723: *
724: * 1-norm-based estimate of Difl.
725: *
726: 50 CONTINUE
727: CALL ZLACN2( MN2, WORK( MN2+1 ), WORK, DIF( 2 ), KASE,
728: $ ISAVE )
729: IF( KASE.NE.0 ) THEN
730: IF( KASE.EQ.1 ) THEN
731: *
732: * Solve generalized Sylvester equation
733: *
734: CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
735: $ WORK, N2, B( I, I ), LDB, B, LDB,
736: $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
737: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
738: $ IERR )
739: ELSE
740: *
741: * Solve the transposed variant.
742: *
743: CALL ZTGSYL( 'C', IJB, N2, N1, A( I, I ), LDA, A, LDA,
744: $ WORK, N2, B, LDB, B( I, I ), LDB,
745: $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
746: $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK,
747: $ IERR )
748: END IF
749: GO TO 50
750: END IF
751: DIF( 2 ) = DSCALE / DIF( 2 )
752: END IF
753: END IF
754: *
755: * If B(K,K) is complex, make it real and positive (normalization
756: * of the generalized Schur form) and Store the generalized
757: * eigenvalues of reordered pair (A, B)
758: *
759: DO 60 K = 1, N
760: DSCALE = ABS( B( K, K ) )
761: IF( DSCALE.GT.SAFMIN ) THEN
762: TEMP1 = DCONJG( B( K, K ) / DSCALE )
763: TEMP2 = B( K, K ) / DSCALE
764: B( K, K ) = DSCALE
765: CALL ZSCAL( N-K, TEMP1, B( K, K+1 ), LDB )
766: CALL ZSCAL( N-K+1, TEMP1, A( K, K ), LDA )
767: IF( WANTQ )
768: $ CALL ZSCAL( N, TEMP2, Q( 1, K ), 1 )
769: ELSE
770: B( K, K ) = DCMPLX( ZERO, ZERO )
771: END IF
772: *
773: ALPHA( K ) = A( K, K )
774: BETA( K ) = B( K, K )
775: *
776: 60 CONTINUE
777: *
778: 70 CONTINUE
779: *
780: WORK( 1 ) = LWMIN
781: IWORK( 1 ) = LIWMIN
782: *
783: RETURN
784: *
785: * End of ZTGSEN
786: *
787: 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