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