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