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