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