LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
ctgsen.f
Go to the documentation of this file.
1 *> \brief \b CTGSEN
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CTGSEN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsen.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsen.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsen.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CTGSEN( 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 * REAL PL, PR
30 * ..
31 * .. Array Arguments ..
32 * LOGICAL SELECT( * )
33 * INTEGER IWORK( * )
34 * REAL DIF( * )
35 * COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
36 * $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> CTGSEN 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 *> CTGSEN 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 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 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 array, dimension (N)
151 *> \endverbatim
152 *>
153 *> \param[out] BETA
154 *> \verbatim
155 *> BETA is COMPLEX 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 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 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 REAL
209 *> \endverbatim
210 *>
211 *> \param[out] PR
212 *> \verbatim
213 *> PR is REAL
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 REAL 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 CLACN2.
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 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 complexOTHERcomputational
296 *
297 *> \par Further Details:
298 * =====================
299 *>
300 *> \verbatim
301 *>
302 *> CTGSEN 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', 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 *> conjuguate 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 CLATDF), then the parameter
398 *> IDIFJB (see below) should be changed from 3 to 4 (routine CLATDF
399 *> (IJOB = 2 will be used)). See CTGSYL 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 ctgsen( 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  REAL pl, pr
446 * ..
447 * .. Array Arguments ..
448  LOGICAL select( * )
449  INTEGER iwork( * )
450  REAL dif( * )
451  COMPLEX 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  REAL zero, one
461  parameter( zero = 0.0e+0, one = 1.0e+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  REAL dscale, dsum, rdscal, safmin
468  COMPLEX temp1, temp2
469 * ..
470 * .. Local Arrays ..
471  INTEGER isave( 3 )
472 * ..
473 * .. External Subroutines ..
474  REAL slamch
475  EXTERNAL clacn2, clacpy, classq, cscal, ctgexc, ctgsyl,
476  $ slamch, xerbla
477 * ..
478 * .. Intrinsic Functions ..
479  INTRINSIC abs, cmplx, conjg, max, sqrt
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( 'CTGSEN', -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
518  DO 10 k = 1, n
519  alpha( k ) = a( k, k )
520  beta( k ) = b( k, k )
521  IF( k.LT.n ) THEN
522  IF( SELECT( k ) )
523  $ m = m + 1
524  ELSE
525  IF( SELECT( n ) )
526  $ m = m + 1
527  END IF
528  10 continue
529 *
530  IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
531  lwmin = max( 1, 2*m*(n-m) )
532  liwmin = max( 1, n+2 )
533  ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 ) THEN
534  lwmin = max( 1, 4*m*(n-m) )
535  liwmin = max( 1, 2*m*(n-m), n+2 )
536  ELSE
537  lwmin = 1
538  liwmin = 1
539  END IF
540 *
541  work( 1 ) = lwmin
542  iwork( 1 ) = liwmin
543 *
544  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
545  info = -21
546  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
547  info = -23
548  END IF
549 *
550  IF( info.NE.0 ) THEN
551  CALL xerbla( 'CTGSEN', -info )
552  return
553  ELSE IF( lquery ) THEN
554  return
555  END IF
556 *
557 * Quick return if possible.
558 *
559  IF( m.EQ.n .OR. m.EQ.0 ) THEN
560  IF( wantp ) THEN
561  pl = one
562  pr = one
563  END IF
564  IF( wantd ) THEN
565  dscale = zero
566  dsum = one
567  DO 20 i = 1, n
568  CALL classq( n, a( 1, i ), 1, dscale, dsum )
569  CALL classq( n, b( 1, i ), 1, dscale, dsum )
570  20 continue
571  dif( 1 ) = dscale*sqrt( dsum )
572  dif( 2 ) = dif( 1 )
573  END IF
574  go to 70
575  END IF
576 *
577 * Get machine constant
578 *
579  safmin = slamch( 'S' )
580 *
581 * Collect the selected blocks at the top-left corner of (A, B).
582 *
583  ks = 0
584  DO 30 k = 1, n
585  swap = SELECT( k )
586  IF( swap ) THEN
587  ks = ks + 1
588 *
589 * Swap the K-th block to position KS. Compute unitary Q
590 * and Z that will swap adjacent diagonal blocks in (A, B).
591 *
592  IF( k.NE.ks )
593  $ CALL ctgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
594  $ ldz, k, ks, ierr )
595 *
596  IF( ierr.GT.0 ) THEN
597 *
598 * Swap is rejected: exit.
599 *
600  info = 1
601  IF( wantp ) THEN
602  pl = zero
603  pr = zero
604  END IF
605  IF( wantd ) THEN
606  dif( 1 ) = zero
607  dif( 2 ) = zero
608  END IF
609  go to 70
610  END IF
611  END IF
612  30 continue
613  IF( wantp ) THEN
614 *
615 * Solve generalized Sylvester equation for R and L:
616 * A11 * R - L * A22 = A12
617 * B11 * R - L * B22 = B12
618 *
619  n1 = m
620  n2 = n - m
621  i = n1 + 1
622  CALL clacpy( 'Full', n1, n2, a( 1, i ), lda, work, n1 )
623  CALL clacpy( 'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
624  $ n1 )
625  ijb = 0
626  CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
627  $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
628  $ dscale, dif( 1 ), work( n1*n2*2+1 ),
629  $ lwork-2*n1*n2, iwork, ierr )
630 *
631 * Estimate the reciprocal of norms of "projections" onto
632 * left and right eigenspaces
633 *
634  rdscal = zero
635  dsum = one
636  CALL classq( n1*n2, work, 1, rdscal, dsum )
637  pl = rdscal*sqrt( dsum )
638  IF( pl.EQ.zero ) THEN
639  pl = one
640  ELSE
641  pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
642  END IF
643  rdscal = zero
644  dsum = one
645  CALL classq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
646  pr = rdscal*sqrt( dsum )
647  IF( pr.EQ.zero ) THEN
648  pr = one
649  ELSE
650  pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
651  END IF
652  END IF
653  IF( wantd ) THEN
654 *
655 * Compute estimates Difu and Difl.
656 *
657  IF( wantd1 ) THEN
658  n1 = m
659  n2 = n - m
660  i = n1 + 1
661  ijb = idifjb
662 *
663 * Frobenius norm-based Difu estimate.
664 *
665  CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
666  $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
667  $ n1, dscale, dif( 1 ), work( n1*n2*2+1 ),
668  $ lwork-2*n1*n2, iwork, ierr )
669 *
670 * Frobenius norm-based Difl estimate.
671 *
672  CALL ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
673  $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
674  $ n2, dscale, dif( 2 ), work( n1*n2*2+1 ),
675  $ lwork-2*n1*n2, iwork, ierr )
676  ELSE
677 *
678 * Compute 1-norm-based estimates of Difu and Difl using
679 * reversed communication with CLACN2. In each step a
680 * generalized Sylvester equation or a transposed variant
681 * is solved.
682 *
683  kase = 0
684  n1 = m
685  n2 = n - m
686  i = n1 + 1
687  ijb = 0
688  mn2 = 2*n1*n2
689 *
690 * 1-norm-based estimate of Difu.
691 *
692  40 continue
693  CALL clacn2( mn2, work( mn2+1 ), work, dif( 1 ), kase,
694  $ isave )
695  IF( kase.NE.0 ) THEN
696  IF( kase.EQ.1 ) THEN
697 *
698 * Solve generalized Sylvester equation
699 *
700  CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,
701  $ work, n1, b, ldb, b( i, i ), ldb,
702  $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
703  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
704  $ ierr )
705  ELSE
706 *
707 * Solve the transposed variant.
708 *
709  CALL ctgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ), lda,
710  $ work, n1, b, ldb, b( i, i ), ldb,
711  $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
712  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
713  $ ierr )
714  END IF
715  go to 40
716  END IF
717  dif( 1 ) = dscale / dif( 1 )
718 *
719 * 1-norm-based estimate of Difl.
720 *
721  50 continue
722  CALL clacn2( mn2, work( mn2+1 ), work, dif( 2 ), kase,
723  $ isave )
724  IF( kase.NE.0 ) THEN
725  IF( kase.EQ.1 ) THEN
726 *
727 * Solve generalized Sylvester equation
728 *
729  CALL ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,
730  $ work, n2, b( i, i ), ldb, b, ldb,
731  $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
732  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
733  $ ierr )
734  ELSE
735 *
736 * Solve the transposed variant.
737 *
738  CALL ctgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a, lda,
739  $ work, n2, b, ldb, b( i, i ), ldb,
740  $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
741  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
742  $ ierr )
743  END IF
744  go to 50
745  END IF
746  dif( 2 ) = dscale / dif( 2 )
747  END IF
748  END IF
749 *
750 * If B(K,K) is complex, make it real and positive (normalization
751 * of the generalized Schur form) and Store the generalized
752 * eigenvalues of reordered pair (A, B)
753 *
754  DO 60 k = 1, n
755  dscale = abs( b( k, k ) )
756  IF( dscale.GT.safmin ) THEN
757  temp1 = conjg( b( k, k ) / dscale )
758  temp2 = b( k, k ) / dscale
759  b( k, k ) = dscale
760  CALL cscal( n-k, temp1, b( k, k+1 ), ldb )
761  CALL cscal( n-k+1, temp1, a( k, k ), lda )
762  IF( wantq )
763  $ CALL cscal( n, temp2, q( 1, k ), 1 )
764  ELSE
765  b( k, k ) = cmplx( zero, zero )
766  END IF
767 *
768  alpha( k ) = a( k, k )
769  beta( k ) = b( k, k )
770 *
771  60 continue
772 *
773  70 continue
774 *
775  work( 1 ) = lwmin
776  iwork( 1 ) = liwmin
777 *
778  return
779 *
780 * End of CTGSEN
781 *
782  END