cdrgsx()

subroutine cdrgsx	(	integer	nsize,
		integer	ncmax,
		real	thresh,
		integer	nin,
		integer	nout,
		complex, dimension( lda, * )	a,
		integer	lda,
		complex, dimension( lda, * )	b,
		complex, dimension( lda, * )	ai,
		complex, dimension( lda, * )	bi,
		complex, dimension( lda, * )	z,
		complex, dimension( lda, * )	q,
		complex, dimension( * )	alpha,
		complex, dimension( * )	beta,
		complex, dimension( ldc, * )	c,
		integer	ldc,
		real, dimension( * )	s,
		complex, dimension( * )	work,
		integer	lwork,
		real, dimension( * )	rwork,
		integer, dimension( * )	iwork,
		integer	liwork,
		logical, dimension( * )	bwork,
		integer	info )

CDRGSX

Purpose:

!>
!> CDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
!> problem expert driver CGGESX.
!>
!> CGGES factors A and B as Q*S*Z'  and Q*T*Z' , where ' means conjugate
!> transpose, S and T are  upper triangular (i.e., in generalized Schur
!> form), and Q and Z are unitary. It also computes the generalized
!> eigenvalues (alpha(j),beta(j)), j=1,...,n.  Thus,
!> w(j) = alpha(j)/beta(j) is a root of the characteristic equation
!>
!>                 det( A - w(j) B ) = 0
!>
!> Optionally it also reorders the eigenvalues so that a selected
!> cluster of eigenvalues appears in the leading diagonal block of the
!> Schur forms; computes a reciprocal condition number for the average
!> of the selected eigenvalues; and computes a reciprocal condition
!> number for the right and left deflating subspaces corresponding to
!> the selected eigenvalues.
!>
!> When CDRGSX is called with NSIZE > 0, five (5) types of built-in
!> matrix pairs are used to test the routine CGGESX.
!>
!> When CDRGSX is called with NSIZE = 0, it reads in test matrix data
!> to test CGGESX.
!> (need more details on what kind of read-in data are needed).
!>
!> For each matrix pair, the following tests will be performed and
!> compared with the threshold THRESH except for the tests (7) and (9):
!>
!> (1)   | A - Q S Z' | / ( |A| n ulp )
!>
!> (2)   | B - Q T Z' | / ( |B| n ulp )
!>
!> (3)   | I - QQ' | / ( n ulp )
!>
!> (4)   | I - ZZ' | / ( n ulp )
!>
!> (5)   if A is in Schur form (i.e. triangular form)
!>
!> (6)   maximum over j of D(j)  where:
!>
!>                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
!>           D(j) = ------------------------ + -----------------------
!>                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
!>
!> (7)   if sorting worked and SDIM is the number of eigenvalues
!>       which were selected.
!>
!> (8)   the estimated value DIF does not differ from the true values of
!>       Difu and Difl more than a factor 10*THRESH. If the estimate DIF
!>       equals zero the corresponding true values of Difu and Difl
!>       should be less than EPS*norm(A, B). If the true value of Difu
!>       and Difl equal zero, the estimate DIF should be less than
!>       EPS*norm(A, B).
!>
!> (9)   If INFO = N+3 is returned by CGGESX, the reordering 
!>       and we check that DIF = PL = PR = 0 and that the true value of
!>       Difu and Difl is < EPS*norm(A, B). We count the events when
!>       INFO=N+3.
!>
!> For read-in test matrices, the same tests are run except that the
!> exact value for DIF (and PL) is input data.  Additionally, there is
!> one more test run for read-in test matrices:
!>
!> (10)  the estimated value PL does not differ from the true value of
!>       PLTRU more than a factor THRESH. If the estimate PL equals
!>       zero the corresponding true value of PLTRU should be less than
!>       EPS*norm(A, B). If the true value of PLTRU equal zero, the
!>       estimate PL should be less than EPS*norm(A, B).
!>
!> Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1)
!> matrix pairs are generated and tested. NSIZE should be kept small.
!>
!> SVD (routine CGESVD) is used for computing the true value of DIF_u
!> and DIF_l when testing the built-in test problems.
!>
!> Built-in Test Matrices
!> ======================
!>
!> All built-in test matrices are the 2 by 2 block of triangular
!> matrices
!>
!>          A = [ A11 A12 ]    and      B = [ B11 B12 ]
!>              [     A22 ]                 [     B22 ]
!>
!> where for different type of A11 and A22 are given as the following.
!> A12 and B12 are chosen so that the generalized Sylvester equation
!>
!>          A11*R - L*A22 = -A12
!>          B11*R - L*B22 = -B12
!>
!> have prescribed solution R and L.
!>
!> Type 1:  A11 = J_m(1,-1) and A_22 = J_k(1-a,1).
!>          B11 = I_m, B22 = I_k
!>          where J_k(a,b) is the k-by-k Jordan block with ``a'' on
!>          diagonal and ``b'' on superdiagonal.
!>
!> Type 2:  A11 = (a_ij) = ( 2(.5-sin(i)) ) and
!>          B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m
!>          A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and
!>          B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k
!>
!> Type 3:  A11, A22 and B11, B22 are chosen as for Type 2, but each
!>          second diagonal block in A_11 and each third diagonal block
!>          in A_22 are made as 2 by 2 blocks.
!>
!> Type 4:  A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) )
!>             for i=1,...,m,  j=1,...,m and
!>          A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) )
!>             for i=m+1,...,k,  j=m+1,...,k
!>
!> Type 5:  (A,B) and have potentially close or common eigenvalues and
!>          very large departure from block diagonality A_11 is chosen
!>          as the m x m leading submatrix of A_1:
!>                  |  1  b                            |
!>                  | -b  1                            |
!>                  |        1+d  b                    |
!>                  |         -b 1+d                   |
!>           A_1 =  |                  d  1            |
!>                  |                 -1  d            |
!>                  |                        -d  1     |
!>                  |                        -1 -d     |
!>                  |                               1  |
!>          and A_22 is chosen as the k x k leading submatrix of A_2:
!>                  | -1  b                            |
!>                  | -b -1                            |
!>                  |       1-d  b                     |
!>                  |       -b  1-d                    |
!>           A_2 =  |                 d 1+b            |
!>                  |               -1-b d             |
!>                  |                       -d  1+b    |
!>                  |                      -1+b  -d    |
!>                  |                              1-d |
!>          and matrix B are chosen as identity matrices (see SLATM5).
!>
!> 

Parameters

[in]	NSIZE	!> NSIZE is INTEGER !> The maximum size of the matrices to use. NSIZE >= 0. !> If NSIZE = 0, no built-in tests matrices are used, but !> read-in test matrices are used to test SGGESX. !>
[in]	NCMAX	!> NCMAX is INTEGER !> Maximum allowable NMAX for generating Kroneker matrix !> in call to CLAKF2 !>
[in]	THRESH	!> THRESH is REAL !> A test will count as if the , computed as !> described above, exceeds THRESH. Note that the error !> is scaled to be O(1), so THRESH should be a reasonably !> small multiple of 1, e.g., 10 or 100. In particular, !> it should not depend on the precision (single vs. double) !> or the size of the matrix. THRESH >= 0. !>
[in]	NIN	!> NIN is INTEGER !> The FORTRAN unit number for reading in the data file of !> problems to solve. !>
[in]	NOUT	!> NOUT is INTEGER !> The FORTRAN unit number for printing out error messages !> (e.g., if a routine returns INFO not equal to 0.) !>
[out]	A	!> A is COMPLEX array, dimension (LDA, NSIZE) !> Used to store the matrix whose eigenvalues are to be !> computed. On exit, A contains the last matrix actually used. !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of A, B, AI, BI, Z and Q, !> LDA >= max( 1, NSIZE ). For the read-in test, !> LDA >= max( 1, N ), N is the size of the test matrices. !>
[out]	B	!> B is COMPLEX array, dimension (LDA, NSIZE) !> Used to store the matrix whose eigenvalues are to be !> computed. On exit, B contains the last matrix actually used. !>
[out]	AI	!> AI is COMPLEX array, dimension (LDA, NSIZE) !> Copy of A, modified by CGGESX. !>
[out]	BI	!> BI is COMPLEX array, dimension (LDA, NSIZE) !> Copy of B, modified by CGGESX. !>
[out]	Z	!> Z is COMPLEX array, dimension (LDA, NSIZE) !> Z holds the left Schur vectors computed by CGGESX. !>
[out]	Q	!> Q is COMPLEX array, dimension (LDA, NSIZE) !> Q holds the right Schur vectors computed by CGGESX. !>
[out]	ALPHA	!> ALPHA is COMPLEX array, dimension (NSIZE) !>
[out]	BETA	!> BETA is COMPLEX array, dimension (NSIZE) !> !> On exit, ALPHA/BETA are the eigenvalues. !>
[out]	C	!> C is COMPLEX array, dimension (LDC, LDC) !> Store the matrix generated by subroutine CLAKF2, this is the !> matrix formed by Kronecker products used for estimating !> DIF. !>
[in]	LDC	!> LDC is INTEGER !> The leading dimension of C. LDC >= max(1, LDA*LDA/2 ). !>
[out]	S	!> S is REAL array, dimension (LDC) !> Singular values of C !>
[out]	WORK	!> WORK is COMPLEX array, dimension (LWORK) !>
[in]	LWORK	!> LWORK is INTEGER !> The dimension of the array WORK. LWORK >= 3*NSIZE*NSIZE/2 !>
[out]	RWORK	!> RWORK is REAL array, !> dimension (5*NSIZE*NSIZE/2 - 4) !>
[out]	IWORK	!> IWORK is INTEGER array, dimension (LIWORK) !>
[in]	LIWORK	!> LIWORK is INTEGER !> The dimension of the array IWORK. LIWORK >= NSIZE + 2. !>
[out]	BWORK	!> BWORK is LOGICAL array, dimension (NSIZE) !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -i, the i-th argument had an illegal value. !> > 0: A routine returned an error code. !>

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Definition at line 346 of file cdrgsx.f.

*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
     $                   NOUT, NSIZE
      REAL               THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            BWORK( * )
      INTEGER            IWORK( * )
      REAL               RWORK( * ), S( * )
      COMPLEX            A( LDA, * ), AI( LDA, * ), ALPHA( * ),
     $                   B( LDA, * ), BETA( * ), BI( LDA, * ),
     $                   C( LDC, * ), Q( LDA, * ), WORK( * ),
     $                   Z( LDA, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, TEN
      parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1 )
      COMPLEX            CZERO
      parameter( czero = ( 0.0e+0, 0.0e+0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            ILABAD
      CHARACTER          SENSE
      INTEGER            BDSPAC, I, IFUNC, J, LINFO, MAXWRK, MINWRK, MM,
     $                   MN2, NERRS, NPTKNT, NTEST, NTESTT, PRTYPE, QBA,
     $                   QBB
      REAL               ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
     $                   TEMP2, THRSH2, ULP, ULPINV, WEIGHT
      COMPLEX            X
*     ..
*     .. Local Arrays ..
      REAL               DIFEST( 2 ), PL( 2 ), RESULT( 10 )
*     ..
*     .. External Functions ..
      LOGICAL            CLCTSX
      INTEGER            ILAENV
      REAL               CLANGE, SLAMCH
      EXTERNAL           clctsx, ilaenv, clange, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           alasvm, cgesvd, cget51, cggesx, clacpy, clakf2,
     $                   claset, clatm5, xerbla
*     ..
*     .. Scalars in Common ..
      LOGICAL            FS
      INTEGER            K, M, MPLUSN, N
*     ..
*     .. Common blocks ..
      COMMON             / mn / m, n, mplusn, k, fs
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, aimag, max, real, sqrt
*     ..
*     .. Statement Functions ..
      REAL               ABS1
*     ..
*     .. Statement Function definitions ..
      abs1( x ) = abs( real( x ) ) + abs( aimag( x ) )
*     ..
*     .. Executable Statements ..
*
*     Check for errors
*
      IF( nsize.LT.0 ) THEN
         info = -1
      ELSE IF( thresh.LT.zero ) THEN
         info = -2
      ELSE IF( nin.LE.0 ) THEN
         info = -3
      ELSE IF( nout.LE.0 ) THEN
         info = -4
      ELSE IF( lda.LT.1 .OR. lda.LT.nsize ) THEN
         info = -6
      ELSE IF( ldc.LT.1 .OR. ldc.LT.nsize*nsize / 2 ) THEN
         info = -15
      ELSE IF( liwork.LT.nsize+2 ) THEN
         info = -21
      END IF
*
*     Compute workspace
*      (Note: Comments in the code beginning "Workspace:" describe the
*       minimal amount of workspace needed at that point in the code,
*       as well as the preferred amount for good performance.
*       NB refers to the optimal block size for the immediately
*       following subroutine, as returned by ILAENV.)
*
      minwrk = 1
      IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
         minwrk = 3*nsize*nsize / 2
*
*        workspace for cggesx
*
         maxwrk = nsize*( 1+ilaenv( 1, 'CGEQRF', ' ', nsize, 1, nsize,
     $            0 ) )
         maxwrk = max( maxwrk, nsize*( 1+ilaenv( 1, 'CUNGQR', ' ',
     $            nsize, 1, nsize, -1 ) ) )
*
*        workspace for cgesvd
*
         bdspac = 3*nsize*nsize / 2
         maxwrk = max( maxwrk, nsize*nsize*
     $            ( 1+ilaenv( 1, 'CGEBRD', ' ', nsize*nsize / 2,
     $            nsize*nsize / 2, -1, -1 ) ) )
         maxwrk = max( maxwrk, bdspac )
*
         maxwrk = max( maxwrk, minwrk )
*
         work( 1 ) = maxwrk
      END IF
*
      IF( lwork.LT.minwrk )
     $   info = -18
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CDRGSX', -info )
         RETURN
      END IF
*
*     Important constants
*
      ulp = slamch( 'P' )
      ulpinv = one / ulp
      smlnum = slamch( 'S' ) / ulp
      bignum = one / smlnum
      thrsh2 = ten*thresh
      ntestt = 0
      nerrs = 0
*
*     Go to the tests for read-in matrix pairs
*
      ifunc = 0
      IF( nsize.EQ.0 )
     $   GO TO 70
*
*     Test the built-in matrix pairs.
*     Loop over different functions (IFUNC) of CGGESX, types (PRTYPE)
*     of test matrices, different size (M+N)
*
      prtype = 0
      qba = 3
      qbb = 4
      weight = sqrt( ulp )
*
      DO 60 ifunc = 0, 3
         DO 50 prtype = 1, 5
            DO 40 m = 1, nsize - 1
               DO 30 n = 1, nsize - m
*
                  weight = one / weight
                  mplusn = m + n
*
*                 Generate test matrices
*
                  fs = .true.
                  k = 0
*
                  CALL claset( 'Full', mplusn, mplusn, czero, czero, ai,
     $                         lda )
                  CALL claset( 'Full', mplusn, mplusn, czero, czero, bi,
     $                         lda )
*
                  CALL clatm5( prtype, m, n, ai, lda, ai( m+1, m+1 ),
     $                         lda, ai( 1, m+1 ), lda, bi, lda,
     $                         bi( m+1, m+1 ), lda, bi( 1, m+1 ), lda,
     $                         q, lda, z, lda, weight, qba, qbb )
*
*                 Compute the Schur factorization and swapping the
*                 m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
*                 Swapping is accomplished via the function CLCTSX
*                 which is supplied below.
*
                  IF( ifunc.EQ.0 ) THEN
                     sense = 'N'
                  ELSE IF( ifunc.EQ.1 ) THEN
                     sense = 'E'
                  ELSE IF( ifunc.EQ.2 ) THEN
                     sense = 'V'
                  ELSE IF( ifunc.EQ.3 ) THEN
                     sense = 'B'
                  END IF
*
                  CALL clacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
                  CALL clacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
*
                  CALL cggesx( 'V', 'V', 'S', clctsx, sense, mplusn, ai,
     $                         lda, bi, lda, mm, alpha, beta, q, lda, z,
     $                         lda, pl, difest, work, lwork, rwork,
     $                         iwork, liwork, bwork, linfo )
*
                  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
                     result( 1 ) = ulpinv
                     WRITE( nout, fmt = 9999 )'CGGESX', linfo, mplusn,
     $                  prtype
                     info = linfo
                     GO TO 30
                  END IF
*
*                 Compute the norm(A, B)
*
                  CALL clacpy( 'Full', mplusn, mplusn, ai, lda, work,
     $                         mplusn )
                  CALL clacpy( 'Full', mplusn, mplusn, bi, lda,
     $                         work( mplusn*mplusn+1 ), mplusn )
                  abnrm = clange( 'Fro', mplusn, 2*mplusn, work, mplusn,
     $                    rwork )
*
*                 Do tests (1) to (4)
*
                  result( 2 ) = zero
                  CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
     $                         lda, work, rwork, result( 1 ) )
                  CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
     $                         lda, work, rwork, result( 2 ) )
                  CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
     $                         lda, work, rwork, result( 3 ) )
                  CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
     $                         lda, work, rwork, result( 4 ) )
                  ntest = 4
*
*                 Do tests (5) and (6): check Schur form of A and
*                 compare eigenvalues with diagonals.
*
                  temp1 = zero
                  result( 5 ) = zero
                  result( 6 ) = zero
*
                  DO 10 j = 1, mplusn
                     ilabad = .false.
                     temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
     $                       max( smlnum, abs1( alpha( j ) ),
     $                       abs1( ai( j, j ) ) )+
     $                       abs1( beta( j )-bi( j, j ) ) /
     $                       max( smlnum, abs1( beta( j ) ),
     $                       abs1( bi( j, j ) ) ) ) / ulp
                     IF( j.LT.mplusn ) THEN
                        IF( ai( j+1, j ).NE.zero ) THEN
                           ilabad = .true.
                           result( 5 ) = ulpinv
                        END IF
                     END IF
                     IF( j.GT.1 ) THEN
                        IF( ai( j, j-1 ).NE.zero ) THEN
                           ilabad = .true.
                           result( 5 ) = ulpinv
                        END IF
                     END IF
                     temp1 = max( temp1, temp2 )
                     IF( ilabad ) THEN
                        WRITE( nout, fmt = 9997 )j, mplusn, prtype
                     END IF
   10             CONTINUE
                  result( 6 ) = temp1
                  ntest = ntest + 2
*
*                 Test (7) (if sorting worked)
*
                  result( 7 ) = zero
                  IF( linfo.EQ.mplusn+3 ) THEN
                     result( 7 ) = ulpinv
                  ELSE IF( mm.NE.n ) THEN
                     result( 7 ) = ulpinv
                  END IF
                  ntest = ntest + 1
*
*                 Test (8): compare the estimated value DIF and its
*                 value. first, compute the exact DIF.
*
                  result( 8 ) = zero
                  mn2 = mm*( mplusn-mm )*2
                  IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
*
*                    Note: for either following two cases, there are
*                    almost same number of test cases fail the test.
*
                     CALL clakf2( mm, mplusn-mm, ai, lda,
     $                            ai( mm+1, mm+1 ), bi,
     $                            bi( mm+1, mm+1 ), c, ldc )
*
                     CALL cgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
     $                            1, work( 2 ), 1, work( 3 ), lwork-2,
     $                            rwork, info )
                     diftru = s( mn2 )
*
                     IF( difest( 2 ).EQ.zero ) THEN
                        IF( diftru.GT.abnrm*ulp )
     $                     result( 8 ) = ulpinv
                     ELSE IF( diftru.EQ.zero ) THEN
                        IF( difest( 2 ).GT.abnrm*ulp )
     $                     result( 8 ) = ulpinv
                     ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
     $                        ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
                        result( 8 ) = max( diftru / difest( 2 ),
     $                                difest( 2 ) / diftru )
                     END IF
                     ntest = ntest + 1
                  END IF
*
*                 Test (9)
*
                  result( 9 ) = zero
                  IF( linfo.EQ.( mplusn+2 ) ) THEN
                     IF( diftru.GT.abnrm*ulp )
     $                  result( 9 ) = ulpinv
                     IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
     $                  result( 9 ) = ulpinv
                     IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
     $                  result( 9 ) = ulpinv
                     ntest = ntest + 1
                  END IF
*
                  ntestt = ntestt + ntest
*
*                 Print out tests which fail.
*
                  DO 20 j = 1, 9
                     IF( result( j ).GE.thresh ) THEN
*
*                       If this is the first test to fail,
*                       print a header to the data file.
*
                        IF( nerrs.EQ.0 ) THEN
                           WRITE( nout, fmt = 9996 )'CGX'
*
*                          Matrix types
*
                           WRITE( nout, fmt = 9994 )
*
*                          Tests performed
*
                           WRITE( nout, fmt = 9993 )'unitary', '''',
     $                        'transpose', ( '''', i = 1, 4 )
*
                        END IF
                        nerrs = nerrs + 1
                        IF( result( j ).LT.10000.0 ) THEN
                           WRITE( nout, fmt = 9992 )mplusn, prtype,
     $                        weight, m, j, result( j )
                        ELSE
                           WRITE( nout, fmt = 9991 )mplusn, prtype,
     $                        weight, m, j, result( j )
                        END IF
                     END IF
   20             CONTINUE
*
   30          CONTINUE
   40       CONTINUE
   50    CONTINUE
   60 CONTINUE
*
      GO TO 150
*
   70 CONTINUE
*
*     Read in data from file to check accuracy of condition estimation
*     Read input data until N=0
*
      nptknt = 0
*
   80 CONTINUE
      READ( nin, fmt = *, END = 140 )mplusn
      IF( mplusn.EQ.0 )
     $   GO TO 140
      READ( nin, fmt = *, END = 140 )n
      DO 90 i = 1, mplusn
         READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
   90 CONTINUE
      DO 100 i = 1, mplusn
         READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
  100 CONTINUE
      READ( nin, fmt = * )pltru, diftru
*
      nptknt = nptknt + 1
      fs = .true.
      k = 0
      m = mplusn - n
*
      CALL clacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
      CALL clacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
*
*     Compute the Schur factorization while swapping the
*     m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
*
      CALL cggesx( 'V', 'V', 'S', clctsx, 'B', mplusn, ai, lda, bi, lda,
     $             mm, alpha, beta, q, lda, z, lda, pl, difest, work,
     $             lwork, rwork, iwork, liwork, bwork, linfo )
*
      IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
         result( 1 ) = ulpinv
         WRITE( nout, fmt = 9998 )'CGGESX', linfo, mplusn, nptknt
         GO TO 130
      END IF
*
*     Compute the norm(A, B)
*        (should this be norm of (A,B) or (AI,BI)?)
*
      CALL clacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
      CALL clacpy( 'Full', mplusn, mplusn, bi, lda,
     $             work( mplusn*mplusn+1 ), mplusn )
      abnrm = clange( 'Fro', mplusn, 2*mplusn, work, mplusn, rwork )
*
*     Do tests (1) to (4)
*
      CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
     $             rwork, result( 1 ) )
      CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
     $             rwork, result( 2 ) )
      CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
     $             rwork, result( 3 ) )
      CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
     $             rwork, result( 4 ) )
*
*     Do tests (5) and (6): check Schur form of A and compare
*     eigenvalues with diagonals.
*
      ntest = 6
      temp1 = zero
      result( 5 ) = zero
      result( 6 ) = zero
*
      DO 110 j = 1, mplusn
         ilabad = .false.
         temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
     $           max( smlnum, abs1( alpha( j ) ), abs1( ai( j, j ) ) )+
     $           abs1( beta( j )-bi( j, j ) ) /
     $           max( smlnum, abs1( beta( j ) ), abs1( bi( j, j ) ) ) )
     $            / ulp
         IF( j.LT.mplusn ) THEN
            IF( ai( j+1, j ).NE.zero ) THEN
               ilabad = .true.
               result( 5 ) = ulpinv
            END IF
         END IF
         IF( j.GT.1 ) THEN
            IF( ai( j, j-1 ).NE.zero ) THEN
               ilabad = .true.
               result( 5 ) = ulpinv
            END IF
         END IF
         temp1 = max( temp1, temp2 )
         IF( ilabad ) THEN
            WRITE( nout, fmt = 9997 )j, mplusn, nptknt
         END IF
  110 CONTINUE
      result( 6 ) = temp1
*
*     Test (7) (if sorting worked)  <--------- need to be checked.
*
      ntest = 7
      result( 7 ) = zero
      IF( linfo.EQ.mplusn+3 )
     $   result( 7 ) = ulpinv
*
*     Test (8): compare the estimated value of DIF and its true value.
*
      ntest = 8
      result( 8 ) = zero
      IF( difest( 2 ).EQ.zero ) THEN
         IF( diftru.GT.abnrm*ulp )
     $      result( 8 ) = ulpinv
      ELSE IF( diftru.EQ.zero ) THEN
         IF( difest( 2 ).GT.abnrm*ulp )
     $      result( 8 ) = ulpinv
      ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
     $         ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
         result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
      END IF
*
*     Test (9)
*
      ntest = 9
      result( 9 ) = zero
      IF( linfo.EQ.( mplusn+2 ) ) THEN
         IF( diftru.GT.abnrm*ulp )
     $      result( 9 ) = ulpinv
         IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
     $      result( 9 ) = ulpinv
         IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
     $      result( 9 ) = ulpinv
      END IF
*
*     Test (10): compare the estimated value of PL and it true value.
*
      ntest = 10
      result( 10 ) = zero
      IF( pl( 1 ).EQ.zero ) THEN
         IF( pltru.GT.abnrm*ulp )
     $      result( 10 ) = ulpinv
      ELSE IF( pltru.EQ.zero ) THEN
         IF( pl( 1 ).GT.abnrm*ulp )
     $      result( 10 ) = ulpinv
      ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
     $         ( pltru*thresh.LT.pl( 1 ) ) ) THEN
         result( 10 ) = ulpinv
      END IF
*
      ntestt = ntestt + ntest
*
*     Print out tests which fail.
*
      DO 120 j = 1, ntest
         IF( result( j ).GE.thresh ) THEN
*
*           If this is the first test to fail,
*           print a header to the data file.
*
            IF( nerrs.EQ.0 ) THEN
               WRITE( nout, fmt = 9996 )'CGX'
*
*              Matrix types
*
               WRITE( nout, fmt = 9995 )
*
*              Tests performed
*
               WRITE( nout, fmt = 9993 )'unitary', '''', 'transpose',
     $            ( '''', i = 1, 4 )
*
            END IF
            nerrs = nerrs + 1
            IF( result( j ).LT.10000.0 ) THEN
               WRITE( nout, fmt = 9990 )nptknt, mplusn, j, result( j )
            ELSE
               WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
            END IF
         END IF
*
  120 CONTINUE
*
  130 CONTINUE
      GO TO 80
  140 CONTINUE
*
  150 CONTINUE
*
*     Summary
*
      CALL alasvm( 'CGX', nout, nerrs, ntestt, 0 )
*
      work( 1 ) = maxwrk
*
      RETURN
*
 9999 FORMAT( ' CDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
     $      i6, ', JTYPE=', i6, ')' )
*
 9998 FORMAT( ' CDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
     $      i6, ', Input Example #', i2, ')' )
*
 9997 FORMAT( ' CDRGSX: S not in Schur form at eigenvalue ', i6, '.',
     $      / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
*
 9996 FORMAT( / 1x, a3, ' -- Complex Expert Generalized Schur form',
     $      ' problem driver' )
*
 9995 FORMAT( 'Input Example' )
*
 9994 FORMAT( ' Matrix types: ', /
     $      '  1:  A is a block diagonal matrix of Jordan blocks ',
     $      'and B is the identity ', / '      matrix, ',
     $      / '  2:  A and B are upper triangular matrices, ',
     $      / '  3:  A and B are as type 2, but each second diagonal ',
     $      'block in A_11 and ', /
     $      '      each third diagonal block in A_22 are 2x2 blocks,',
     $      / '  4:  A and B are block diagonal matrices, ',
     $      / '  5:  (A,B) has potentially close or common ',
     $      'eigenvalues.', / )
*
 9993 FORMAT( / ' Tests performed:  (S is Schur, T is triangular, ',
     $      'Q and Z are ', a, ',', / 19x,
     $      ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
     $      / '  1 = | A - Q S Z', a,
     $      ' | / ( |A| n ulp )      2 = | B - Q T Z', a,
     $      ' | / ( |B| n ulp )', / '  3 = | I - QQ', a,
     $      ' | / ( n ulp )             4 = | I - ZZ', a,
     $      ' | / ( n ulp )', / '  5 = 1/ULP  if A is not in ',
     $      'Schur form S', / '  6 = difference between (alpha,beta)',
     $      ' and diagonals of (S,T)', /
     $      '  7 = 1/ULP  if SDIM is not the correct number of ',
     $      'selected eigenvalues', /
     $      '  8 = 1/ULP  if DIFEST/DIFTRU > 10*THRESH or ',
     $      'DIFTRU/DIFEST > 10*THRESH',
     $      / '  9 = 1/ULP  if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
     $      'when reordering fails', /
     $      ' 10 = 1/ULP  if PLEST/PLTRU > THRESH or ',
     $      'PLTRU/PLEST > THRESH', /
     $      '    ( Test 10 is only for input examples )', / )
 9992 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
     $      ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
 9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
     $      ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, e10.3 )
 9990 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
     $      ' result ', i2, ' is', 0p, f8.2 )
 9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
     $      ' result ', i2, ' is', 1p, e10.3 )
*
*     End of CDRGSX
*

Here is the call graph for this function:

Here is the caller graph for this function: