Annotation of rpl/lapack/lapack/dtgsen.f, revision 1.8
1.1 bertrand 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: *
1.5 bertrand 5: * -- LAPACK routine (version 3.2.2) --
1.1 bertrand 6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * January 2007
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' * (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))
1.5 bertrand 182: * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
1.1 bertrand 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'*(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' 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', then the
221: * reordered generalized real Schur form of (C, D) is given by
222: *
223: * (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',
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', In1) ]
246: * [ kron(In2, B11) -kron(B22', In1) ].
247: *
248: * Here, Inx is the identity matrix of size nx and A22' 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
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