ctgsen()

subroutine ctgsen	(	integer	ijob,
		logical	wantq,
		logical	wantz,
		logical, dimension( * )	select,
		integer	n,
		complex, dimension( lda, * )	a,
		integer	lda,
		complex, dimension( ldb, * )	b,
		integer	ldb,
		complex, dimension( * )	alpha,
		complex, dimension( * )	beta,
		complex, dimension( ldq, * )	q,
		integer	ldq,
		complex, dimension( ldz, * )	z,
		integer	ldz,
		integer	m,
		real	pl,
		real	pr,
		real, dimension( * )	dif,
		complex, dimension( * )	work,
		integer	lwork,
		integer, dimension( * )	iwork,
		integer	liwork,
		integer	info )

CTGSEN

Download CTGSEN + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!> CTGSEN reorders the generalized Schur decomposition of a complex
!> matrix pair (A, B) (in terms of an unitary equivalence trans-
!> formation Q**H * (A, B) * Z), so that a selected cluster of eigenvalues
!> appears in the leading diagonal blocks of the pair (A,B). The leading
!> columns of Q and Z form unitary bases of the corresponding left and
!> right eigenspaces (deflating subspaces). (A, B) must be in
!> generalized Schur canonical form, that is, A and B are both upper
!> triangular.
!>
!> CTGSEN also computes the generalized eigenvalues
!>
!>          w(j)= ALPHA(j) / BETA(j)
!>
!> of the reordered matrix pair (A, B).
!>
!> Optionally, the routine computes estimates of reciprocal condition
!> numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
!> (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
!> between the matrix pairs (A11, B11) and (A22,B22) that correspond to
!> the selected cluster and the eigenvalues outside the cluster, resp.,
!> and norms of  onto left and right eigenspaces w.r.t.
!> the selected cluster in the (1,1)-block.
!>
!> 

Parameters

[in]	IJOB	!> IJOB is INTEGER !> Specifies whether condition numbers are required for the !> cluster of eigenvalues (PL and PR) or the deflating subspaces !> (Difu and Difl): !> =0: Only reorder w.r.t. SELECT. No extras. !> =1: Reciprocal of norms of onto left and right !> eigenspaces w.r.t. the selected cluster (PL and PR). !> =2: Upper bounds on Difu and Difl. F-norm-based estimate !> (DIF(1:2)). !> =3: Estimate of Difu and Difl. 1-norm-based estimate !> (DIF(1:2)). !> About 5 times as expensive as IJOB = 2. !> =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic !> version to get it all. !> =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above) !>
[in]	WANTQ	!> WANTQ is LOGICAL !> .TRUE. : update the left transformation matrix Q; !> .FALSE.: do not update Q. !>
[in]	WANTZ	!> WANTZ is LOGICAL !> .TRUE. : update the right transformation matrix Z; !> .FALSE.: do not update Z. !>
[in]	SELECT	!> SELECT is LOGICAL array, dimension (N) !> SELECT specifies the eigenvalues in the selected cluster. To !> select an eigenvalue w(j), SELECT(j) must be set to !> .TRUE.. !>
[in]	N	!> N is INTEGER !> The order of the matrices A and B. N >= 0. !>
[in,out]	A	!> A is COMPLEX array, dimension(LDA,N) !> On entry, the upper triangular matrix A, in generalized !> Schur canonical form. !> On exit, A is overwritten by the reordered matrix A. !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,N). !>
[in,out]	B	!> B is COMPLEX array, dimension(LDB,N) !> On entry, the upper triangular matrix B, in generalized !> Schur canonical form. !> On exit, B is overwritten by the reordered matrix B. !>
[in]	LDB	!> LDB is INTEGER !> The leading dimension of the array B. LDB >= max(1,N). !>
[out]	ALPHA	!> ALPHA is COMPLEX array, dimension (N) !>
[out]	BETA	!> BETA is COMPLEX array, dimension (N) !> !> The diagonal elements of A and B, respectively, !> when the pair (A,B) has been reduced to generalized Schur !> form. ALPHA(i)/BETA(i) i=1,...,N are the generalized !> eigenvalues. !>
[in,out]	Q	!> Q is COMPLEX array, dimension (LDQ,N) !> On entry, if WANTQ = .TRUE., Q is an N-by-N matrix. !> On exit, Q has been postmultiplied by the left unitary !> transformation matrix which reorder (A, B); The leading M !> columns of Q form orthonormal bases for the specified pair of !> left eigenspaces (deflating subspaces). !> If WANTQ = .FALSE., Q is not referenced. !>
[in]	LDQ	!> LDQ is INTEGER !> The leading dimension of the array Q. LDQ >= 1. !> If WANTQ = .TRUE., LDQ >= N. !>
[in,out]	Z	!> Z is COMPLEX array, dimension (LDZ,N) !> On entry, if WANTZ = .TRUE., Z is an N-by-N matrix. !> On exit, Z has been postmultiplied by the left unitary !> transformation matrix which reorder (A, B); The leading M !> columns of Z form orthonormal bases for the specified pair of !> left eigenspaces (deflating subspaces). !> If WANTZ = .FALSE., Z is not referenced. !>
[in]	LDZ	!> LDZ is INTEGER !> The leading dimension of the array Z. LDZ >= 1. !> If WANTZ = .TRUE., LDZ >= N. !>
[out]	M	!> M is INTEGER !> The dimension of the specified pair of left and right !> eigenspaces, (deflating subspaces) 0 <= M <= N. !>
[out]	PL	!> PL is REAL !>
[out]	PR	!> PR is REAL !> !> If IJOB = 1, 4 or 5, PL, PR are lower bounds on the !> reciprocal of the norm of onto left and right !> eigenspace with respect to the selected cluster. !> 0 < PL, PR <= 1. !> If M = 0 or M = N, PL = PR = 1. !> If IJOB = 0, 2 or 3 PL, PR are not referenced. !>
[out]	DIF	!> DIF is REAL array, dimension (2). !> If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl. !> If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on !> Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based !> estimates of Difu and Difl, computed using reversed !> communication with CLACN2. !> If M = 0 or N, DIF(1:2) = F-norm([A, B]). !> If IJOB = 0 or 1, DIF is not referenced. !>
[out]	WORK	!> WORK is COMPLEX array, dimension (MAX(1,LWORK)) !> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. !>
[in]	LWORK	!> LWORK is INTEGER !> The dimension of the array WORK. LWORK >= 1 !> If IJOB = 1, 2 or 4, LWORK >= 2*M*(N-M) !> If IJOB = 3 or 5, LWORK >= 4*M*(N-M) !> !> If LWORK = -1, then a workspace query is assumed; the routine !> only calculates the optimal size of the WORK array, returns !> this value as the first entry of the WORK array, and no error !> message related to LWORK is issued by XERBLA. !>
[out]	IWORK	!> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) !> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. !>
[in]	LIWORK	!> LIWORK is INTEGER !> The dimension of the array IWORK. LIWORK >= 1. !> If IJOB = 1, 2 or 4, LIWORK >= N+2; !> If IJOB = 3 or 5, LIWORK >= MAX(N+2, 2*M*(N-M)); !> !> If LIWORK = -1, then a workspace query is assumed; the !> routine only calculates the optimal size of the IWORK array, !> returns this value as the first entry of the IWORK array, and !> no error message related to LIWORK is issued by XERBLA. !>
[out]	INFO	!> INFO is INTEGER !> =0: Successful exit. !> <0: If INFO = -i, the i-th argument had an illegal value. !> =1: Reordering of (A, B) failed because the transformed !> matrix pair (A, B) would be too far from generalized !> Schur form; the problem is very ill-conditioned. !> (A, B) may have been partially reordered. !> If requested, 0 is returned in DIF(*), PL and PR. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Further Details:

!>
!>  CTGSEN first collects the selected eigenvalues by computing unitary
!>  U and W that move them to the top left corner of (A, B). In other
!>  words, the selected eigenvalues are the eigenvalues of (A11, B11) in
!>
!>              U**H*(A, B)*W = (A11 A12) (B11 B12) n1
!>                              ( 0  A22),( 0  B22) n2
!>                                n1  n2    n1  n2
!>
!>  where N = n1+n2 and U**H means the conjugate transpose of U. The first
!>  n1 columns of U and W span the specified pair of left and right
!>  eigenspaces (deflating subspaces) of (A, B).
!>
!>  If (A, B) has been obtained from the generalized real Schur
!>  decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then the
!>  reordered generalized Schur form of (C, D) is given by
!>
!>           (C, D) = (Q*U)*(U**H *(A, B)*W)*(Z*W)**H,
!>
!>  and the first n1 columns of Q*U and Z*W span the corresponding
!>  deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
!>
!>  Note that if the selected eigenvalue is sufficiently ill-conditioned,
!>  then its value may differ significantly from its value before
!>  reordering.
!>
!>  The reciprocal condition numbers of the left and right eigenspaces
!>  spanned by the first n1 columns of U and W (or Q*U and Z*W) may
!>  be returned in DIF(1:2), corresponding to Difu and Difl, resp.
!>
!>  The Difu and Difl are defined as:
!>
!>       Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
!>  and
!>       Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
!>
!>  where sigma-min(Zu) is the smallest singular value of the
!>  (2*n1*n2)-by-(2*n1*n2) matrix
!>
!>       Zu = [ kron(In2, A11)  -kron(A22**H, In1) ]
!>            [ kron(In2, B11)  -kron(B22**H, In1) ].
!>
!>  Here, Inx is the identity matrix of size nx and A22**H is the
!>  conjugate transpose of A22. kron(X, Y) is the Kronecker product between
!>  the matrices X and Y.
!>
!>  When DIF(2) is small, small changes in (A, B) can cause large changes
!>  in the deflating subspace. An approximate (asymptotic) bound on the
!>  maximum angular error in the computed deflating subspaces is
!>
!>       EPS * norm((A, B)) / DIF(2),
!>
!>  where EPS is the machine precision.
!>
!>  The reciprocal norm of the projectors on the left and right
!>  eigenspaces associated with (A11, B11) may be returned in PL and PR.
!>  They are computed as follows. First we compute L and R so that
!>  P*(A, B)*Q is block diagonal, where
!>
!>       P = ( I -L ) n1           Q = ( I R ) n1
!>           ( 0  I ) n2    and        ( 0 I ) n2
!>             n1 n2                    n1 n2
!>
!>  and (L, R) is the solution to the generalized Sylvester equation
!>
!>       A11*R - L*A22 = -A12
!>       B11*R - L*B22 = -B12
!>
!>  Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).
!>  An approximate (asymptotic) bound on the average absolute error of
!>  the selected eigenvalues is
!>
!>       EPS * norm((A, B)) / PL.
!>
!>  There are also global error bounds which valid for perturbations up
!>  to a certain restriction:  A lower bound (x) on the smallest
!>  F-norm(E,F) for which an eigenvalue of (A11, B11) may move and
!>  coalesce with an eigenvalue of (A22, B22) under perturbation (E,F),
!>  (i.e. (A + E, B + F), is
!>
!>   x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
!>
!>  An approximate bound on x can be computed from DIF(1:2), PL and PR.
!>
!>  If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed
!>  (L', R') and unperturbed (L, R) left and right deflating subspaces
!>  associated with the selected cluster in the (1,1)-blocks can be
!>  bounded as
!>
!>   max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
!>   max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
!>
!>  See LAPACK User's Guide section 4.11 or the following references
!>  for more information.
!>
!>  Note that if the default method for computing the Frobenius-norm-
!>  based estimate DIF is not wanted (see CLATDF), then the parameter
!>  IDIFJB (see below) should be changed from 3 to 4 (routine CLATDF
!>  (IJOB = 2 will be used)). See CTGSYL for more details.
!> 

Contributors:

Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

References:

[1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the Generalized Real Schur Form of a Regular Matrix Pair (A, B), in M.S. Moonen et al (eds), Linear Algebra for Large Scale and Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
[2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified Eigenvalues of a Regular Matrix Pair (A, B) and Condition Estimation: Theory, Algorithms and Software, Report UMINF - 94.04, Department of Computing Science, Umea University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87. To appear in Numerical Algorithms, 1996.
[3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software for Solving the Generalized Sylvester Equation and Estimating the Separation between Regular Matrix Pairs, Report UMINF - 93.23, Department of Computing Science, Umea University, S-901 87 Umea, Sweden, December 1993, Revised April 1994, Also as LAPACK working Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.

Definition at line 428 of file ctgsen.f.

*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      LOGICAL            WANTQ, WANTZ
      INTEGER            IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK,
     $                   M, N
      REAL               PL, PR
*     ..
*     .. Array Arguments ..
      LOGICAL            SELECT( * )
      INTEGER            IWORK( * )
      REAL               DIF( * )
      COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
     $                   BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            IDIFJB
      parameter( idifjb = 3 )
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            LQUERY, SWAP, WANTD, WANTD1, WANTD2, WANTP
      INTEGER            I, IERR, IJB, K, KASE, KS, LIWMIN, LWMIN, MN2,
     $                   N1, N2
      REAL               DSCALE, DSUM, RDSCAL, SAFMIN
      COMPLEX            TEMP1, TEMP2
*     ..
*     .. Local Arrays ..
      INTEGER            ISAVE( 3 )
*     ..
*     .. External Functions ..
      REAL               SROUNDUP_LWORK
      EXTERNAL           sroundup_lwork
*     ..
*     .. External Subroutines ..
      REAL               SLAMCH
      EXTERNAL           clacn2, clacpy, classq, cscal, ctgexc,
     $                   ctgsyl,
     $                   slamch, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, cmplx, conjg, max, sqrt
*     ..
*     .. Executable Statements ..
*
*     Decode and test the input parameters
*
      info = 0
      lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
*
      IF( ijob.LT.0 .OR. ijob.GT.5 ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -5
      ELSE IF( lda.LT.max( 1, n ) ) THEN
         info = -7
      ELSE IF( ldb.LT.max( 1, n ) ) THEN
         info = -9
      ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
         info = -13
      ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
         info = -15
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CTGSEN', -info )
         RETURN
      END IF
*
      ierr = 0
*
      wantp = ijob.EQ.1 .OR. ijob.GE.4
      wantd1 = ijob.EQ.2 .OR. ijob.EQ.4
      wantd2 = ijob.EQ.3 .OR. ijob.EQ.5
      wantd = wantd1 .OR. wantd2
*
*     Set M to the dimension of the specified pair of deflating
*     subspaces.
*
      m = 0
      IF( .NOT.lquery .OR. ijob.NE.0 ) THEN
      DO 10 k = 1, n
         alpha( k ) = a( k, k )
         beta( k ) = b( k, k )
         IF( k.LT.n ) THEN
            IF( SELECT( k ) )
     $         m = m + 1
         ELSE
            IF( SELECT( n ) )
     $         m = m + 1
         END IF
   10 CONTINUE
      END IF
*
      IF( ijob.EQ.1 .OR. ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
         lwmin = max( 1, 2*m*(n-m) )
         liwmin = max( 1, n+2 )
      ELSE IF( ijob.EQ.3 .OR. ijob.EQ.5 ) THEN
         lwmin = max( 1, 4*m*(n-m) )
         liwmin = max( 1, 2*m*(n-m), n+2 )
      ELSE
         lwmin = 1
         liwmin = 1
      END IF
*
      work( 1 ) =  sroundup_lwork(lwmin)
      iwork( 1 ) = liwmin
*
      IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
         info = -21
      ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
         info = -23
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CTGSEN', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
*     Quick return if possible.
*
      IF( m.EQ.n .OR. m.EQ.0 ) THEN
         IF( wantp ) THEN
            pl = one
            pr = one
         END IF
         IF( wantd ) THEN
            dscale = zero
            dsum = one
            DO 20 i = 1, n
               CALL classq( n, a( 1, i ), 1, dscale, dsum )
               CALL classq( n, b( 1, i ), 1, dscale, dsum )
   20       CONTINUE
            dif( 1 ) = dscale*sqrt( dsum )
            dif( 2 ) = dif( 1 )
         END IF
         GO TO 70
      END IF
*
*     Get machine constant
*
      safmin = slamch( 'S' )
*
*     Collect the selected blocks at the top-left corner of (A, B).
*
      ks = 0
      DO 30 k = 1, n
         swap = SELECT( k )
         IF( swap ) THEN
            ks = ks + 1
*
*           Swap the K-th block to position KS. Compute unitary Q
*           and Z that will swap adjacent diagonal blocks in (A, B).
*
            IF( k.NE.ks )
     $         CALL ctgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq,
     $                      z,
     $                      ldz, k, ks, ierr )
*
            IF( ierr.GT.0 ) THEN
*
*              Swap is rejected: exit.
*
               info = 1
               IF( wantp ) THEN
                  pl = zero
                  pr = zero
               END IF
               IF( wantd ) THEN
                  dif( 1 ) = zero
                  dif( 2 ) = zero
               END IF
               GO TO 70
            END IF
         END IF
   30 CONTINUE
      IF( wantp ) THEN
*
*        Solve generalized Sylvester equation for R and L:
*                   A11 * R - L * A22 = A12
*                   B11 * R - L * B22 = B12
*
         n1 = m
         n2 = n - m
         i = n1 + 1
         CALL clacpy( 'Full', n1, n2, a( 1, i ), lda, work, n1 )
         CALL clacpy( 'Full', n1, n2, b( 1, i ), ldb,
     $                work( n1*n2+1 ),
     $                n1 )
         ijb = 0
         CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,
     $                n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ), n1,
     $                dscale, dif( 1 ), work( n1*n2*2+1 ),
     $                lwork-2*n1*n2, iwork, ierr )
*
*        Estimate the reciprocal of norms of "projections" onto
*        left and right eigenspaces
*
         rdscal = zero
         dsum = one
         CALL classq( n1*n2, work, 1, rdscal, dsum )
         pl = rdscal*sqrt( dsum )
         IF( pl.EQ.zero ) THEN
            pl = one
         ELSE
            pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
         END IF
         rdscal = zero
         dsum = one
         CALL classq( n1*n2, work( n1*n2+1 ), 1, rdscal, dsum )
         pr = rdscal*sqrt( dsum )
         IF( pr.EQ.zero ) THEN
            pr = one
         ELSE
            pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
         END IF
      END IF
      IF( wantd ) THEN
*
*        Compute estimates Difu and Difl.
*
         IF( wantd1 ) THEN
            n1 = m
            n2 = n - m
            i = n1 + 1
            ijb = idifjb
*
*           Frobenius norm-based Difu estimate.
*
            CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,
     $                   work,
     $                   n1, b, ldb, b( i, i ), ldb, work( n1*n2+1 ),
     $                   n1, dscale, dif( 1 ), work( n1*n2*2+1 ),
     $                   lwork-2*n1*n2, iwork, ierr )
*
*           Frobenius norm-based Difl estimate.
*
            CALL ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,
     $                   work,
     $                   n2, b( i, i ), ldb, b, ldb, work( n1*n2+1 ),
     $                   n2, dscale, dif( 2 ), work( n1*n2*2+1 ),
     $                   lwork-2*n1*n2, iwork, ierr )
         ELSE
*
*           Compute 1-norm-based estimates of Difu and Difl using
*           reversed communication with CLACN2. In each step a
*           generalized Sylvester equation or a transposed variant
*           is solved.
*
            kase = 0
            n1 = m
            n2 = n - m
            i = n1 + 1
            ijb = 0
            mn2 = 2*n1*n2
*
*           1-norm-based estimate of Difu.
*
   40       CONTINUE
            CALL clacn2( mn2, work( mn2+1 ), work, dif( 1 ), kase,
     $                   isave )
            IF( kase.NE.0 ) THEN
               IF( kase.EQ.1 ) THEN
*
*                 Solve generalized Sylvester equation
*
                  CALL ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ),
     $                         lda,
     $                         work, n1, b, ldb, b( i, i ), ldb,
     $                         work( n1*n2+1 ), n1, dscale, dif( 1 ),
     $                         work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
     $                         ierr )
               ELSE
*
*                 Solve the transposed variant.
*
                  CALL ctgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ),
     $                         lda,
     $                         work, n1, b, ldb, b( i, i ), ldb,
     $                         work( n1*n2+1 ), n1, dscale, dif( 1 ),
     $                         work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
     $                         ierr )
               END IF
               GO TO 40
            END IF
            dif( 1 ) = dscale / dif( 1 )
*
*           1-norm-based estimate of Difl.
*
   50       CONTINUE
            CALL clacn2( mn2, work( mn2+1 ), work, dif( 2 ), kase,
     $                   isave )
            IF( kase.NE.0 ) THEN
               IF( kase.EQ.1 ) THEN
*
*                 Solve generalized Sylvester equation
*
                  CALL ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a,
     $                         lda,
     $                         work, n2, b( i, i ), ldb, b, ldb,
     $                         work( n1*n2+1 ), n2, dscale, dif( 2 ),
     $                         work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
     $                         ierr )
               ELSE
*
*                 Solve the transposed variant.
*
                  CALL ctgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a,
     $                         lda,
     $                         work, n2, b, ldb, b( i, i ), ldb,
     $                         work( n1*n2+1 ), n2, dscale, dif( 2 ),
     $                         work( n1*n2*2+1 ), lwork-2*n1*n2, iwork,
     $                         ierr )
               END IF
               GO TO 50
            END IF
            dif( 2 ) = dscale / dif( 2 )
         END IF
      END IF
*
*     If B(K,K) is complex, make it real and positive (normalization
*     of the generalized Schur form) and Store the generalized
*     eigenvalues of reordered pair (A, B)
*
      DO 60 k = 1, n
         dscale = abs( b( k, k ) )
         IF( dscale.GT.safmin ) THEN
            temp1 = conjg( b( k, k ) / dscale )
            temp2 = b( k, k ) / dscale
            b( k, k ) = dscale
            CALL cscal( n-k, temp1, b( k, k+1 ), ldb )
            CALL cscal( n-k+1, temp1, a( k, k ), lda )
            IF( wantq )
     $         CALL cscal( n, temp2, q( 1, k ), 1 )
         ELSE
            b( k, k ) = cmplx( zero, zero )
         END IF
*
         alpha( k ) = a( k, k )
         beta( k ) = b( k, k )
*
   60 CONTINUE
*
   70 CONTINUE
*
      work( 1 ) =  sroundup_lwork(lwmin)
      iwork( 1 ) = liwmin
*
      RETURN
*
*     End of CTGSEN
*

Here is the call graph for this function:

Here is the caller graph for this function: