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