LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
stgsen.f
Go to the documentation of this file.
1 *> \brief \b STGSEN
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download STGSEN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsen.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsen.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsen.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STGSEN( 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 * REAL PL, PR
30 * ..
31 * .. Array Arguments ..
32 * LOGICAL SELECT( * )
33 * INTEGER IWORK( * )
34 * REAL 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 *> STGSEN 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 SGGES), 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 *> STGSEN 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, STGSEN 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 REAL 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 REAL 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 REAL array, dimension (N)
156 *> \endverbatim
157 *>
158 *> \param[out] ALPHAI
159 *> \verbatim
160 *> ALPHAI is REAL array, dimension (N)
161 *> \endverbatim
162 *>
163 *> \param[out] BETA
164 *> \verbatim
165 *> BETA is REAL 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 REAL 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 REAL 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 REAL
224 *> \endverbatim
225 *>
226 *> \param[out] PR
227 *> \verbatim
228 *> PR is REAL
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 REAL 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 REAL array, dimension (MAX(1,LWORK))
252 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
253 *> \endverbatim
254 *>
255 *> \param[in] LWORK
256 *> \verbatim
257 *> LWORK is INTEGER
258 *> The dimension of the array WORK. LWORK >= 4*N+16.
259 *> If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)).
260 *> If IJOB = 3 or 5, LWORK >= MAX(4*N+16, 4*M*(N-M)).
261 *>
262 *> If LWORK = -1, then a workspace query is assumed; the routine
263 *> only calculates the optimal size of the WORK array, returns
264 *> this value as the first entry of the WORK array, and no error
265 *> message related to LWORK is issued by XERBLA.
266 *> \endverbatim
267 *>
268 *> \param[out] IWORK
269 *> \verbatim
270 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
271 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
272 *> \endverbatim
273 *>
274 *> \param[in] LIWORK
275 *> \verbatim
276 *> LIWORK is INTEGER
277 *> The dimension of the array IWORK. LIWORK >= 1.
278 *> If IJOB = 1, 2 or 4, LIWORK >= N+6.
279 *> If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-M), N+6).
280 *>
281 *> If LIWORK = -1, then a workspace query is assumed; the
282 *> routine only calculates the optimal size of the IWORK array,
283 *> returns this value as the first entry of the IWORK array, and
284 *> no error message related to LIWORK is issued by XERBLA.
285 *> \endverbatim
286 *>
287 *> \param[out] INFO
288 *> \verbatim
289 *> INFO is INTEGER
290 *> =0: Successful exit.
291 *> <0: If INFO = -i, the i-th argument had an illegal value.
292 *> =1: Reordering of (A, B) failed because the transformed
293 *> matrix pair (A, B) would be too far from generalized
294 *> Schur form; the problem is very ill-conditioned.
295 *> (A, B) may have been partially reordered.
296 *> If requested, 0 is returned in DIF(*), PL and PR.
297 *> \endverbatim
298 *
299 * Authors:
300 * ========
301 *
302 *> \author Univ. of Tennessee
303 *> \author Univ. of California Berkeley
304 *> \author Univ. of Colorado Denver
305 *> \author NAG Ltd.
306 *
307 *> \date June 2016
308 *
309 *> \ingroup realOTHERcomputational
310 *
311 *> \par Further Details:
312 * =====================
313 *>
314 *> \verbatim
315 *>
316 *> STGSEN first collects the selected eigenvalues by computing
317 *> orthogonal U and W that move them to the top left corner of (A, B).
318 *> In other words, the selected eigenvalues are the eigenvalues of
319 *> (A11, B11) in:
320 *>
321 *> U**T*(A, B)*W = (A11 A12) (B11 B12) n1
322 *> ( 0 A22),( 0 B22) n2
323 *> n1 n2 n1 n2
324 *>
325 *> where N = n1+n2 and U**T means the transpose of U. The first n1 columns
326 *> of U and W span the specified pair of left and right eigenspaces
327 *> (deflating subspaces) of (A, B).
328 *>
329 *> If (A, B) has been obtained from the generalized real Schur
330 *> decomposition of a matrix pair (C, D) = Q*(A, B)*Z**T, then the
331 *> reordered generalized real Schur form of (C, D) is given by
332 *>
333 *> (C, D) = (Q*U)*(U**T*(A, B)*W)*(Z*W)**T,
334 *>
335 *> and the first n1 columns of Q*U and Z*W span the corresponding
336 *> deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
337 *>
338 *> Note that if the selected eigenvalue is sufficiently ill-conditioned,
339 *> then its value may differ significantly from its value before
340 *> reordering.
341 *>
342 *> The reciprocal condition numbers of the left and right eigenspaces
343 *> spanned by the first n1 columns of U and W (or Q*U and Z*W) may
344 *> be returned in DIF(1:2), corresponding to Difu and Difl, resp.
345 *>
346 *> The Difu and Difl are defined as:
347 *>
348 *> Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
349 *> and
350 *> Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
351 *>
352 *> where sigma-min(Zu) is the smallest singular value of the
353 *> (2*n1*n2)-by-(2*n1*n2) matrix
354 *>
355 *> Zu = [ kron(In2, A11) -kron(A22**T, In1) ]
356 *> [ kron(In2, B11) -kron(B22**T, In1) ].
357 *>
358 *> Here, Inx is the identity matrix of size nx and A22**T is the
359 *> transpose of A22. kron(X, Y) is the Kronecker product between
360 *> the matrices X and Y.
361 *>
362 *> When DIF(2) is small, small changes in (A, B) can cause large changes
363 *> in the deflating subspace. An approximate (asymptotic) bound on the
364 *> maximum angular error in the computed deflating subspaces is
365 *>
366 *> EPS * norm((A, B)) / DIF(2),
367 *>
368 *> where EPS is the machine precision.
369 *>
370 *> The reciprocal norm of the projectors on the left and right
371 *> eigenspaces associated with (A11, B11) may be returned in PL and PR.
372 *> They are computed as follows. First we compute L and R so that
373 *> P*(A, B)*Q is block diagonal, where
374 *>
375 *> P = ( I -L ) n1 Q = ( I R ) n1
376 *> ( 0 I ) n2 and ( 0 I ) n2
377 *> n1 n2 n1 n2
378 *>
379 *> and (L, R) is the solution to the generalized Sylvester equation
380 *>
381 *> A11*R - L*A22 = -A12
382 *> B11*R - L*B22 = -B12
383 *>
384 *> Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
385 *> An approximate (asymptotic) bound on the average absolute error of
386 *> the selected eigenvalues is
387 *>
388 *> EPS * norm((A, B)) / PL.
389 *>
390 *> There are also global error bounds which valid for perturbations up
391 *> to a certain restriction: A lower bound (x) on the smallest
392 *> F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
393 *> coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
394 *> (i.e. (A + E, B + F), is
395 *>
396 *> x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
397 *>
398 *> An approximate bound on x can be computed from DIF(1:2), PL and PR.
399 *>
400 *> If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
401 *> (L', R') and unperturbed (L, R) left and right deflating subspaces
402 *> associated with the selected cluster in the (1,1)-blocks can be
403 *> bounded as
404 *>
405 *> max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
406 *> max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
407 *>
408 *> See LAPACK User's Guide section 4.11 or the following references
409 *> for more information.
410 *>
411 *> Note that if the default method for computing the Frobenius-norm-
412 *> based estimate DIF is not wanted (see SLATDF), then the parameter
413 *> IDIFJB (see below) should be changed from 3 to 4 (routine SLATDF
414 *> (IJOB = 2 will be used)). See STGSYL for more details.
415 *> \endverbatim
416 *
417 *> \par Contributors:
418 * ==================
419 *>
420 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
421 *> Umea University, S-901 87 Umea, Sweden.
422 *
423 *> \par References:
424 * ================
425 *>
426 *> \verbatim
427 *>
428 *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
429 *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
430 *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and
431 *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
432 *>
433 *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
434 *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition
435 *> Estimation: Theory, Algorithms and Software,
436 *> Report UMINF - 94.04, Department of Computing Science, Umea
437 *> University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
438 *> Note 87. To appear in Numerical Algorithms, 1996.
439 *>
440 *> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
441 *> for Solving the Generalized Sylvester Equation and Estimating the
442 *> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
443 *> Department of Computing Science, Umea University, S-901 87 Umea,
444 *> Sweden, December 1993, Revised April 1994, Also as LAPACK Working
445 *> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
446 *> 1996.
447 *> \endverbatim
448 *>
449 * =====================================================================
450  SUBROUTINE stgsen( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
451  $ alphar, alphai, beta, q, ldq, z, ldz, m, pl,
452  $ pr, dif, work, lwork, iwork, liwork, info )
453 *
454 * -- LAPACK computational routine (version 3.6.1) --
455 * -- LAPACK is a software package provided by Univ. of Tennessee, --
456 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
457 * June 2016
458 *
459 * .. Scalar Arguments ..
460  LOGICAL WANTQ, WANTZ
461  INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
462  $ m, n
463  REAL PL, PR
464 * ..
465 * .. Array Arguments ..
466  LOGICAL SELECT( * )
467  INTEGER IWORK( * )
468  REAL A( lda, * ), ALPHAI( * ), ALPHAR( * ),
469  $ b( ldb, * ), beta( * ), dif( * ), q( ldq, * ),
470  $ work( * ), z( ldz, * )
471 * ..
472 *
473 * =====================================================================
474 *
475 * .. Parameters ..
476  INTEGER IDIFJB
477  parameter ( idifjb = 3 )
478  REAL ZERO, ONE
479  parameter ( zero = 0.0e+0, one = 1.0e+0 )
480 * ..
481 * .. Local Scalars ..
482  LOGICAL LQUERY, PAIR, SWAP, WANTD, WANTD1, WANTD2,
483  $ wantp
484  INTEGER I, IERR, IJB, K, KASE, KK, KS, LIWMIN, LWMIN,
485  $ mn2, n1, n2
486  REAL DSCALE, DSUM, EPS, RDSCAL, SMLNUM
487 * ..
488 * .. Local Arrays ..
489  INTEGER ISAVE( 3 )
490 * ..
491 * .. External Subroutines ..
492  EXTERNAL slacn2, slacpy, slag2, slassq, stgexc, stgsyl,
493  $ xerbla
494 * ..
495 * .. External Functions ..
496  REAL SLAMCH
497  EXTERNAL slamch
498 * ..
499 * .. Intrinsic Functions ..
500  INTRINSIC max, sign, sqrt
501 * ..
502 * .. Executable Statements ..
503 *
504 * Decode and test the input parameters
505 *
506  info = 0
507  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
508 *
509  IF( ijob.LT.0 .OR. ijob.GT.5 ) THEN
510  info = -1
511  ELSE IF( n.LT.0 ) THEN
512  info = -5
513  ELSE IF( lda.LT.max( 1, n ) ) THEN
514  info = -7
515  ELSE IF( ldb.LT.max( 1, n ) ) THEN
516  info = -9
517  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
518  info = -14
519  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
520  info = -16
521  END IF
522 *
523  IF( info.NE.0 ) THEN
524  CALL xerbla( 'STGSEN', -info )
525  RETURN
526  END IF
527 *
528 * Get machine constants
529 *
530  eps = slamch( 'P' )
531  smlnum = slamch( 'S' ) / eps
532  ierr = 0
533 *
534  wantp = ijob.EQ.1 .OR. ijob.GE.4
535  wantd1 = ijob.EQ.2 .OR. ijob.EQ.4
536  wantd2 = ijob.EQ.3 .OR. ijob.EQ.5
537  wantd = wantd1 .OR. wantd2
538 *
539 * Set M to the dimension of the specified pair of deflating
540 * subspaces.
541 *
542  m = 0
543  pair = .false.
544  IF( .NOT.lquery .OR. ijob.NE.0 ) THEN
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  END IF
565 *
566  IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
567  lwmin = max( 1, 4*n+16, 2*m*(n-m) )
568  liwmin = max( 1, n+6 )
569  ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 ) THEN
570  lwmin = max( 1, 4*n+16, 4*m*(n-m) )
571  liwmin = max( 1, 2*m*(n-m), n+6 )
572  ELSE
573  lwmin = max( 1, 4*n+16 )
574  liwmin = 1
575  END IF
576 *
577  work( 1 ) = lwmin
578  iwork( 1 ) = liwmin
579 *
580  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
581  info = -22
582  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
583  info = -24
584  END IF
585 *
586  IF( info.NE.0 ) THEN
587  CALL xerbla( 'STGSEN', -info )
588  RETURN
589  ELSE IF( lquery ) THEN
590  RETURN
591  END IF
592 *
593 * Quick return if possible.
594 *
595  IF( m.EQ.n .OR. m.EQ.0 ) THEN
596  IF( wantp ) THEN
597  pl = one
598  pr = one
599  END IF
600  IF( wantd ) THEN
601  dscale = zero
602  dsum = one
603  DO 20 i = 1, n
604  CALL slassq( n, a( 1, i ), 1, dscale, dsum )
605  CALL slassq( n, b( 1, i ), 1, dscale, dsum )
606  20 CONTINUE
607  dif( 1 ) = dscale*sqrt( dsum )
608  dif( 2 ) = dif( 1 )
609  END IF
610  GO TO 60
611  END IF
612 *
613 * Collect the selected blocks at the top-left corner of (A, B).
614 *
615  ks = 0
616  pair = .false.
617  DO 30 k = 1, n
618  IF( pair ) THEN
619  pair = .false.
620  ELSE
621 *
622  swap = SELECT( k )
623  IF( k.LT.n ) THEN
624  IF( a( k+1, k ).NE.zero ) THEN
625  pair = .true.
626  swap = swap .OR. SELECT( k+1 )
627  END IF
628  END IF
629 *
630  IF( swap ) THEN
631  ks = ks + 1
632 *
633 * Swap the K-th block to position KS.
634 * Perform the reordering of diagonal blocks in (A, B)
635 * by orthogonal transformation matrices and update
636 * Q and Z accordingly (if requested):
637 *
638  kk = k
639  IF( k.NE.ks )
640  $ CALL stgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq,
641  $ z, ldz, kk, ks, work, lwork, ierr )
642 *
643  IF( ierr.GT.0 ) THEN
644 *
645 * Swap is rejected: exit.
646 *
647  info = 1
648  IF( wantp ) THEN
649  pl = zero
650  pr = zero
651  END IF
652  IF( wantd ) THEN
653  dif( 1 ) = zero
654  dif( 2 ) = zero
655  END IF
656  GO TO 60
657  END IF
658 *
659  IF( pair )
660  $ ks = ks + 1
661  END IF
662  END IF
663  30 CONTINUE
664  IF( wantp ) THEN
665 *
666 * Solve generalized Sylvester equation for R and L
667 * and compute PL and PR.
668 *
669  n1 = m
670  n2 = n - m
671  i = n1 + 1
672  ijb = 0
673  CALL slacpy( 'Full', n1, n2, a( 1, i ), lda, work, n1 )
674  CALL slacpy( 'Full', n1, n2, b( 1, i ), ldb, work( n1*n2+1 ),
675  $ n1 )
676  CALL stgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
677  $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
678  $ dscale, dif( 1 ), work( n1*n2*2+1 ),
679  $ lwork-2*n1*n2, iwork, ierr )
680 *
681 * Estimate the reciprocal of norms of "projections" onto left
682 * and right eigenspaces.
683 *
684  rdscal = zero
685  dsum = one
686  CALL slassq( n1*n2, work, 1, rdscal, dsum )
687  pl = rdscal*sqrt( dsum )
688  IF( pl.EQ.zero ) THEN
689  pl = one
690  ELSE
691  pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
692  END IF
693  rdscal = zero
694  dsum = one
695  CALL slassq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
696  pr = rdscal*sqrt( dsum )
697  IF( pr.EQ.zero ) THEN
698  pr = one
699  ELSE
700  pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
701  END IF
702  END IF
703 *
704  IF( wantd ) THEN
705 *
706 * Compute estimates of Difu and Difl.
707 *
708  IF( wantd1 ) THEN
709  n1 = m
710  n2 = n - m
711  i = n1 + 1
712  ijb = idifjb
713 *
714 * Frobenius norm-based Difu-estimate.
715 *
716  CALL stgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
717  $ n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
718  $ n1, dscale, dif( 1 ), work( 2*n1*n2+1 ),
719  $ lwork-2*n1*n2, iwork, ierr )
720 *
721 * Frobenius norm-based Difl-estimate.
722 *
723  CALL stgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,
724  $ n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
725  $ n2, dscale, dif( 2 ), work( 2*n1*n2+1 ),
726  $ lwork-2*n1*n2, iwork, ierr )
727  ELSE
728 *
729 *
730 * Compute 1-norm-based estimates of Difu and Difl using
731 * reversed communication with SLACN2. In each step a
732 * generalized Sylvester equation or a transposed variant
733 * is solved.
734 *
735  kase = 0
736  n1 = m
737  n2 = n - m
738  i = n1 + 1
739  ijb = 0
740  mn2 = 2*n1*n2
741 *
742 * 1-norm-based estimate of Difu.
743 *
744  40 CONTINUE
745  CALL slacn2( mn2, work( mn2+1 ), work, iwork, dif( 1 ),
746  $ kase, isave )
747  IF( kase.NE.0 ) THEN
748  IF( kase.EQ.1 ) THEN
749 *
750 * Solve generalized Sylvester equation.
751 *
752  CALL stgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,
753  $ work, n1, b, ldb, b( i, i ), ldb,
754  $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
755  $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
756  $ ierr )
757  ELSE
758 *
759 * Solve the transposed variant.
760 *
761  CALL stgsyl( 'T', ijb, n1, n2, a, lda, a( i, i ), lda,
762  $ work, n1, b, ldb, b( i, i ), ldb,
763  $ work( n1*n2+1 ), n1, dscale, dif( 1 ),
764  $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
765  $ ierr )
766  END IF
767  GO TO 40
768  END IF
769  dif( 1 ) = dscale / dif( 1 )
770 *
771 * 1-norm-based estimate of Difl.
772 *
773  50 CONTINUE
774  CALL slacn2( mn2, work( mn2+1 ), work, iwork, dif( 2 ),
775  $ kase, isave )
776  IF( kase.NE.0 ) THEN
777  IF( kase.EQ.1 ) THEN
778 *
779 * Solve generalized Sylvester equation.
780 *
781  CALL stgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,
782  $ work, n2, b( i, i ), ldb, b, ldb,
783  $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
784  $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
785  $ ierr )
786  ELSE
787 *
788 * Solve the transposed variant.
789 *
790  CALL stgsyl( 'T', ijb, n2, n1, a( i, i ), lda, a, lda,
791  $ work, n2, b( i, i ), ldb, b, ldb,
792  $ work( n1*n2+1 ), n2, dscale, dif( 2 ),
793  $ work( 2*n1*n2+1 ), lwork-2*n1*n2, iwork,
794  $ ierr )
795  END IF
796  GO TO 50
797  END IF
798  dif( 2 ) = dscale / dif( 2 )
799 *
800  END IF
801  END IF
802 *
803  60 CONTINUE
804 *
805 * Compute generalized eigenvalues of reordered pair (A, B) and
806 * normalize the generalized Schur form.
807 *
808  pair = .false.
809  DO 70 k = 1, n
810  IF( pair ) THEN
811  pair = .false.
812  ELSE
813 *
814  IF( k.LT.n ) THEN
815  IF( a( k+1, k ).NE.zero ) THEN
816  pair = .true.
817  END IF
818  END IF
819 *
820  IF( pair ) THEN
821 *
822 * Compute the eigenvalue(s) at position K.
823 *
824  work( 1 ) = a( k, k )
825  work( 2 ) = a( k+1, k )
826  work( 3 ) = a( k, k+1 )
827  work( 4 ) = a( k+1, k+1 )
828  work( 5 ) = b( k, k )
829  work( 6 ) = b( k+1, k )
830  work( 7 ) = b( k, k+1 )
831  work( 8 ) = b( k+1, k+1 )
832  CALL slag2( work, 2, work( 5 ), 2, smlnum*eps, beta( k ),
833  $ beta( k+1 ), alphar( k ), alphar( k+1 ),
834  $ alphai( k ) )
835  alphai( k+1 ) = -alphai( k )
836 *
837  ELSE
838 *
839  IF( sign( one, b( k, k ) ).LT.zero ) THEN
840 *
841 * If B(K,K) is negative, make it positive
842 *
843  DO 80 i = 1, n
844  a( k, i ) = -a( k, i )
845  b( k, i ) = -b( k, i )
846  IF( wantq ) q( i, k ) = -q( i, k )
847  80 CONTINUE
848  END IF
849 *
850  alphar( k ) = a( k, k )
851  alphai( k ) = zero
852  beta( k ) = b( k, k )
853 *
854  END IF
855  END IF
856  70 CONTINUE
857 *
858  work( 1 ) = lwmin
859  iwork( 1 ) = liwmin
860 *
861  RETURN
862 *
863 * End of STGSEN
864 *
865  END
subroutine stgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
STGSEN
Definition: stgsen.f:453
subroutine slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.
Definition: slassq.f:105
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine stgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
STGSYL
Definition: stgsyl.f:301
subroutine slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: slacn2.f:138
subroutine slag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition: slag2.f:158
subroutine stgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO)
STGEXC
Definition: stgexc.f:222