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