Move LAPACK trunk into position.

1 files changed, 905 insertions, 0 deletions

diff --git a/TESTING/EIG/sdrgsx.f b/TESTING/EIG/sdrgsx.f
new file mode 100644
index 00000000..ee79a68e
--- /dev/null
+++ b/TESTING/EIG/sdrgsx.f
@@ -0,0 +1,905 @@
+      SUBROUTINE SDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
+     $                   AI, BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S,
+     $                   WORK, LWORK, IWORK, LIWORK, BWORK, INFO )
+*
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      INTEGER            INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
+     $                   NOUT, NSIZE
+      REAL               THRESH
+*     ..
+*     .. Array Arguments ..
+      LOGICAL            BWORK( * )
+      INTEGER            IWORK( * )
+      REAL               A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
+     $                   ALPHAR( * ), B( LDA, * ), BETA( * ),
+     $                   BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ),
+     $                   WORK( * ), Z( LDA, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
+*  problem expert driver SGGESX.
+*
+*  SGGESX factors A and B as Q S Z' and Q T Z', where ' means
+*  transpose, T is upper triangular, S is in generalized Schur form
+*  (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
+*  the 2x2 blocks corresponding to complex conjugate pairs of
+*  generalized eigenvalues), and Q and Z are orthogonal.  It also
+*  computes the generalized eigenvalues (alpha(1),beta(1)), ...,
+*  (alpha(n),beta(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 SDRGSX is called with NSIZE > 0, five (5) types of built-in
+*  matrix pairs are used to test the routine SGGESX.
+*
+*  When SDRGSX is called with NSIZE = 0, it reads in test matrix data
+*  to test SGGESX.
+*
+*  For each matrix pair, the following tests will be performed and
+*  compared with the threshhold 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. quasi-triangular form)
+*
+*  (6)   maximum over j of D(j)  where:
+*
+*        if alpha(j) is real:
+*                      |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
+*            D(j) = ------------------------ + -----------------------
+*                   max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
+*
+*        if alpha(j) is complex:
+*                                  | det( s S - w T ) |
+*            D(j) = ---------------------------------------------------
+*                   ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
+*
+*            and S and T are here the 2 x 2 diagonal blocks of S and T
+*            corresponding to the j-th and j+1-th eigenvalues.
+*
+*  (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 SGGESX, the reordering "failed"
+*        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 above 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 SGESVD) 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).
+*
+*
+*  Arguments
+*  =========
+*
+*  NSIZE   (input) 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.
+*
+*  NCMAX   (input) INTEGER
+*          Maximum allowable NMAX for generating Kroneker matrix
+*          in call to SLAKF2
+*
+*  THRESH  (input) REAL
+*          A test will count as "failed" if the "error", 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.
+*
+*  NIN     (input) INTEGER
+*          The FORTRAN unit number for reading in the data file of
+*          problems to solve.
+*
+*  NOUT    (input) INTEGER
+*          The FORTRAN unit number for printing out error messages
+*          (e.g., if a routine returns IINFO not equal to 0.)
+*
+*  A       (workspace) REAL array, dimension (LDA, NSIZE)
+*          Used to store the matrix whose eigenvalues are to be
+*          computed.  On exit, A contains the last matrix actually used.
+*
+*  LDA     (input) 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.
+*
+*  B       (workspace) REAL array, dimension (LDA, NSIZE)
+*          Used to store the matrix whose eigenvalues are to be
+*          computed.  On exit, B contains the last matrix actually used.
+*
+*  AI      (workspace) REAL array, dimension (LDA, NSIZE)
+*          Copy of A, modified by SGGESX.
+*
+*  BI      (workspace) REAL array, dimension (LDA, NSIZE)
+*          Copy of B, modified by SGGESX.
+*
+*  Z       (workspace) REAL array, dimension (LDA, NSIZE)
+*          Z holds the left Schur vectors computed by SGGESX.
+*
+*  Q       (workspace) REAL array, dimension (LDA, NSIZE)
+*          Q holds the right Schur vectors computed by SGGESX.
+*
+*  ALPHAR  (workspace) REAL array, dimension (NSIZE)
+*  ALPHAI  (workspace) REAL array, dimension (NSIZE)
+*  BETA    (workspace) REAL array, dimension (NSIZE)
+*          On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
+*
+*  C       (workspace) REAL array, dimension (LDC, LDC)
+*          Store the matrix generated by subroutine SLAKF2, this is the
+*          matrix formed by Kronecker products used for estimating
+*          DIF.
+*
+*  LDC     (input) INTEGER
+*          The leading dimension of C. LDC >= max(1, LDA*LDA/2 ).
+*
+*  S       (workspace) REAL array, dimension (LDC)
+*          Singular values of C
+*
+*  WORK    (workspace) REAL array, dimension (LWORK)
+*
+*  LWORK   (input) INTEGER
+*          The dimension of the array WORK.
+*          LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) )
+*
+*  IWORK   (workspace) INTEGER array, dimension (LIWORK)
+*
+*  LIWORK  (input) INTEGER
+*          The dimension of the array IWORK. LIWORK >= NSIZE + 6.
+*
+*  BWORK   (workspace) LOGICAL array, dimension (LDA)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value.
+*          > 0:  A routine returned an error code.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO, ONE, TEN
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            ILABAD
+      CHARACTER          SENSE
+      INTEGER            BDSPAC, I, I1, IFUNC, IINFO, 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
+*     ..
+*     .. Local Arrays ..
+      REAL               DIFEST( 2 ), PL( 2 ), RESULT( 10 )
+*     ..
+*     .. External Functions ..
+      LOGICAL            SLCTSX
+      INTEGER            ILAENV
+      REAL               SLAMCH, SLANGE
+      EXTERNAL           SLCTSX, ILAENV, SLAMCH, SLANGE
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALASVM, SGESVD, SGET51, SGET53, SGGESX, SLABAD,
+     $                   SLACPY, SLAKF2, SLASET, SLATM5, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, SQRT
+*     ..
+*     .. Scalars in Common ..
+      LOGICAL            FS
+      INTEGER            K, M, MPLUSN, N
+*     ..
+*     .. Common blocks ..
+      COMMON             / MN / M, N, MPLUSN, K, FS
+*     ..
+*     .. 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 = -17
+      ELSE IF( LIWORK.LT.NSIZE+6 ) 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
+c        MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2-2 )
+         MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 )
+*
+*        workspace for sggesx
+*
+         MAXWRK = 9*( NSIZE+1 ) + NSIZE*
+     $            ILAENV( 1, 'SGEQRF', ' ', NSIZE, 1, NSIZE, 0 )
+         MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE*
+     $            ILAENV( 1, 'SORGQR', ' ', NSIZE, 1, NSIZE, -1 ) )
+*
+*        workspace for sgesvd
+*
+         BDSPAC = 5*NSIZE*NSIZE / 2
+         MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE*
+     $            ILAENV( 1, 'SGEBRD', ' ', 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 = -19
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'SDRGSX', -INFO )
+         RETURN
+      END IF
+*
+*     Important constants
+*
+      ULP = SLAMCH( 'P' )
+      ULPINV = ONE / ULP
+      SMLNUM = SLAMCH( 'S' ) / ULP
+      BIGNUM = ONE / SMLNUM
+      CALL SLABAD( SMLNUM, BIGNUM )
+      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 SGGESX, 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 SLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI,
+     $                         LDA )
+                  CALL SLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI,
+     $                         LDA )
+*
+                  CALL SLATM5( 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 SLCTSX
+*                 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 SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
+                  CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
+*
+                  CALL SGGESX( 'V', 'V', 'S', SLCTSX, SENSE, MPLUSN, AI,
+     $                         LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA,
+     $                         Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK,
+     $                         IWORK, LIWORK, BWORK, LINFO )
+*
+                  IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
+                     RESULT( 1 ) = ULPINV
+                     WRITE( NOUT, FMT = 9999 )'SGGESX', LINFO, MPLUSN,
+     $                  PRTYPE
+                     INFO = LINFO
+                     GO TO 30
+                  END IF
+*
+*                 Compute the norm(A, B)
+*
+                  CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
+     $                         MPLUSN )
+                  CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
+     $                         WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
+                  ABNRM = SLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
+     $                    WORK )
+*
+*                 Do tests (1) to (4)
+*
+                  CALL SGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
+     $                         LDA, WORK, RESULT( 1 ) )
+                  CALL SGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
+     $                         LDA, WORK, RESULT( 2 ) )
+                  CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
+     $                         LDA, WORK, RESULT( 3 ) )
+                  CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
+     $                         LDA, WORK, 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.
+                     IF( ALPHAI( J ).EQ.ZERO ) THEN
+                        TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
+     $                          MAX( SMLNUM, ABS( ALPHAR( J ) ),
+     $                          ABS( AI( J, J ) ) )+
+     $                          ABS( BETA( J )-BI( J, J ) ) /
+     $                          MAX( SMLNUM, ABS( BETA( J ) ),
+     $                          ABS( 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
+                     ELSE
+                        IF( ALPHAI( J ).GT.ZERO ) THEN
+                           I1 = J
+                        ELSE
+                           I1 = J - 1
+                        END IF
+                        IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
+                           ILABAD = .TRUE.
+                        ELSE IF( I1.LT.MPLUSN-1 ) THEN
+                           IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
+                              ILABAD = .TRUE.
+                              RESULT( 5 ) = ULPINV
+                           END IF
+                        ELSE IF( I1.GT.1 ) THEN
+                           IF( AI( I1, I1-1 ).NE.ZERO ) THEN
+                              ILABAD = .TRUE.
+                              RESULT( 5 ) = ULPINV
+                           END IF
+                        END IF
+                        IF( .NOT.ILABAD ) THEN
+                           CALL SGET53( AI( I1, I1 ), LDA, BI( I1, I1 ),
+     $                                  LDA, BETA( J ), ALPHAR( J ),
+     $                                  ALPHAI( J ), TEMP2, IINFO )
+                           IF( IINFO.GE.3 ) THEN
+                              WRITE( NOUT, FMT = 9997 )IINFO, J,
+     $                           MPLUSN, PRTYPE
+                              INFO = ABS( IINFO )
+                           END IF
+                        ELSE
+                           TEMP2 = ULPINV
+                        END IF
+                     END IF
+                     TEMP1 = MAX( TEMP1, TEMP2 )
+                     IF( ILABAD ) THEN
+                        WRITE( NOUT, FMT = 9996 )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 causes, there are
+*                    almost same number of test cases fail the test.
+*
+                     CALL SLAKF2( MM, MPLUSN-MM, AI, LDA,
+     $                            AI( MM+1, MM+1 ), BI,
+     $                            BI( MM+1, MM+1 ), C, LDC )
+*
+                     CALL SGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
+     $                            1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
+     $                            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 = 9995 )'SGX'
+*
+*                          Matrix types
+*
+                           WRITE( NOUT, FMT = 9993 )
+*
+*                          Tests performed
+*
+                           WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
+     $                        'transpose', ( '''', I = 1, 4 )
+*
+                        END IF
+                        NERRS = NERRS + 1
+                        IF( RESULT( J ).LT.10000.0 ) THEN
+                           WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
+     $                        WEIGHT, M, J, RESULT( J )
+                        ELSE
+                           WRITE( NOUT, FMT = 9990 )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 SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
+      CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
+*
+*     Compute the Schur factorization while swaping the
+*     m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
+*
+      CALL SGGESX( 'V', 'V', 'S', SLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
+     $             MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST,
+     $             WORK, LWORK, IWORK, LIWORK, BWORK, LINFO )
+*
+      IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
+         RESULT( 1 ) = ULPINV
+         WRITE( NOUT, FMT = 9998 )'SGGESX', LINFO, MPLUSN, NPTKNT
+         GO TO 130
+      END IF
+*
+*     Compute the norm(A, B)
+*        (should this be norm of (A,B) or (AI,BI)?)
+*
+      CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
+      CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
+     $             WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
+      ABNRM = SLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK )
+*
+*     Do tests (1) to (4)
+*
+      CALL SGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
+     $             RESULT( 1 ) )
+      CALL SGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
+     $             RESULT( 2 ) )
+      CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
+     $             RESULT( 3 ) )
+      CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
+     $             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.
+         IF( ALPHAI( J ).EQ.ZERO ) THEN
+            TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
+     $              MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J,
+     $              J ) ) )+ABS( BETA( J )-BI( J, J ) ) /
+     $              MAX( SMLNUM, ABS( BETA( J ) ), ABS( 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
+         ELSE
+            IF( ALPHAI( J ).GT.ZERO ) THEN
+               I1 = J
+            ELSE
+               I1 = J - 1
+            END IF
+            IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
+               ILABAD = .TRUE.
+            ELSE IF( I1.LT.MPLUSN-1 ) THEN
+               IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
+                  ILABAD = .TRUE.
+                  RESULT( 5 ) = ULPINV
+               END IF
+            ELSE IF( I1.GT.1 ) THEN
+               IF( AI( I1, I1-1 ).NE.ZERO ) THEN
+                  ILABAD = .TRUE.
+                  RESULT( 5 ) = ULPINV
+               END IF
+            END IF
+            IF( .NOT.ILABAD ) THEN
+               CALL SGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA,
+     $                      BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2,
+     $                      IINFO )
+               IF( IINFO.GE.3 ) THEN
+                  WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT
+                  INFO = ABS( IINFO )
+               END IF
+            ELSE
+               TEMP2 = ULPINV
+            END IF
+         END IF
+         TEMP1 = MAX( TEMP1, TEMP2 )
+         IF( ILABAD ) THEN
+            WRITE( NOUT, FMT = 9996 )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 = 9995 )'SGX'
+*
+*              Matrix types
+*
+               WRITE( NOUT, FMT = 9994 )
+*
+*              Tests performed
+*
+               WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
+     $            'transpose', ( '''', I = 1, 4 )
+*
+            END IF
+            NERRS = NERRS + 1
+            IF( RESULT( J ).LT.10000.0 ) THEN
+               WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
+            ELSE
+               WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J )
+            END IF
+         END IF
+*
+  120 CONTINUE
+*
+  130 CONTINUE
+      GO TO 80
+  140 CONTINUE
+*
+  150 CONTINUE
+*
+*     Summary
+*
+      CALL ALASVM( 'SGX', NOUT, NERRS, NTESTT, 0 )
+*
+      WORK( 1 ) = MAXWRK
+*
+      RETURN
+*
+ 9999 FORMAT( ' SDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
+     $      I6, ', JTYPE=', I6, ')' )
+*
+ 9998 FORMAT( ' SDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
+     $      I6, ', Input Example #', I2, ')' )
+*
+ 9997 FORMAT( ' SDRGSX: SGET53 returned INFO=', I1, ' for eigenvalue ',
+     $      I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
+*
+ 9996 FORMAT( ' SDRGSX: S not in Schur form at eigenvalue ', I6, '.',
+     $      / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
+*
+ 9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form',
+     $      ' problem driver' )
+*
+ 9994 FORMAT( 'Input Example' )
+*
+ 9993 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 diaongal block in A_22 are 2x2 blocks,',
+     $      / '  4:  A and B are block diagonal matrices, ',
+     $      / '  5:  (A,B) has potentially close or common ',
+     $      'eigenvalues.', / )
+*
+ 9992 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 )', / )
+ 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.4,
+     $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
+ 9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.4,
+     $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, E10.4 )
+ 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
+     $      ' result ', I2, ' is', 0P, F8.2 )
+ 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
+     $      ' result ', I2, ' is', 1P, E10.3 )
+*
+*     End of SDRGSX
+*
+      END