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