1: SUBROUTINE DTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
2: $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
3: $ IWORK, INFO )
4: *
5: * -- LAPACK routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * .. Scalar Arguments ..
11: CHARACTER HOWMNY, JOB
12: INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
13: * ..
14: * .. Array Arguments ..
15: LOGICAL SELECT( * )
16: INTEGER IWORK( * )
17: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
18: $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * DTGSNA estimates reciprocal condition numbers for specified
25: * eigenvalues and/or eigenvectors of a matrix pair (A, B) in
26: * generalized real Schur canonical form (or of any matrix pair
27: * (Q*A*Z', Q*B*Z') with orthogonal matrices Q and Z, where
28: * Z' denotes the transpose of Z.
29: *
30: * (A, B) must be in generalized real Schur form (as returned by DGGES),
31: * i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
32: * blocks. B is upper triangular.
33: *
34: *
35: * Arguments
36: * =========
37: *
38: * JOB (input) CHARACTER*1
39: * Specifies whether condition numbers are required for
40: * eigenvalues (S) or eigenvectors (DIF):
41: * = 'E': for eigenvalues only (S);
42: * = 'V': for eigenvectors only (DIF);
43: * = 'B': for both eigenvalues and eigenvectors (S and DIF).
44: *
45: * HOWMNY (input) CHARACTER*1
46: * = 'A': compute condition numbers for all eigenpairs;
47: * = 'S': compute condition numbers for selected eigenpairs
48: * specified by the array SELECT.
49: *
50: * SELECT (input) LOGICAL array, dimension (N)
51: * If HOWMNY = 'S', SELECT specifies the eigenpairs for which
52: * condition numbers are required. To select condition numbers
53: * for the eigenpair corresponding to a real eigenvalue w(j),
54: * SELECT(j) must be set to .TRUE.. To select condition numbers
55: * corresponding to a complex conjugate pair of eigenvalues w(j)
56: * and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
57: * set to .TRUE..
58: * If HOWMNY = 'A', SELECT is not referenced.
59: *
60: * N (input) INTEGER
61: * The order of the square matrix pair (A, B). N >= 0.
62: *
63: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
64: * The upper quasi-triangular matrix A in the pair (A,B).
65: *
66: * LDA (input) INTEGER
67: * The leading dimension of the array A. LDA >= max(1,N).
68: *
69: * B (input) DOUBLE PRECISION array, dimension (LDB,N)
70: * The upper triangular matrix B in the pair (A,B).
71: *
72: * LDB (input) INTEGER
73: * The leading dimension of the array B. LDB >= max(1,N).
74: *
75: * VL (input) DOUBLE PRECISION array, dimension (LDVL,M)
76: * If JOB = 'E' or 'B', VL must contain left eigenvectors of
77: * (A, B), corresponding to the eigenpairs specified by HOWMNY
78: * and SELECT. The eigenvectors must be stored in consecutive
79: * columns of VL, as returned by DTGEVC.
80: * If JOB = 'V', VL is not referenced.
81: *
82: * LDVL (input) INTEGER
83: * The leading dimension of the array VL. LDVL >= 1.
84: * If JOB = 'E' or 'B', LDVL >= N.
85: *
86: * VR (input) DOUBLE PRECISION array, dimension (LDVR,M)
87: * If JOB = 'E' or 'B', VR must contain right eigenvectors of
88: * (A, B), corresponding to the eigenpairs specified by HOWMNY
89: * and SELECT. The eigenvectors must be stored in consecutive
90: * columns ov VR, as returned by DTGEVC.
91: * If JOB = 'V', VR is not referenced.
92: *
93: * LDVR (input) INTEGER
94: * The leading dimension of the array VR. LDVR >= 1.
95: * If JOB = 'E' or 'B', LDVR >= N.
96: *
97: * S (output) DOUBLE PRECISION array, dimension (MM)
98: * If JOB = 'E' or 'B', the reciprocal condition numbers of the
99: * selected eigenvalues, stored in consecutive elements of the
100: * array. For a complex conjugate pair of eigenvalues two
101: * consecutive elements of S are set to the same value. Thus
102: * S(j), DIF(j), and the j-th columns of VL and VR all
103: * correspond to the same eigenpair (but not in general the
104: * j-th eigenpair, unless all eigenpairs are selected).
105: * If JOB = 'V', S is not referenced.
106: *
107: * DIF (output) DOUBLE PRECISION array, dimension (MM)
108: * If JOB = 'V' or 'B', the estimated reciprocal condition
109: * numbers of the selected eigenvectors, stored in consecutive
110: * elements of the array. For a complex eigenvector two
111: * consecutive elements of DIF are set to the same value. If
112: * the eigenvalues cannot be reordered to compute DIF(j), DIF(j)
113: * is set to 0; this can only occur when the true value would be
114: * very small anyway.
115: * If JOB = 'E', DIF is not referenced.
116: *
117: * MM (input) INTEGER
118: * The number of elements in the arrays S and DIF. MM >= M.
119: *
120: * M (output) INTEGER
121: * The number of elements of the arrays S and DIF used to store
122: * the specified condition numbers; for each selected real
123: * eigenvalue one element is used, and for each selected complex
124: * conjugate pair of eigenvalues, two elements are used.
125: * If HOWMNY = 'A', M is set to N.
126: *
127: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
128: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
129: *
130: * LWORK (input) INTEGER
131: * The dimension of the array WORK. LWORK >= max(1,N).
132: * If JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16.
133: *
134: * If LWORK = -1, then a workspace query is assumed; the routine
135: * only calculates the optimal size of the WORK array, returns
136: * this value as the first entry of the WORK array, and no error
137: * message related to LWORK is issued by XERBLA.
138: *
139: * IWORK (workspace) INTEGER array, dimension (N + 6)
140: * If JOB = 'E', IWORK is not referenced.
141: *
142: * INFO (output) INTEGER
143: * =0: Successful exit
144: * <0: If INFO = -i, the i-th argument had an illegal value
145: *
146: *
147: * Further Details
148: * ===============
149: *
150: * The reciprocal of the condition number of a generalized eigenvalue
151: * w = (a, b) is defined as
152: *
153: * S(w) = (|u'Av|**2 + |u'Bv|**2)**(1/2) / (norm(u)*norm(v))
154: *
155: * where u and v are the left and right eigenvectors of (A, B)
156: * corresponding to w; |z| denotes the absolute value of the complex
157: * number, and norm(u) denotes the 2-norm of the vector u.
158: * The pair (a, b) corresponds to an eigenvalue w = a/b (= u'Av/u'Bv)
159: * of the matrix pair (A, B). If both a and b equal zero, then (A B) is
160: * singular and S(I) = -1 is returned.
161: *
162: * An approximate error bound on the chordal distance between the i-th
163: * computed generalized eigenvalue w and the corresponding exact
164: * eigenvalue lambda is
165: *
166: * chord(w, lambda) <= EPS * norm(A, B) / S(I)
167: *
168: * where EPS is the machine precision.
169: *
170: * The reciprocal of the condition number DIF(i) of right eigenvector u
171: * and left eigenvector v corresponding to the generalized eigenvalue w
172: * is defined as follows:
173: *
174: * a) If the i-th eigenvalue w = (a,b) is real
175: *
176: * Suppose U and V are orthogonal transformations such that
177: *
178: * U'*(A, B)*V = (S, T) = ( a * ) ( b * ) 1
179: * ( 0 S22 ),( 0 T22 ) n-1
180: * 1 n-1 1 n-1
181: *
182: * Then the reciprocal condition number DIF(i) is
183: *
184: * Difl((a, b), (S22, T22)) = sigma-min( Zl ),
185: *
186: * where sigma-min(Zl) denotes the smallest singular value of the
187: * 2(n-1)-by-2(n-1) matrix
188: *
189: * Zl = [ kron(a, In-1) -kron(1, S22) ]
190: * [ kron(b, In-1) -kron(1, T22) ] .
191: *
192: * Here In-1 is the identity matrix of size n-1. kron(X, Y) is the
193: * Kronecker product between the matrices X and Y.
194: *
195: * Note that if the default method for computing DIF(i) is wanted
196: * (see DLATDF), then the parameter DIFDRI (see below) should be
197: * changed from 3 to 4 (routine DLATDF(IJOB = 2 will be used)).
198: * See DTGSYL for more details.
199: *
200: * b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair,
201: *
202: * Suppose U and V are orthogonal transformations such that
203: *
204: * U'*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2
205: * ( 0 S22 ),( 0 T22) n-2
206: * 2 n-2 2 n-2
207: *
208: * and (S11, T11) corresponds to the complex conjugate eigenvalue
209: * pair (w, conjg(w)). There exist unitary matrices U1 and V1 such
210: * that
211: *
212: * U1'*S11*V1 = ( s11 s12 ) and U1'*T11*V1 = ( t11 t12 )
213: * ( 0 s22 ) ( 0 t22 )
214: *
215: * where the generalized eigenvalues w = s11/t11 and
216: * conjg(w) = s22/t22.
217: *
218: * Then the reciprocal condition number DIF(i) is bounded by
219: *
220: * min( d1, max( 1, |real(s11)/real(s22)| )*d2 )
221: *
222: * where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where
223: * Z1 is the complex 2-by-2 matrix
224: *
225: * Z1 = [ s11 -s22 ]
226: * [ t11 -t22 ],
227: *
228: * This is done by computing (using real arithmetic) the
229: * roots of the characteristical polynomial det(Z1' * Z1 - lambda I),
230: * where Z1' denotes the conjugate transpose of Z1 and det(X) denotes
231: * the determinant of X.
232: *
233: * and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an
234: * upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2)
235: *
236: * Z2 = [ kron(S11', In-2) -kron(I2, S22) ]
237: * [ kron(T11', In-2) -kron(I2, T22) ]
238: *
239: * Note that if the default method for computing DIF is wanted (see
240: * DLATDF), then the parameter DIFDRI (see below) should be changed
241: * from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). See DTGSYL
242: * for more details.
243: *
244: * For each eigenvalue/vector specified by SELECT, DIF stores a
245: * Frobenius norm-based estimate of Difl.
246: *
247: * An approximate error bound for the i-th computed eigenvector VL(i) or
248: * VR(i) is given by
249: *
250: * EPS * norm(A, B) / DIF(i).
251: *
252: * See ref. [2-3] for more details and further references.
253: *
254: * Based on contributions by
255: * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
256: * Umea University, S-901 87 Umea, Sweden.
257: *
258: * References
259: * ==========
260: *
261: * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
262: * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
263: * M.S. Moonen et al (eds), Linear Algebra for Large Scale and
264: * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
265: *
266: * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
267: * Eigenvalues of a Regular Matrix Pair (A, B) and Condition
268: * Estimation: Theory, Algorithms and Software,
269: * Report UMINF - 94.04, Department of Computing Science, Umea
270: * University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
271: * Note 87. To appear in Numerical Algorithms, 1996.
272: *
273: * [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
274: * for Solving the Generalized Sylvester Equation and Estimating the
275: * Separation between Regular Matrix Pairs, Report UMINF - 93.23,
276: * Department of Computing Science, Umea University, S-901 87 Umea,
277: * Sweden, December 1993, Revised April 1994, Also as LAPACK Working
278: * Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
279: * No 1, 1996.
280: *
281: * =====================================================================
282: *
283: * .. Parameters ..
284: INTEGER DIFDRI
285: PARAMETER ( DIFDRI = 3 )
286: DOUBLE PRECISION ZERO, ONE, TWO, FOUR
287: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
288: $ FOUR = 4.0D+0 )
289: * ..
290: * .. Local Scalars ..
291: LOGICAL LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS
292: INTEGER I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2
293: DOUBLE PRECISION ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND,
294: $ EPS, LNRM, RNRM, ROOT1, ROOT2, SCALE, SMLNUM,
295: $ TMPII, TMPIR, TMPRI, TMPRR, UHAV, UHAVI, UHBV,
296: $ UHBVI
297: * ..
298: * .. Local Arrays ..
299: DOUBLE PRECISION DUMMY( 1 ), DUMMY1( 1 )
300: * ..
301: * .. External Functions ..
302: LOGICAL LSAME
303: DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
304: EXTERNAL LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
305: * ..
306: * .. External Subroutines ..
307: EXTERNAL DGEMV, DLACPY, DLAG2, DTGEXC, DTGSYL, XERBLA
308: * ..
309: * .. Intrinsic Functions ..
310: INTRINSIC MAX, MIN, SQRT
311: * ..
312: * .. Executable Statements ..
313: *
314: * Decode and test the input parameters
315: *
316: WANTBH = LSAME( JOB, 'B' )
317: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
318: WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
319: *
320: SOMCON = LSAME( HOWMNY, 'S' )
321: *
322: INFO = 0
323: LQUERY = ( LWORK.EQ.-1 )
324: *
325: IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
326: INFO = -1
327: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
328: INFO = -2
329: ELSE IF( N.LT.0 ) THEN
330: INFO = -4
331: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
332: INFO = -6
333: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
334: INFO = -8
335: ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
336: INFO = -10
337: ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
338: INFO = -12
339: ELSE
340: *
341: * Set M to the number of eigenpairs for which condition numbers
342: * are required, and test MM.
343: *
344: IF( SOMCON ) THEN
345: M = 0
346: PAIR = .FALSE.
347: DO 10 K = 1, N
348: IF( PAIR ) THEN
349: PAIR = .FALSE.
350: ELSE
351: IF( K.LT.N ) THEN
352: IF( A( K+1, K ).EQ.ZERO ) THEN
353: IF( SELECT( K ) )
354: $ M = M + 1
355: ELSE
356: PAIR = .TRUE.
357: IF( SELECT( K ) .OR. SELECT( K+1 ) )
358: $ M = M + 2
359: END IF
360: ELSE
361: IF( SELECT( N ) )
362: $ M = M + 1
363: END IF
364: END IF
365: 10 CONTINUE
366: ELSE
367: M = N
368: END IF
369: *
370: IF( N.EQ.0 ) THEN
371: LWMIN = 1
372: ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
373: LWMIN = 2*N*( N + 2 ) + 16
374: ELSE
375: LWMIN = N
376: END IF
377: WORK( 1 ) = LWMIN
378: *
379: IF( MM.LT.M ) THEN
380: INFO = -15
381: ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
382: INFO = -18
383: END IF
384: END IF
385: *
386: IF( INFO.NE.0 ) THEN
387: CALL XERBLA( 'DTGSNA', -INFO )
388: RETURN
389: ELSE IF( LQUERY ) THEN
390: RETURN
391: END IF
392: *
393: * Quick return if possible
394: *
395: IF( N.EQ.0 )
396: $ RETURN
397: *
398: * Get machine constants
399: *
400: EPS = DLAMCH( 'P' )
401: SMLNUM = DLAMCH( 'S' ) / EPS
402: KS = 0
403: PAIR = .FALSE.
404: *
405: DO 20 K = 1, N
406: *
407: * Determine whether A(k,k) begins a 1-by-1 or 2-by-2 block.
408: *
409: IF( PAIR ) THEN
410: PAIR = .FALSE.
411: GO TO 20
412: ELSE
413: IF( K.LT.N )
414: $ PAIR = A( K+1, K ).NE.ZERO
415: END IF
416: *
417: * Determine whether condition numbers are required for the k-th
418: * eigenpair.
419: *
420: IF( SOMCON ) THEN
421: IF( PAIR ) THEN
422: IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
423: $ GO TO 20
424: ELSE
425: IF( .NOT.SELECT( K ) )
426: $ GO TO 20
427: END IF
428: END IF
429: *
430: KS = KS + 1
431: *
432: IF( WANTS ) THEN
433: *
434: * Compute the reciprocal condition number of the k-th
435: * eigenvalue.
436: *
437: IF( PAIR ) THEN
438: *
439: * Complex eigenvalue pair.
440: *
441: RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
442: $ DNRM2( N, VR( 1, KS+1 ), 1 ) )
443: LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
444: $ DNRM2( N, VL( 1, KS+1 ), 1 ) )
445: CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
446: $ WORK, 1 )
447: TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
448: TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
449: CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS+1 ), 1,
450: $ ZERO, WORK, 1 )
451: TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
452: TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
453: UHAV = TMPRR + TMPII
454: UHAVI = TMPIR - TMPRI
455: CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
456: $ WORK, 1 )
457: TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
458: TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
459: CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS+1 ), 1,
460: $ ZERO, WORK, 1 )
461: TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 )
462: TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
463: UHBV = TMPRR + TMPII
464: UHBVI = TMPIR - TMPRI
465: UHAV = DLAPY2( UHAV, UHAVI )
466: UHBV = DLAPY2( UHBV, UHBVI )
467: COND = DLAPY2( UHAV, UHBV )
468: S( KS ) = COND / ( RNRM*LNRM )
469: S( KS+1 ) = S( KS )
470: *
471: ELSE
472: *
473: * Real eigenvalue.
474: *
475: RNRM = DNRM2( N, VR( 1, KS ), 1 )
476: LNRM = DNRM2( N, VL( 1, KS ), 1 )
477: CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO,
478: $ WORK, 1 )
479: UHAV = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
480: CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO,
481: $ WORK, 1 )
482: UHBV = DDOT( N, WORK, 1, VL( 1, KS ), 1 )
483: COND = DLAPY2( UHAV, UHBV )
484: IF( COND.EQ.ZERO ) THEN
485: S( KS ) = -ONE
486: ELSE
487: S( KS ) = COND / ( RNRM*LNRM )
488: END IF
489: END IF
490: END IF
491: *
492: IF( WANTDF ) THEN
493: IF( N.EQ.1 ) THEN
494: DIF( KS ) = DLAPY2( A( 1, 1 ), B( 1, 1 ) )
495: GO TO 20
496: END IF
497: *
498: * Estimate the reciprocal condition number of the k-th
499: * eigenvectors.
500: IF( PAIR ) THEN
501: *
502: * Copy the 2-by 2 pencil beginning at (A(k,k), B(k, k)).
503: * Compute the eigenvalue(s) at position K.
504: *
505: WORK( 1 ) = A( K, K )
506: WORK( 2 ) = A( K+1, K )
507: WORK( 3 ) = A( K, K+1 )
508: WORK( 4 ) = A( K+1, K+1 )
509: WORK( 5 ) = B( K, K )
510: WORK( 6 ) = B( K+1, K )
511: WORK( 7 ) = B( K, K+1 )
512: WORK( 8 ) = B( K+1, K+1 )
513: CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA,
514: $ DUMMY1( 1 ), ALPHAR, DUMMY( 1 ), ALPHAI )
515: ALPRQT = ONE
516: C1 = TWO*( ALPHAR*ALPHAR+ALPHAI*ALPHAI+BETA*BETA )
517: C2 = FOUR*BETA*BETA*ALPHAI*ALPHAI
518: ROOT1 = C1 + SQRT( C1*C1-4.0D0*C2 )
519: ROOT2 = C2 / ROOT1
520: ROOT1 = ROOT1 / TWO
521: COND = MIN( SQRT( ROOT1 ), SQRT( ROOT2 ) )
522: END IF
523: *
524: * Copy the matrix (A, B) to the array WORK and swap the
525: * diagonal block beginning at A(k,k) to the (1,1) position.
526: *
527: CALL DLACPY( 'Full', N, N, A, LDA, WORK, N )
528: CALL DLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
529: IFST = K
530: ILST = 1
531: *
532: CALL DTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), N,
533: $ DUMMY, 1, DUMMY1, 1, IFST, ILST,
534: $ WORK( N*N*2+1 ), LWORK-2*N*N, IERR )
535: *
536: IF( IERR.GT.0 ) THEN
537: *
538: * Ill-conditioned problem - swap rejected.
539: *
540: DIF( KS ) = ZERO
541: ELSE
542: *
543: * Reordering successful, solve generalized Sylvester
544: * equation for R and L,
545: * A22 * R - L * A11 = A12
546: * B22 * R - L * B11 = B12,
547: * and compute estimate of Difl((A11,B11), (A22, B22)).
548: *
549: N1 = 1
550: IF( WORK( 2 ).NE.ZERO )
551: $ N1 = 2
552: N2 = N - N1
553: IF( N2.EQ.0 ) THEN
554: DIF( KS ) = COND
555: ELSE
556: I = N*N + 1
557: IZ = 2*N*N + 1
558: CALL DTGSYL( 'N', DIFDRI, N2, N1, WORK( N*N1+N1+1 ),
559: $ N, WORK, N, WORK( N1+1 ), N,
560: $ WORK( N*N1+N1+I ), N, WORK( I ), N,
561: $ WORK( N1+I ), N, SCALE, DIF( KS ),
562: $ WORK( IZ+1 ), LWORK-2*N*N, IWORK, IERR )
563: *
564: IF( PAIR )
565: $ DIF( KS ) = MIN( MAX( ONE, ALPRQT )*DIF( KS ),
566: $ COND )
567: END IF
568: END IF
569: IF( PAIR )
570: $ DIF( KS+1 ) = DIF( KS )
571: END IF
572: IF( PAIR )
573: $ KS = KS + 1
574: *
575: 20 CONTINUE
576: WORK( 1 ) = LWMIN
577: RETURN
578: *
579: * End of DTGSNA
580: *
581: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtgsna.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack