1: *> \brief \b DTGSEN
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DTGSEN + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsen.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsen.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsen.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DTGSEN( 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: * DOUBLE PRECISION PL, PR
30: * ..
31: * .. Array Arguments ..
32: * LOGICAL SELECT( * )
33: * INTEGER IWORK( * )
34: * DOUBLE PRECISION 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: *> DTGSEN 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 DGGES), 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: *> DTGSEN 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, DTGSEN 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
156: *> \endverbatim
157: *>
158: *> \param[out] ALPHAI
159: *> \verbatim
160: *> ALPHAI is DOUBLE PRECISION array, dimension (N)
161: *> \endverbatim
162: *>
163: *> \param[out] BETA
164: *> \verbatim
165: *> BETA is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
224: *> \endverbatim
225: *>
226: *> \param[out] PR
227: *> \verbatim
228: *> PR is DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleOTHERcomputational
308: *
309: *> \par Further Details:
310: * =====================
311: *>
312: *> \verbatim
313: *>
314: *> DTGSEN 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 DLATDF), then the parameter
411: *> IDIFJB (see below) should be changed from 3 to 4 (routine DLATDF
412: *> (IJOB = 2 will be used)). See DTGSYL 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 DTGSEN( 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: DOUBLE PRECISION PL, PR
461: * ..
462: * .. Array Arguments ..
463: LOGICAL SELECT( * )
464: INTEGER IWORK( * )
465: DOUBLE PRECISION 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: DOUBLE PRECISION ZERO, ONE
476: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+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: DOUBLE PRECISION DSCALE, DSUM, EPS, RDSCAL, SMLNUM
484: * ..
485: * .. Local Arrays ..
486: INTEGER ISAVE( 3 )
487: * ..
488: * .. External Subroutines ..
489: EXTERNAL DLACN2, DLACPY, DLAG2, DLASSQ, DTGEXC, DTGSYL,
490: $ XERBLA
491: * ..
492: * .. External Functions ..
493: DOUBLE PRECISION DLAMCH
494: EXTERNAL DLAMCH
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( 'DTGSEN', -INFO )
522: RETURN
523: END IF
524: *
525: * Get machine constants
526: *
527: EPS = DLAMCH( 'P' )
528: SMLNUM = DLAMCH( '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( 'DTGSEN', -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 DLASSQ( N, A( 1, I ), 1, DSCALE, DSUM )
602: CALL DLASSQ( 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 DTGEXC( 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 DLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
671: CALL DLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
672: $ N1 )
673: CALL DTGSYL( '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 DLASSQ( 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 DLASSQ( 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 DTGSYL( '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 DTGSYL( '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 DLACN2. 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 DLACN2( 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 DTGSYL( '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 DTGSYL( '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 DLACN2( 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 DTGSYL( '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 DTGSYL( '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 80 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 DLAG2( 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 70 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: 70 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: 80 CONTINUE
854: *
855: WORK( 1 ) = LWMIN
856: IWORK( 1 ) = LIWMIN
857: *
858: RETURN
859: *
860: * End of DTGSEN
861: *
862: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtgsen.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack