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