LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 June 2016
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.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  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  IF( .NOT.lquery .OR. ijob.NE.0 ) THEN
519  DO 10 k = 1, n
520  alpha( k ) = a( k, k )
521  beta( k ) = b( k, k )
522  IF( k.LT.n ) THEN
523  IF( SELECT( k ) )
524  $ m = m + 1
525  ELSE
526  IF( SELECT( n ) )
527  $ m = m + 1
528  END IF
529  10 CONTINUE
530  END IF
531 *
532  IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
533  lwmin = max( 1, 2*m*(n-m) )
534  liwmin = max( 1, n+2 )
535  ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 ) THEN
536  lwmin = max( 1, 4*m*(n-m) )
537  liwmin = max( 1, 2*m*(n-m), n+2 )
538  ELSE
539  lwmin = 1
540  liwmin = 1
541  END IF
542 *
543  work( 1 ) = lwmin
544  iwork( 1 ) = liwmin
545 *
546  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
547  info = -21
548  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
549  info = -23
550  END IF
551 *
552  IF( info.NE.0 ) THEN
553  CALL xerbla( 'CTGSEN', -info )
554  RETURN
555  ELSE IF( lquery ) THEN
556  RETURN
557  END IF
558 *
559 * Quick return if possible.
560 *
561  IF( m.EQ.n .OR. m.EQ.0 ) THEN
562  IF( wantp ) THEN
563  pl = one
564  pr = one
565  END IF
566  IF( wantd ) THEN
567  dscale = zero
568  dsum = one
569  DO 20 i = 1, n
570  CALL classq( n, a( 1, i ), 1, dscale, dsum )
571  CALL classq( n, b( 1, i ), 1, dscale, dsum )
572  20 CONTINUE
573  dif( 1 ) = dscale*sqrt( dsum )
574  dif( 2 ) = dif( 1 )
575  END IF
576  GO TO 70
577  END IF
578 *
579 * Get machine constant
580 *
581  safmin = slamch( 'S' )
582 *
583 * Collect the selected blocks at the top-left corner of (A, B).
584 *
585  ks = 0
586  DO 30 k = 1, n
587  swap = SELECT( k )
588  IF( swap ) THEN
589  ks = ks + 1
590 *
591 * Swap the K-th block to position KS. Compute unitary Q
592 * and Z that will swap adjacent diagonal blocks in (A, B).
593 *
594  IF( k.NE.ks )
595  $ CALL ctgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,
596  $ ldz, k, ks, ierr )
597 *
598  IF( ierr.GT.0 ) THEN
599 *
600 * Swap is rejected: exit.
601 *
602  info = 1
603  IF( wantp ) THEN
604  pl = zero
605  pr = zero
606  END IF
607  IF( wantd ) THEN
608  dif( 1 ) = zero
609  dif( 2 ) = zero
610  END IF
611  GO TO 70
612  END IF
613  END IF
614  30 CONTINUE
615  IF( wantp ) THEN
616 *
617 * Solve generalized Sylvester equation for R and L:
618 * A11 * R - L * A22 = A12
619 * B11 * R - L * B22 = B12
620 *
621  n1 = m
622  n2 = n - m
623  i = n1 + 1
624  CALL clacpy( 'Full', n1, n2, a( 1, i ), lda, work, n1 )
625  CALL clacpy( 'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
626  $ n1 )
627  ijb = 0
628  CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
629  $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
630  $ dscale, dif( 1 ), work( n1*n2*2+1 ),
631  $ lwork-2*n1*n2, iwork, ierr )
632 *
633 * Estimate the reciprocal of norms of "projections" onto
634 * left and right eigenspaces
635 *
636  rdscal = zero
637  dsum = one
638  CALL classq( n1*n2, work, 1, rdscal, dsum )
639  pl = rdscal*sqrt( dsum )
640  IF( pl.EQ.zero ) THEN
641  pl = one
642  ELSE
643  pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
644  END IF
645  rdscal = zero
646  dsum = one
647  CALL classq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
648  pr = rdscal*sqrt( dsum )
649  IF( pr.EQ.zero ) THEN
650  pr = one
651  ELSE
652  pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
653  END IF
654  END IF
655  IF( wantd ) THEN
656 *
657 * Compute estimates Difu and Difl.
658 *
659  IF( wantd1 ) THEN
660  n1 = m
661  n2 = n - m
662  i = n1 + 1
663  ijb = idifjb
664 *
665 * Frobenius norm-based Difu estimate.
666 *
667  CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
668  $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
669  $ n1, dscale, dif( 1 ), work( n1*n2*2+1 ),
670  $ lwork-2*n1*n2, iwork, ierr )
671 *
672 * Frobenius norm-based Difl estimate.
673 *
674  CALL ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
675  $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
676  $ n2, dscale, dif( 2 ), work( n1*n2*2+1 ),
677  $ lwork-2*n1*n2, iwork, ierr )
678  ELSE
679 *
680 * Compute 1-norm-based estimates of Difu and Difl using
681 * reversed communication with CLACN2. In each step a
682 * generalized Sylvester equation or a transposed variant
683 * is solved.
684 *
685  kase = 0
686  n1 = m
687  n2 = n - m
688  i = n1 + 1
689  ijb = 0
690  mn2 = 2*n1*n2
691 *
692 * 1-norm-based estimate of Difu.
693 *
694  40 CONTINUE
695  CALL clacn2( mn2, work( mn2+1 ), work, dif( 1 ), kase,
696  $ isave )
697  IF( kase.NE.0 ) THEN
698  IF( kase.EQ.1 ) THEN
699 *
700 * Solve generalized Sylvester equation
701 *
702  CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,
703  $ work, n1, b, ldb, b( i, i ), ldb,
704  $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
705  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
706  $ ierr )
707  ELSE
708 *
709 * Solve the transposed variant.
710 *
711  CALL ctgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ), lda,
712  $ work, n1, b, ldb, b( i, i ), ldb,
713  $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
714  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
715  $ ierr )
716  END IF
717  GO TO 40
718  END IF
719  dif( 1 ) = dscale / dif( 1 )
720 *
721 * 1-norm-based estimate of Difl.
722 *
723  50 CONTINUE
724  CALL clacn2( mn2, work( mn2+1 ), work, dif( 2 ), kase,
725  $ isave )
726  IF( kase.NE.0 ) THEN
727  IF( kase.EQ.1 ) THEN
728 *
729 * Solve generalized Sylvester equation
730 *
731  CALL ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,
732  $ work, n2, b( i, i ), ldb, b, ldb,
733  $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
734  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
735  $ ierr )
736  ELSE
737 *
738 * Solve the transposed variant.
739 *
740  CALL ctgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a, lda,
741  $ work, n2, b, ldb, b( i, i ), ldb,
742  $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
743  $ work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
744  $ ierr )
745  END IF
746  GO TO 50
747  END IF
748  dif( 2 ) = dscale / dif( 2 )
749  END IF
750  END IF
751 *
752 * If B(K,K) is complex, make it real and positive (normalization
753 * of the generalized Schur form) and Store the generalized
754 * eigenvalues of reordered pair (A, B)
755 *
756  DO 60 k = 1, n
757  dscale = abs( b( k, k ) )
758  IF( dscale.GT.safmin ) THEN
759  temp1 = conjg( b( k, k ) / dscale )
760  temp2 = b( k, k ) / dscale
761  b( k, k ) = dscale
762  CALL cscal( n-k, temp1, b( k, k+1 ), ldb )
763  CALL cscal( n-k+1, temp1, a( k, k ), lda )
764  IF( wantq )
765  $ CALL cscal( n, temp2, q( 1, k ), 1 )
766  ELSE
767  b( k, k ) = cmplx( zero, zero )
768  END IF
769 *
770  alpha( k ) = a( k, k )
771  beta( k ) = b( k, k )
772 *
773  60 CONTINUE
774 *
775  70 CONTINUE
776 *
777  work( 1 ) = lwmin
778  iwork( 1 ) = liwmin
779 *
780  RETURN
781 *
782 * End of CTGSEN
783 *
784  END
subroutine classq(N, X, INCX, SCALE, SUMSQ)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
subroutine ctgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
CTGSEN
Definition: ctgsen.f:435
subroutine ctgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
CTGSYL
Definition: ctgsyl.f:297
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine ctgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
CTGEXC
Definition: ctgexc.f:202
subroutine clacn2(N, V, X, EST, KASE, ISAVE)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: clacn2.f:135