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: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       CHARACTER          HOWMNY, JOB
012:       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
013: *     ..
014: *     .. Array Arguments ..
015:       LOGICAL            SELECT( * )
016:       INTEGER            IWORK( * )
017:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), DIF( * ), S( * ),
018:      $                   VL( LDVL, * ), VR( LDVR, * ), WORK( * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  DTGSNA estimates reciprocal condition numbers for specified
025: *  eigenvalues and/or eigenvectors of a matrix pair (A, B) in
026: *  generalized real Schur canonical form (or of any matrix pair
027: *  (Q*A*Z', Q*B*Z') with orthogonal matrices Q and Z, where
028: *  Z' denotes the transpose of Z.
029: *
030: *  (A, B) must be in generalized real Schur form (as returned by DGGES),
031: *  i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
032: *  blocks. B is upper triangular.
033: *
034: *
035: *  Arguments
036: *  =========
037: *
038: *  JOB     (input) CHARACTER*1
039: *          Specifies whether condition numbers are required for
040: *          eigenvalues (S) or eigenvectors (DIF):
041: *          = 'E': for eigenvalues only (S);
042: *          = 'V': for eigenvectors only (DIF);
043: *          = 'B': for both eigenvalues and eigenvectors (S and DIF).
044: *
045: *  HOWMNY  (input) CHARACTER*1
046: *          = 'A': compute condition numbers for all eigenpairs;
047: *          = 'S': compute condition numbers for selected eigenpairs
048: *                 specified by the array SELECT.
049: *
050: *  SELECT  (input) LOGICAL array, dimension (N)
051: *          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
052: *          condition numbers are required. To select condition numbers
053: *          for the eigenpair corresponding to a real eigenvalue w(j),
054: *          SELECT(j) must be set to .TRUE.. To select condition numbers
055: *          corresponding to a complex conjugate pair of eigenvalues w(j)
056: *          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
057: *          set to .TRUE..
058: *          If HOWMNY = 'A', SELECT is not referenced.
059: *
060: *  N       (input) INTEGER
061: *          The order of the square matrix pair (A, B). N >= 0.
062: *
063: *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
064: *          The upper quasi-triangular matrix A in the pair (A,B).
065: *
066: *  LDA     (input) INTEGER
067: *          The leading dimension of the array A. LDA >= max(1,N).
068: *
069: *  B       (input) DOUBLE PRECISION array, dimension (LDB,N)
070: *          The upper triangular matrix B in the pair (A,B).
071: *
072: *  LDB     (input) INTEGER
073: *          The leading dimension of the array B. LDB >= max(1,N).
074: *
075: *  VL      (input) DOUBLE PRECISION array, dimension (LDVL,M)
076: *          If JOB = 'E' or 'B', VL must contain left eigenvectors of
077: *          (A, B), corresponding to the eigenpairs specified by HOWMNY
078: *          and SELECT. The eigenvectors must be stored in consecutive
079: *          columns of VL, as returned by DTGEVC.
080: *          If JOB = 'V', VL is not referenced.
081: *
082: *  LDVL    (input) INTEGER
083: *          The leading dimension of the array VL. LDVL >= 1.
084: *          If JOB = 'E' or 'B', LDVL >= N.
085: *
086: *  VR      (input) DOUBLE PRECISION array, dimension (LDVR,M)
087: *          If JOB = 'E' or 'B', VR must contain right eigenvectors of
088: *          (A, B), corresponding to the eigenpairs specified by HOWMNY
089: *          and SELECT. The eigenvectors must be stored in consecutive
090: *          columns ov VR, as returned by DTGEVC.
091: *          If JOB = 'V', VR is not referenced.
092: *
093: *  LDVR    (input) INTEGER
094: *          The leading dimension of the array VR. LDVR >= 1.
095: *          If JOB = 'E' or 'B', LDVR >= N.
096: *
097: *  S       (output) DOUBLE PRECISION array, dimension (MM)
098: *          If JOB = 'E' or 'B', the reciprocal condition numbers of the
099: *          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
582: