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