001:       SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB,
002:      $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL,
003:      $                   PR, DIF, WORK, LWORK, IWORK, LIWORK, 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: *     January 2007
009: *
010: *     Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
011: *
012: *     .. Scalar Arguments ..
013:       LOGICAL            WANTQ, WANTZ
014:       INTEGER            IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
015:      $                   M, N
016:       DOUBLE PRECISION   PL, PR
017: *     ..
018: *     .. Array Arguments ..
019:       LOGICAL            SELECT( * )
020:       INTEGER            IWORK( * )
021:       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
022:      $                   B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ),
023:      $                   WORK( * ), Z( LDZ, * )
024: *     ..
025: *
026: *  Purpose
027: *  =======
028: *
029: *  DTGSEN reorders the generalized real Schur decomposition of a real
030: *  matrix pair (A, B) (in terms of an orthonormal equivalence trans-
031: *  formation Q' * (A, B) * Z), so that a selected cluster of eigenvalues
032: *  appears in the leading diagonal blocks of the upper quasi-triangular
033: *  matrix A and the upper triangular B. The leading columns of Q and
034: *  Z form orthonormal bases of the corresponding left and right eigen-
035: *  spaces (deflating subspaces). (A, B) must be in generalized real
036: *  Schur canonical form (as returned by DGGES), i.e. A is block upper
037: *  triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper
038: *  triangular.
039: *
040: *  DTGSEN also computes the generalized eigenvalues
041: *
042: *              w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)
043: *
044: *  of the reordered matrix pair (A, B).
045: *
046: *  Optionally, DTGSEN computes the estimates of reciprocal condition
047: *  numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
048: *  (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
049: *  between the matrix pairs (A11, B11) and (A22,B22) that correspond to
050: *  the selected cluster and the eigenvalues outside the cluster, resp.,
051: *  and norms of "projections" onto left and right eigenspaces w.r.t.
052: *  the selected cluster in the (1,1)-block.
053: *
054: *  Arguments
055: *  =========
056: *
057: *  IJOB    (input) INTEGER
058: *          Specifies whether condition numbers are required for the
059: *          cluster of eigenvalues (PL and PR) or the deflating subspaces
060: *          (Difu and Difl):
061: *           =0: Only reorder w.r.t. SELECT. No extras.
062: *           =1: Reciprocal of norms of "projections" onto left and right
063: *               eigenspaces w.r.t. the selected cluster (PL and PR).
064: *           =2: Upper bounds on Difu and Difl. F-norm-based estimate
065: *               (DIF(1:2)).
066: *           =3: Estimate of Difu and Difl. 1-norm-based estimate
067: *               (DIF(1:2)).
068: *               About 5 times as expensive as IJOB = 2.
069: *           =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic
070: *               version to get it all.
071: *           =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above)
072: *
073: *  WANTQ   (input) LOGICAL
074: *          .TRUE. : update the left transformation matrix Q;
075: *          .FALSE.: do not update Q.
076: *
077: *  WANTZ   (input) LOGICAL
078: *          .TRUE. : update the right transformation matrix Z;
079: *          .FALSE.: do not update Z.
080: *
081: *  SELECT  (input) LOGICAL array, dimension (N)
082: *          SELECT specifies the eigenvalues in the selected cluster.
083: *          To select a real eigenvalue w(j), SELECT(j) must be set to
084: *          .TRUE.. To select a complex conjugate pair of eigenvalues
085: *          w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
086: *          either SELECT(j) or SELECT(j+1) or both must be set to
087: *          .TRUE.; a complex conjugate pair of eigenvalues must be
088: *          either both included in the cluster or both excluded.
089: *
090: *  N       (input) INTEGER
091: *          The order of the matrices A and B. N >= 0.
092: *
093: *  A       (input/output) DOUBLE PRECISION array, dimension(LDA,N)
094: *          On entry, the upper quasi-triangular matrix A, with (A, B) in
095: *          generalized real Schur canonical form.
096: *          On exit, A is overwritten by the reordered matrix A.
097: *
098: *  LDA     (input) INTEGER
099: *          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))
182: *          IF IJOB = 0, IWORK is not referenced.  Otherwise,
183: *          on exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
184: *
185: *  LIWORK  (input) INTEGER
186: *          The dimension of the array IWORK. LIWORK >= 1.
187: *          If IJOB = 1, 2 or 4, LIWORK >=  N+6.
188: *          If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-M), N+6).
189: *
190: *          If LIWORK = -1, then a workspace query is assumed; the
191: *          routine only calculates the optimal size of the IWORK array,
192: *          returns this value as the first entry of the IWORK array, and
193: *          no error message related to LIWORK is issued by XERBLA.
194: *
195: *  INFO    (output) INTEGER
196: *            =0: Successful exit.
197: *            <0: If INFO = -i, the i-th argument had an illegal value.
198: *            =1: Reordering of (A, B) failed because the transformed
199: *                matrix pair (A, B) would be too far from generalized
200: *                Schur form; the problem is very ill-conditioned.
201: *                (A, B) may have been partially reordered.
202: *                If requested, 0 is returned in DIF(*), PL and PR.
203: *
204: *  Further Details
205: *  ===============
206: *
207: *  DTGSEN first collects the selected eigenvalues by computing
208: *  orthogonal U and W that move them to the top left corner of (A, B).
209: *  In other words, the selected eigenvalues are the eigenvalues of
210: *  (A11, B11) in:
211: *
212: *                U'*(A, B)*W = (A11 A12) (B11 B12) n1
213: *                              ( 0  A22),( 0  B22) n2
214: *                                n1  n2    n1  n2
215: *
216: *  where N = n1+n2 and U' means the transpose of U. The first n1 columns
217: *  of U and W span the specified pair of left and right eigenspaces
218: *  (deflating subspaces) of (A, B).
219: *
220: *  If (A, B) has been obtained from the generalized real Schur
221: *  decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then the
222: *  reordered generalized real Schur form of (C, D) is given by
223: *
224: *           (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',
225: *
226: *  and the first n1 columns of Q*U and Z*W span the corresponding
227: *  deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
228: *
229: *  Note that if the selected eigenvalue is sufficiently ill-conditioned,
230: *  then its value may differ significantly from its value before
231: *  reordering.
232: *
233: *  The reciprocal condition numbers of the left and right eigenspaces
234: *  spanned by the first n1 columns of U and W (or Q*U and Z*W) may
235: *  be returned in DIF(1:2), corresponding to Difu and Difl, resp.
236: *
237: *  The Difu and Difl are defined as:
238: *
239: *       Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
240: *  and
241: *       Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
242: *
243: *  where sigma-min(Zu) is the smallest singular value of the
244: *  (2*n1*n2)-by-(2*n1*n2) matrix
245: *
246: *       Zu = [ kron(In2, A11)  -kron(A22', In1) ]
247: *            [ kron(In2, B11)  -kron(B22', In1) ].
248: *
249: *  Here, Inx is the identity matrix of size nx and A22' is the
250: *  transpose of A22. kron(X, Y) is the Kronecker product between
251: *  the matrices X and Y.
252: *
253: *  When DIF(2) is small, small changes in (A, B) can cause large changes
254: *  in the deflating subspace. An approximate (asymptotic) bound on the
255: *  maximum angular error in the computed deflating subspaces is
256: *
257: *       EPS * norm((A, B)) / DIF(2),
258: *
259: *  where EPS is the machine precision.
260: *
261: *  The reciprocal norm of the projectors on the left and right
262: *  eigenspaces associated with (A11, B11) may be returned in PL and PR.
263: *  They are computed as follows. First we compute L and R so that
264: *  P*(A, B)*Q is block diagonal, where
265: *
266: *       P = ( I -L ) n1           Q = ( I R ) n1
267: *           ( 0  I ) n2    and        ( 0 I ) n2
268: *             n1 n2                    n1 n2
269: *
270: *  and (L, R) is the solution to the generalized Sylvester equation
271: *
272: *       A11*R - L*A22 = -A12
273: *       B11*R - L*B22 = -B12
274: *
275: *  Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
276: *  An approximate (asymptotic) bound on the average absolute error of
277: *  the selected eigenvalues is
278: *
279: *       EPS * norm((A, B)) / PL.
280: *
281: *  There are also global error bounds which valid for perturbations up
282: *  to a certain restriction:  A lower bound (x) on the smallest
283: *  F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
284: *  coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
285: *  (i.e. (A + E, B + F), is
286: *
287: *   x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
288: *
289: *  An approximate bound on x can be computed from DIF(1:2), PL and PR.
290: *
291: *  If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
292: *  (L', R') and unperturbed (L, R) left and right deflating subspaces
293: *  associated with the selected cluster in the (1,1)-blocks can be
294: *  bounded as
295: *
296: *   max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
297: *   max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
298: *
299: *  See LAPACK User's Guide section 4.11 or the following references
300: *  for more information.
301: *
302: *  Note that if the default method for computing the Frobenius-norm-
303: *  based estimate DIF is not wanted (see DLATDF), then the parameter
304: *  IDIFJB (see below) should be changed from 3 to 4 (routine DLATDF
305: *  (IJOB = 2 will be used)). See DTGSYL for more details.
306: *
307: *  Based on contributions by
308: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
309: *     Umea University, S-901 87 Umea, Sweden.
310: *
311: *  References
312: *  ==========
313: *
314: *  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
315: *      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
316: *      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
317: *      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
318: *
319: *  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
320: *      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
321: *      Estimation: Theory, Algorithms and Software,
322: *      Report UMINF - 94.04, Department of Computing Science, Umea
323: *      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
324: *      Note 87. To appear in Numerical Algorithms, 1996.
325: *
326: *  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
327: *      for Solving the Generalized Sylvester Equation and Estimating the
328: *      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
329: *      Department of Computing Science, Umea University, S-901 87 Umea,
330: *      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
331: *      Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
332: *      1996.
333: *
334: *  =====================================================================
335: *
336: *     .. Parameters ..
337:       INTEGER            IDIFJB
338:       PARAMETER          ( IDIFJB = 3 )
339:       DOUBLE PRECISION   ZERO, ONE
340:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
341: *     ..
342: *     .. Local Scalars ..
343:       LOGICAL            LQUERY, PAIR, SWAP, WANTD, WANTD1, WANTD2,
344:      $                   WANTP
345:       INTEGER            I, IERR, IJB, K, KASE, KK, KS, LIWMIN, LWMIN,
346:      $                   MN2, N1, N2
347:       DOUBLE PRECISION   DSCALE, DSUM, EPS, RDSCAL, SMLNUM
348: *     ..
349: *     .. Local Arrays ..
350:       INTEGER            ISAVE( 3 )
351: *     ..
352: *     .. External Subroutines ..
353:       EXTERNAL           DLACN2, DLACPY, DLAG2, DLASSQ, DTGEXC, DTGSYL,
354:      $                   XERBLA
355: *     ..
356: *     .. External Functions ..
357:       DOUBLE PRECISION   DLAMCH
358:       EXTERNAL           DLAMCH
359: *     ..
360: *     .. Intrinsic Functions ..
361:       INTRINSIC          MAX, SIGN, SQRT
362: *     ..
363: *     .. Executable Statements ..
364: *
365: *     Decode and test the input parameters
366: *
367:       INFO = 0
368:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
369: *
370:       IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN
371:          INFO = -1
372:       ELSE IF( N.LT.0 ) THEN
373:          INFO = -5
374:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
375:          INFO = -7
376:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
377:          INFO = -9
378:       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
379:          INFO = -14
380:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
381:          INFO = -16
382:       END IF
383: *
384:       IF( INFO.NE.0 ) THEN
385:          CALL XERBLA( 'DTGSEN', -INFO )
386:          RETURN
387:       END IF
388: *
389: *     Get machine constants
390: *
391:       EPS = DLAMCH( 'P' )
392:       SMLNUM = DLAMCH( 'S' ) / EPS
393:       IERR = 0
394: *
395:       WANTP = IJOB.EQ.1 .OR. IJOB.GE.4
396:       WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4
397:       WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5
398:       WANTD = WANTD1 .OR. WANTD2
399: *
400: *     Set M to the dimension of the specified pair of deflating
401: *     subspaces.
402: *
403:       M = 0
404:       PAIR = .FALSE.
405:       DO 10 K = 1, N
406:          IF( PAIR ) THEN
407:             PAIR = .FALSE.
408:          ELSE
409:             IF( K.LT.N ) THEN
410:                IF( A( K+1, K ).EQ.ZERO ) THEN
411:                   IF( SELECT( K ) )
412:      $               M = M + 1
413:                ELSE
414:                   PAIR = .TRUE.
415:                   IF( SELECT( K ) .OR. SELECT( K+1 ) )
416:      $               M = M + 2
417:                END IF
418:             ELSE
419:                IF( SELECT( N ) )
420:      $            M = M + 1
421:             END IF
422:          END IF
423:    10 CONTINUE
424: *
425:       IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
426:          LWMIN = MAX( 1, 4*N+16, 2*M*( N-M ) )
427:          LIWMIN = MAX( 1, N+6 )
428:       ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN
429:          LWMIN = MAX( 1, 4*N+16, 4*M*( N-M ) )
430:          LIWMIN = MAX( 1, 2*M*( N-M ), N+6 )
431:       ELSE
432:          LWMIN = MAX( 1, 4*N+16 )
433:          LIWMIN = 1
434:       END IF
435: *
436:       WORK( 1 ) = LWMIN
437:       IWORK( 1 ) = LIWMIN
438: *
439:       IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
440:          INFO = -22
441:       ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
442:          INFO = -24
443:       END IF
444: *
445:       IF( INFO.NE.0 ) THEN
446:          CALL XERBLA( 'DTGSEN', -INFO )
447:          RETURN
448:       ELSE IF( LQUERY ) THEN
449:          RETURN
450:       END IF
451: *
452: *     Quick return if possible.
453: *
454:       IF( M.EQ.N .OR. M.EQ.0 ) THEN
455:          IF( WANTP ) THEN
456:             PL = ONE
457:             PR = ONE
458:          END IF
459:          IF( WANTD ) THEN
460:             DSCALE = ZERO
461:             DSUM = ONE
462:             DO 20 I = 1, N
463:                CALL DLASSQ( N, A( 1, I ), 1, DSCALE, DSUM )
464:                CALL DLASSQ( N, B( 1, I ), 1, DSCALE, DSUM )
465:    20       CONTINUE
466:             DIF( 1 ) = DSCALE*SQRT( DSUM )
467:             DIF( 2 ) = DIF( 1 )
468:          END IF
469:          GO TO 60
470:       END IF
471: *
472: *     Collect the selected blocks at the top-left corner of (A, B).
473: *
474:       KS = 0
475:       PAIR = .FALSE.
476:       DO 30 K = 1, N
477:          IF( PAIR ) THEN
478:             PAIR = .FALSE.
479:          ELSE
480: *
481:             SWAP = SELECT( K )
482:             IF( K.LT.N ) THEN
483:                IF( A( K+1, K ).NE.ZERO ) THEN
484:                   PAIR = .TRUE.
485:                   SWAP = SWAP .OR. SELECT( K+1 )
486:                END IF
487:             END IF
488: *
489:             IF( SWAP ) THEN
490:                KS = KS + 1
491: *
492: *              Swap the K-th block to position KS.
493: *              Perform the reordering of diagonal blocks in (A, B)
494: *              by orthogonal transformation matrices and update
495: *              Q and Z accordingly (if requested):
496: *
497:                KK = K
498:                IF( K.NE.KS )
499:      $            CALL DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
500:      $                         Z, LDZ, KK, KS, WORK, LWORK, IERR )
501: *
502:                IF( IERR.GT.0 ) THEN
503: *
504: *                 Swap is rejected: exit.
505: *
506:                   INFO = 1
507:                   IF( WANTP ) THEN
508:                      PL = ZERO
509:                      PR = ZERO
510:                   END IF
511:                   IF( WANTD ) THEN
512:                      DIF( 1 ) = ZERO
513:                      DIF( 2 ) = ZERO
514:                   END IF
515:                   GO TO 60
516:                END IF
517: *
518:                IF( PAIR )
519:      $            KS = KS + 1
520:             END IF
521:          END IF
522:    30 CONTINUE
523:       IF( WANTP ) THEN
524: *
525: *        Solve generalized Sylvester equation for R and L
526: *        and compute PL and PR.
527: *
528:          N1 = M
529:          N2 = N - M
530:          I = N1 + 1
531:          IJB = 0
532:          CALL DLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 )
533:          CALL DLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ),
534:      $                N1 )
535:          CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
536:      $                N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1,
537:      $                DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ),
538:      $                LWORK-2*N1*N2, IWORK, IERR )
539: *
540: *        Estimate the reciprocal of norms of "projections" onto left
541: *        and right eigenspaces.
542: *
543:          RDSCAL = ZERO
544:          DSUM = ONE
545:          CALL DLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM )
546:          PL = RDSCAL*SQRT( DSUM )
547:          IF( PL.EQ.ZERO ) THEN
548:             PL = ONE
549:          ELSE
550:             PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) )
551:          END IF
552:          RDSCAL = ZERO
553:          DSUM = ONE
554:          CALL DLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM )
555:          PR = RDSCAL*SQRT( DSUM )
556:          IF( PR.EQ.ZERO ) THEN
557:             PR = ONE
558:          ELSE
559:             PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) )
560:          END IF
561:       END IF
562: *
563:       IF( WANTD ) THEN
564: *
565: *        Compute estimates of Difu and Difl.
566: *
567:          IF( WANTD1 ) THEN
568:             N1 = M
569:             N2 = N - M
570:             I = N1 + 1
571:             IJB = IDIFJB
572: *
573: *           Frobenius norm-based Difu-estimate.
574: *
575:             CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK,
576:      $                   N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ),
577:      $                   N1, DSCALE, DIF( 1 ), WORK( 2*N1*N2+1 ),
578:      $                   LWORK-2*N1*N2, IWORK, IERR )
579: *
580: *           Frobenius norm-based Difl-estimate.
581: *
582:             CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK,
583:      $                   N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ),
584:      $                   N2, DSCALE, DIF( 2 ), WORK( 2*N1*N2+1 ),
585:      $                   LWORK-2*N1*N2, IWORK, IERR )
586:          ELSE
587: *
588: *
589: *           Compute 1-norm-based estimates of Difu and Difl using
590: *           reversed communication with DLACN2. In each step a
591: *           generalized Sylvester equation or a transposed variant
592: *           is solved.
593: *
594:             KASE = 0
595:             N1 = M
596:             N2 = N - M
597:             I = N1 + 1
598:             IJB = 0
599:             MN2 = 2*N1*N2
600: *
601: *           1-norm-based estimate of Difu.
602: *
603:    40       CONTINUE
604:             CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 1 ),
605:      $                   KASE, ISAVE )
606:             IF( KASE.NE.0 ) THEN
607:                IF( KASE.EQ.1 ) THEN
608: *
609: *                 Solve generalized Sylvester equation.
610: *
611:                   CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA,
612:      $                         WORK, N1, B, LDB, B( I, I ), LDB,
613:      $                         WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
614:      $                         WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
615:      $                         IERR )
616:                ELSE
617: *
618: *                 Solve the transposed variant.
619: *
620:                   CALL DTGSYL( 'T', IJB, N1, N2, A, LDA, A( I, I ), LDA,
621:      $                         WORK, N1, B, LDB, B( I, I ), LDB,
622:      $                         WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ),
623:      $                         WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
624:      $                         IERR )
625:                END IF
626:                GO TO 40
627:             END IF
628:             DIF( 1 ) = DSCALE / DIF( 1 )
629: *
630: *           1-norm-based estimate of Difl.
631: *
632:    50       CONTINUE
633:             CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 2 ),
634:      $                   KASE, ISAVE )
635:             IF( KASE.NE.0 ) THEN
636:                IF( KASE.EQ.1 ) THEN
637: *
638: *                 Solve generalized Sylvester equation.
639: *
640:                   CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA,
641:      $                         WORK, N2, B( I, I ), LDB, B, LDB,
642:      $                         WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
643:      $                         WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
644:      $                         IERR )
645:                ELSE
646: *
647: *                 Solve the transposed variant.
648: *
649:                   CALL DTGSYL( 'T', IJB, N2, N1, A( I, I ), LDA, A, LDA,
650:      $                         WORK, N2, B( I, I ), LDB, B, LDB,
651:      $                         WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ),
652:      $                         WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK,
653:      $                         IERR )
654:                END IF
655:                GO TO 50
656:             END IF
657:             DIF( 2 ) = DSCALE / DIF( 2 )
658: *
659:          END IF
660:       END IF
661: *
662:    60 CONTINUE
663: *
664: *     Compute generalized eigenvalues of reordered pair (A, B) and
665: *     normalize the generalized Schur form.
666: *
667:       PAIR = .FALSE.
668:       DO 80 K = 1, N
669:          IF( PAIR ) THEN
670:             PAIR = .FALSE.
671:          ELSE
672: *
673:             IF( K.LT.N ) THEN
674:                IF( A( K+1, K ).NE.ZERO ) THEN
675:                   PAIR = .TRUE.
676:                END IF
677:             END IF
678: *
679:             IF( PAIR ) THEN
680: *
681: *             Compute the eigenvalue(s) at position K.
682: *
683:                WORK( 1 ) = A( K, K )
684:                WORK( 2 ) = A( K+1, K )
685:                WORK( 3 ) = A( K, K+1 )
686:                WORK( 4 ) = A( K+1, K+1 )
687:                WORK( 5 ) = B( K, K )
688:                WORK( 6 ) = B( K+1, K )
689:                WORK( 7 ) = B( K, K+1 )
690:                WORK( 8 ) = B( K+1, K+1 )
691:                CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA( K ),
692:      $                     BETA( K+1 ), ALPHAR( K ), ALPHAR( K+1 ),
693:      $                     ALPHAI( K ) )
694:                ALPHAI( K+1 ) = -ALPHAI( K )
695: *
696:             ELSE
697: *
698:                IF( SIGN( ONE, B( K, K ) ).LT.ZERO ) THEN
699: *
700: *                 If B(K,K) is negative, make it positive
701: *
702:                   DO 70 I = 1, N
703:                      A( K, I ) = -A( K, I )
704:                      B( K, I ) = -B( K, I )
705:                      IF( WANTQ ) Q( I, K ) = -Q( I, K )
706:    70             CONTINUE
707:                END IF
708: *
709:                ALPHAR( K ) = A( K, K )
710:                ALPHAI( K ) = ZERO
711:                BETA( K ) = B( K, K )
712: *
713:             END IF
714:          END IF
715:    80 CONTINUE
716: *
717:       WORK( 1 ) = LWMIN
718:       IWORK( 1 ) = LIWMIN
719: *
720:       RETURN
721: *
722: *     End of DTGSEN
723: *
724:       END
725: