LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
ztgsen.f
Go to the documentation of this file.
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 November 2011
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.4.0) --
437 * -- LAPACK is a software package provided by Univ. of Tennessee, --
438 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
439 * November 2011
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  DO 10 k = 1, n
522  alpha( k ) = a( k, k )
523  beta( k ) = b( k, k )
524  IF( k.LT.n ) THEN
525  IF( SELECT( k ) )
526  $ m = m + 1
527  ELSE
528  IF( SELECT( n ) )
529  $ m = m + 1
530  END IF
531  10 continue
532 *
533  IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
534  lwmin = max( 1, 2*m*( n-m ) )
535  liwmin = max( 1, n+2 )
536  ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 ) THEN
537  lwmin = max( 1, 4*m*( n-m ) )
538  liwmin = max( 1, 2*m*( n-m ), n+2 )
539  ELSE
540  lwmin = 1
541  liwmin = 1
542  END IF
543 *
544  work( 1 ) = lwmin
545  iwork( 1 ) = liwmin
546 *
547  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
548  info = -21
549  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
550  info = -23
551  END IF
552 *
553  IF( info.NE.0 ) THEN
554  CALL xerbla( 'ZTGSEN', -info )
555  return
556  ELSE IF( lquery ) THEN
557  return
558  END IF
559 *
560 * Quick return if possible.
561 *
562  IF( m.EQ.n .OR. m.EQ.0 ) THEN
563  IF( wantp ) THEN
564  pl = one
565  pr = one
566  END IF
567  IF( wantd ) THEN
568  dscale = zero
569  dsum = one
570  DO 20 i = 1, n
571  CALL zlassq( n, a( 1, i ), 1, dscale, dsum )
572  CALL zlassq( n, b( 1, i ), 1, dscale, dsum )
573  20 continue
574  dif( 1 ) = dscale*sqrt( dsum )
575  dif( 2 ) = dif( 1 )
576  END IF
577  go to 70
578  END IF
579 *
580 * Get machine constant
581 *
582  safmin = dlamch( 'S' )
583 *
584 * Collect the selected blocks at the top-left corner of (A, B).
585 *
586  ks = 0
587  DO 30 k = 1, n
588  swap = SELECT( k )
589  IF( swap ) THEN
590  ks = ks + 1
591 *
592 * Swap the K-th block to position KS. Compute unitary Q
593 * and Z that will swap adjacent diagonal blocks in (A, B).
594 *
595  IF( k.NE.ks )
596  $ CALL ztgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
597  $ ldz, k, ks, ierr )
598 *
599  IF( ierr.GT.0 ) THEN
600 *
601 * Swap is rejected: exit.
602 *
603  info = 1
604  IF( wantp ) THEN
605  pl = zero
606  pr = zero
607  END IF
608  IF( wantd ) THEN
609  dif( 1 ) = zero
610  dif( 2 ) = zero
611  END IF
612  go to 70
613  END IF
614  END IF
615  30 continue
616  IF( wantp ) THEN
617 *
618 * Solve generalized Sylvester equation for R and L:
619 * A11 * R - L * A22 = A12
620 * B11 * R - L * B22 = B12
621 *
622  n1 = m
623  n2 = n - m
624  i = n1 + 1
625  CALL zlacpy( 'Full', n1, n2, a( 1, i ), lda, work, n1 )
626  CALL zlacpy( 'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
627  $ n1 )
628  ijb = 0
629  CALL ztgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
630  $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
631  $ dscale, dif( 1 ), work( n1*n2*2+1 ),
632  $ lwork-2*n1*n2, iwork, ierr )
633 *
634 * Estimate the reciprocal of norms of "projections" onto
635 * left and right eigenspaces
636 *
637  rdscal = zero
638  dsum = one
639  CALL zlassq( n1*n2, work, 1, rdscal, dsum )
640  pl = rdscal*sqrt( dsum )
641  IF( pl.EQ.zero ) THEN
642  pl = one
643  ELSE
644  pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
645  END IF
646  rdscal = zero
647  dsum = one
648  CALL zlassq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
649  pr = rdscal*sqrt( dsum )
650  IF( pr.EQ.zero ) THEN
651  pr = one
652  ELSE
653  pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
654  END IF
655  END IF
656  IF( wantd ) THEN
657 *
658 * Compute estimates Difu and Difl.
659 *
660  IF( wantd1 ) THEN
661  n1 = m
662  n2 = n - m
663  i = n1 + 1
664  ijb = idifjb
665 *
666 * Frobenius norm-based Difu estimate.
667 *
668  CALL ztgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
669  $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
670  $ n1, dscale, dif( 1 ), work( n1*n2*2+1 ),
671  $ lwork-2*n1*n2, iwork, ierr )
672 *
673 * Frobenius norm-based Difl estimate.
674 *
675  CALL ztgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
676  $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
677  $ n2, dscale, dif( 2 ), work( n1*n2*2+1 ),
678  $ lwork-2*n1*n2, iwork, ierr )
679  ELSE
680 *
681 * Compute 1-norm-based estimates of Difu and Difl using
682 * reversed communication with ZLACN2. In each step a
683 * generalized Sylvester equation or a transposed variant
684 * is solved.
685 *
686  kase = 0
687  n1 = m
688  n2 = n - m
689  i = n1 + 1
690  ijb = 0
691  mn2 = 2*n1*n2
692 *
693 * 1-norm-based estimate of Difu.
694 *
695  40 continue
696  CALL zlacn2( mn2, work( mn2+1 ), work, dif( 1 ), kase,
697  $ isave )
698  IF( kase.NE.0 ) THEN
699  IF( kase.EQ.1 ) THEN
700 *
701 * Solve generalized Sylvester equation
702 *
703  CALL ztgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,
704  $ work, n1, b, ldb, b( i, i ), ldb,
705  $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
706  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
707  $ ierr )
708  ELSE
709 *
710 * Solve the transposed variant.
711 *
712  CALL ztgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ), lda,
713  $ work, n1, b, ldb, b( i, i ), ldb,
714  $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
715  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
716  $ ierr )
717  END IF
718  go to 40
719  END IF
720  dif( 1 ) = dscale / dif( 1 )
721 *
722 * 1-norm-based estimate of Difl.
723 *
724  50 continue
725  CALL zlacn2( mn2, work( mn2+1 ), work, dif( 2 ), kase,
726  $ isave )
727  IF( kase.NE.0 ) THEN
728  IF( kase.EQ.1 ) THEN
729 *
730 * Solve generalized Sylvester equation
731 *
732  CALL ztgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,
733  $ work, n2, b( i, i ), ldb, b, ldb,
734  $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
735  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
736  $ ierr )
737  ELSE
738 *
739 * Solve the transposed variant.
740 *
741  CALL ztgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a, lda,
742  $ work, n2, b, ldb, b( i, i ), ldb,
743  $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
744  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
745  $ ierr )
746  END IF
747  go to 50
748  END IF
749  dif( 2 ) = dscale / dif( 2 )
750  END IF
751  END IF
752 *
753 * If B(K,K) is complex, make it real and positive (normalization
754 * of the generalized Schur form) and Store the generalized
755 * eigenvalues of reordered pair (A, B)
756 *
757  DO 60 k = 1, n
758  dscale = abs( b( k, k ) )
759  IF( dscale.GT.safmin ) THEN
760  temp1 = dconjg( b( k, k ) / dscale )
761  temp2 = b( k, k ) / dscale
762  b( k, k ) = dscale
763  CALL zscal( n-k, temp1, b( k, k+1 ), ldb )
764  CALL zscal( n-k+1, temp1, a( k, k ), lda )
765  IF( wantq )
766  $ CALL zscal( n, temp2, q( 1, k ), 1 )
767  ELSE
768  b( k, k ) = dcmplx( zero, zero )
769  END IF
770 *
771  alpha( k ) = a( k, k )
772  beta( k ) = b( k, k )
773 *
774  60 continue
775 *
776  70 continue
777 *
778  work( 1 ) = lwmin
779  iwork( 1 ) = liwmin
780 *
781  return
782 *
783 * End of ZTGSEN
784 *
785  END