LAPACK: TESTING/EIG/sdrgsx.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
00002      $                   AI, BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S,
00003      $                   WORK, LWORK, IWORK, LIWORK, BWORK, INFO )
00004 *
00005 *  -- LAPACK test routine (version 3.1) --
00006 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
00011      $                   NOUT, NSIZE
00012       REAL               THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            BWORK( * )
00016       INTEGER            IWORK( * )
00017       REAL               A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
00018      $                   ALPHAR( * ), B( LDA, * ), BETA( * ),
00019      $                   BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ),
00020      $                   WORK( * ), Z( LDA, * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  SDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
00027 *  problem expert driver SGGESX.
00028 *
00029 *  SGGESX factors A and B as Q S Z' and Q T Z', where ' means
00030 *  transpose, T is upper triangular, S is in generalized Schur form
00031 *  (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
00032 *  the 2x2 blocks corresponding to complex conjugate pairs of
00033 *  generalized eigenvalues), and Q and Z are orthogonal.  It also
00034 *  computes the generalized eigenvalues (alpha(1),beta(1)), ...,
00035 *  (alpha(n),beta(n)). Thus, w(j) = alpha(j)/beta(j) is a root of the
00036 *  characteristic equation
00037 *
00038 *      det( A - w(j) B ) = 0
00039 *
00040 *  Optionally it also reorders the eigenvalues so that a selected
00041 *  cluster of eigenvalues appears in the leading diagonal block of the
00042 *  Schur forms; computes a reciprocal condition number for the average
00043 *  of the selected eigenvalues; and computes a reciprocal condition
00044 *  number for the right and left deflating subspaces corresponding to
00045 *  the selected eigenvalues.
00046 *
00047 *  When SDRGSX is called with NSIZE > 0, five (5) types of built-in
00048 *  matrix pairs are used to test the routine SGGESX.
00049 *
00050 *  When SDRGSX is called with NSIZE = 0, it reads in test matrix data
00051 *  to test SGGESX.
00052 *
00053 *  For each matrix pair, the following tests will be performed and
00054 *  compared with the threshhold THRESH except for the tests (7) and (9):
00055 *
00056 *  (1)   | A - Q S Z' | / ( |A| n ulp )
00057 *
00058 *  (2)   | B - Q T Z' | / ( |B| n ulp )
00059 *
00060 *  (3)   | I - QQ' | / ( n ulp )
00061 *
00062 *  (4)   | I - ZZ' | / ( n ulp )
00063 *
00064 *  (5)   if A is in Schur form (i.e. quasi-triangular form)
00065 *
00066 *  (6)   maximum over j of D(j)  where:
00067 *
00068 *        if alpha(j) is real:
00069 *                      |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
00070 *            D(j) = ------------------------ + -----------------------
00071 *                   max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
00072 *
00073 *        if alpha(j) is complex:
00074 *                                  | det( s S - w T ) |
00075 *            D(j) = ---------------------------------------------------
00076 *                   ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
00077 *
00078 *            and S and T are here the 2 x 2 diagonal blocks of S and T
00079 *            corresponding to the j-th and j+1-th eigenvalues.
00080 *
00081 *  (7)   if sorting worked and SDIM is the number of eigenvalues
00082 *        which were selected.
00083 *
00084 *  (8)   the estimated value DIF does not differ from the true values of
00085 *        Difu and Difl more than a factor 10*THRESH. If the estimate DIF
00086 *        equals zero the corresponding true values of Difu and Difl
00087 *        should be less than EPS*norm(A, B). If the true value of Difu
00088 *        and Difl equal zero, the estimate DIF should be less than
00089 *        EPS*norm(A, B).
00090 *
00091 *  (9)   If INFO = N+3 is returned by SGGESX, the reordering "failed"
00092 *        and we check that DIF = PL = PR = 0 and that the true value of
00093 *        Difu and Difl is < EPS*norm(A, B). We count the events when
00094 *        INFO=N+3.
00095 *
00096 *  For read-in test matrices, the above tests are run except that the
00097 *  exact value for DIF (and PL) is input data.  Additionally, there is
00098 *  one more test run for read-in test matrices:
00099 *
00100 *  (10)  the estimated value PL does not differ from the true value of
00101 *        PLTRU more than a factor THRESH. If the estimate PL equals
00102 *        zero the corresponding true value of PLTRU should be less than
00103 *        EPS*norm(A, B). If the true value of PLTRU equal zero, the
00104 *        estimate PL should be less than EPS*norm(A, B).
00105 *
00106 *  Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1)
00107 *  matrix pairs are generated and tested. NSIZE should be kept small.
00108 *
00109 *  SVD (routine SGESVD) is used for computing the true value of DIF_u
00110 *  and DIF_l when testing the built-in test problems.
00111 *
00112 *  Built-in Test Matrices
00113 *  ======================
00114 *
00115 *  All built-in test matrices are the 2 by 2 block of triangular
00116 *  matrices
00117 *
00118 *           A = [ A11 A12 ]    and      B = [ B11 B12 ]
00119 *               [     A22 ]                 [     B22 ]
00120 *
00121 *  where for different type of A11 and A22 are given as the following.
00122 *  A12 and B12 are chosen so that the generalized Sylvester equation
00123 *
00124 *           A11*R - L*A22 = -A12
00125 *           B11*R - L*B22 = -B12
00126 *
00127 *  have prescribed solution R and L.
00128 *
00129 *  Type 1:  A11 = J_m(1,-1) and A_22 = J_k(1-a,1).
00130 *           B11 = I_m, B22 = I_k
00131 *           where J_k(a,b) is the k-by-k Jordan block with ``a'' on
00132 *           diagonal and ``b'' on superdiagonal.
00133 *
00134 *  Type 2:  A11 = (a_ij) = ( 2(.5-sin(i)) ) and
00135 *           B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m
00136 *           A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and
00137 *           B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k
00138 *
00139 *  Type 3:  A11, A22 and B11, B22 are chosen as for Type 2, but each
00140 *           second diagonal block in A_11 and each third diagonal block
00141 *           in A_22 are made as 2 by 2 blocks.
00142 *
00143 *  Type 4:  A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) )
00144 *              for i=1,...,m,  j=1,...,m and
00145 *           A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) )
00146 *              for i=m+1,...,k,  j=m+1,...,k
00147 *
00148 *  Type 5:  (A,B) and have potentially close or common eigenvalues and
00149 *           very large departure from block diagonality A_11 is chosen
00150 *           as the m x m leading submatrix of A_1:
00151 *                   |  1  b                            |
00152 *                   | -b  1                            |
00153 *                   |        1+d  b                    |
00154 *                   |         -b 1+d                   |
00155 *            A_1 =  |                  d  1            |
00156 *                   |                 -1  d            |
00157 *                   |                        -d  1     |
00158 *                   |                        -1 -d     |
00159 *                   |                               1  |
00160 *           and A_22 is chosen as the k x k leading submatrix of A_2:
00161 *                   | -1  b                            |
00162 *                   | -b -1                            |
00163 *                   |       1-d  b                     |
00164 *                   |       -b  1-d                    |
00165 *            A_2 =  |                 d 1+b            |
00166 *                   |               -1-b d             |
00167 *                   |                       -d  1+b    |
00168 *                   |                      -1+b  -d    |
00169 *                   |                              1-d |
00170 *           and matrix B are chosen as identity matrices (see SLATM5).
00171 *
00172 *
00173 *  Arguments
00174 *  =========
00175 *
00176 *  NSIZE   (input) INTEGER
00177 *          The maximum size of the matrices to use. NSIZE >= 0.
00178 *          If NSIZE = 0, no built-in tests matrices are used, but
00179 *          read-in test matrices are used to test SGGESX.
00180 *
00181 *  NCMAX   (input) INTEGER
00182 *          Maximum allowable NMAX for generating Kroneker matrix
00183 *          in call to SLAKF2
00184 *
00185 *  THRESH  (input) REAL
00186 *          A test will count as "failed" if the "error", computed as
00187 *          described above, exceeds THRESH.  Note that the error
00188 *          is scaled to be O(1), so THRESH should be a reasonably
00189 *          small multiple of 1, e.g., 10 or 100.  In particular,
00190 *          it should not depend on the precision (single vs. double)
00191 *          or the size of the matrix.  THRESH >= 0.
00192 *
00193 *  NIN     (input) INTEGER
00194 *          The FORTRAN unit number for reading in the data file of
00195 *          problems to solve.
00196 *
00197 *  NOUT    (input) INTEGER
00198 *          The FORTRAN unit number for printing out error messages
00199 *          (e.g., if a routine returns IINFO not equal to 0.)
00200 *
00201 *  A       (workspace) REAL array, dimension (LDA, NSIZE)
00202 *          Used to store the matrix whose eigenvalues are to be
00203 *          computed.  On exit, A contains the last matrix actually used.
00204 *
00205 *  LDA     (input) INTEGER
00206 *          The leading dimension of A, B, AI, BI, Z and Q,
00207 *          LDA >= max( 1, NSIZE ). For the read-in test,
00208 *          LDA >= max( 1, N ), N is the size of the test matrices.
00209 *
00210 *  B       (workspace) REAL array, dimension (LDA, NSIZE)
00211 *          Used to store the matrix whose eigenvalues are to be
00212 *          computed.  On exit, B contains the last matrix actually used.
00213 *
00214 *  AI      (workspace) REAL array, dimension (LDA, NSIZE)
00215 *          Copy of A, modified by SGGESX.
00216 *
00217 *  BI      (workspace) REAL array, dimension (LDA, NSIZE)
00218 *          Copy of B, modified by SGGESX.
00219 *
00220 *  Z       (workspace) REAL array, dimension (LDA, NSIZE)
00221 *          Z holds the left Schur vectors computed by SGGESX.
00222 *
00223 *  Q       (workspace) REAL array, dimension (LDA, NSIZE)
00224 *          Q holds the right Schur vectors computed by SGGESX.
00225 *
00226 *  ALPHAR  (workspace) REAL array, dimension (NSIZE)
00227 *  ALPHAI  (workspace) REAL array, dimension (NSIZE)
00228 *  BETA    (workspace) REAL array, dimension (NSIZE)
00229 *          On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
00230 *
00231 *  C       (workspace) REAL array, dimension (LDC, LDC)
00232 *          Store the matrix generated by subroutine SLAKF2, this is the
00233 *          matrix formed by Kronecker products used for estimating
00234 *          DIF.
00235 *
00236 *  LDC     (input) INTEGER
00237 *          The leading dimension of C. LDC >= max(1, LDA*LDA/2 ).
00238 *
00239 *  S       (workspace) REAL array, dimension (LDC)
00240 *          Singular values of C
00241 *
00242 *  WORK    (workspace) REAL array, dimension (LWORK)
00243 *
00244 *  LWORK   (input) INTEGER
00245 *          The dimension of the array WORK.
00246 *          LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) )
00247 *
00248 *  IWORK   (workspace) INTEGER array, dimension (LIWORK)
00249 *
00250 *  LIWORK  (input) INTEGER
00251 *          The dimension of the array IWORK. LIWORK >= NSIZE + 6.
00252 *
00253 *  BWORK   (workspace) LOGICAL array, dimension (LDA)
00254 *
00255 *  INFO    (output) INTEGER
00256 *          = 0:  successful exit
00257 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00258 *          > 0:  A routine returned an error code.
00259 *
00260 *  =====================================================================
00261 *
00262 *     .. Parameters ..
00263       REAL               ZERO, ONE, TEN
00264       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1 )
00265 *     ..
00266 *     .. Local Scalars ..
00267       LOGICAL            ILABAD
00268       CHARACTER          SENSE
00269       INTEGER            BDSPAC, I, I1, IFUNC, IINFO, J, LINFO, MAXWRK,
00270      $                   MINWRK, MM, MN2, NERRS, NPTKNT, NTEST, NTESTT,
00271      $                   PRTYPE, QBA, QBB
00272       REAL               ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
00273      $                   TEMP2, THRSH2, ULP, ULPINV, WEIGHT
00274 *     ..
00275 *     .. Local Arrays ..
00276       REAL               DIFEST( 2 ), PL( 2 ), RESULT( 10 )
00277 *     ..
00278 *     .. External Functions ..
00279       LOGICAL            SLCTSX
00280       INTEGER            ILAENV
00281       REAL               SLAMCH, SLANGE
00282       EXTERNAL           SLCTSX, ILAENV, SLAMCH, SLANGE
00283 *     ..
00284 *     .. External Subroutines ..
00285       EXTERNAL           ALASVM, SGESVD, SGET51, SGET53, SGGESX, SLABAD,
00286      $                   SLACPY, SLAKF2, SLASET, SLATM5, XERBLA
00287 *     ..
00288 *     .. Intrinsic Functions ..
00289       INTRINSIC          ABS, MAX, SQRT
00290 *     ..
00291 *     .. Scalars in Common ..
00292       LOGICAL            FS
00293       INTEGER            K, M, MPLUSN, N
00294 *     ..
00295 *     .. Common blocks ..
00296       COMMON             / MN / M, N, MPLUSN, K, FS
00297 *     ..
00298 *     .. Executable Statements ..
00299 *
00300 *     Check for errors
00301 *
00302       IF( NSIZE.LT.0 ) THEN
00303          INFO = -1
00304       ELSE IF( THRESH.LT.ZERO ) THEN
00305          INFO = -2
00306       ELSE IF( NIN.LE.0 ) THEN
00307          INFO = -3
00308       ELSE IF( NOUT.LE.0 ) THEN
00309          INFO = -4
00310       ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN
00311          INFO = -6
00312       ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN
00313          INFO = -17
00314       ELSE IF( LIWORK.LT.NSIZE+6 ) THEN
00315          INFO = -21
00316       END IF
00317 *
00318 *     Compute workspace
00319 *      (Note: Comments in the code beginning "Workspace:" describe the
00320 *       minimal amount of workspace needed at that point in the code,
00321 *       as well as the preferred amount for good performance.
00322 *       NB refers to the optimal block size for the immediately
00323 *       following subroutine, as returned by ILAENV.)
00324 *
00325       MINWRK = 1
00326       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
00327 c        MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2-2 )
00328          MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 )
00329 *
00330 *        workspace for sggesx
00331 *
00332          MAXWRK = 9*( NSIZE+1 ) + NSIZE*
00333      $            ILAENV( 1, 'SGEQRF', ' ', NSIZE, 1, NSIZE, 0 )
00334          MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE*
00335      $            ILAENV( 1, 'SORGQR', ' ', NSIZE, 1, NSIZE, -1 ) )
00336 *
00337 *        workspace for sgesvd
00338 *
00339          BDSPAC = 5*NSIZE*NSIZE / 2
00340          MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE*
00341      $            ILAENV( 1, 'SGEBRD', ' ', NSIZE*NSIZE / 2,
00342      $            NSIZE*NSIZE / 2, -1, -1 ) )
00343          MAXWRK = MAX( MAXWRK, BDSPAC )
00344 *
00345          MAXWRK = MAX( MAXWRK, MINWRK )
00346 *
00347          WORK( 1 ) = MAXWRK
00348       END IF
00349 *
00350       IF( LWORK.LT.MINWRK )
00351      $   INFO = -19
00352 *
00353       IF( INFO.NE.0 ) THEN
00354          CALL XERBLA( 'SDRGSX', -INFO )
00355          RETURN
00356       END IF
00357 *
00358 *     Important constants
00359 *
00360       ULP = SLAMCH( 'P' )
00361       ULPINV = ONE / ULP
00362       SMLNUM = SLAMCH( 'S' ) / ULP
00363       BIGNUM = ONE / SMLNUM
00364       CALL SLABAD( SMLNUM, BIGNUM )
00365       THRSH2 = TEN*THRESH
00366       NTESTT = 0
00367       NERRS = 0
00368 *
00369 *     Go to the tests for read-in matrix pairs
00370 *
00371       IFUNC = 0
00372       IF( NSIZE.EQ.0 )
00373      $   GO TO 70
00374 *
00375 *     Test the built-in matrix pairs.
00376 *     Loop over different functions (IFUNC) of SGGESX, types (PRTYPE)
00377 *     of test matrices, different size (M+N)
00378 *
00379       PRTYPE = 0
00380       QBA = 3
00381       QBB = 4
00382       WEIGHT = SQRT( ULP )
00383 *
00384       DO 60 IFUNC = 0, 3
00385          DO 50 PRTYPE = 1, 5
00386             DO 40 M = 1, NSIZE - 1
00387                DO 30 N = 1, NSIZE - M
00388 *
00389                   WEIGHT = ONE / WEIGHT
00390                   MPLUSN = M + N
00391 *
00392 *                 Generate test matrices
00393 *
00394                   FS = .TRUE.
00395                   K = 0
00396 *
00397                   CALL SLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI,
00398      $                         LDA )
00399                   CALL SLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI,
00400      $                         LDA )
00401 *
00402                   CALL SLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ),
00403      $                         LDA, AI( 1, M+1 ), LDA, BI, LDA,
00404      $                         BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA,
00405      $                         Q, LDA, Z, LDA, WEIGHT, QBA, QBB )
00406 *
00407 *                 Compute the Schur factorization and swapping the
00408 *                 m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
00409 *                 Swapping is accomplished via the function SLCTSX
00410 *                 which is supplied below.
00411 *
00412                   IF( IFUNC.EQ.0 ) THEN
00413                      SENSE = 'N'
00414                   ELSE IF( IFUNC.EQ.1 ) THEN
00415                      SENSE = 'E'
00416                   ELSE IF( IFUNC.EQ.2 ) THEN
00417                      SENSE = 'V'
00418                   ELSE IF( IFUNC.EQ.3 ) THEN
00419                      SENSE = 'B'
00420                   END IF
00421 *
00422                   CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
00423                   CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
00424 *
00425                   CALL SGGESX( 'V', 'V', 'S', SLCTSX, SENSE, MPLUSN, AI,
00426      $                         LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA,
00427      $                         Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK,
00428      $                         IWORK, LIWORK, BWORK, LINFO )
00429 *
00430                   IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
00431                      RESULT( 1 ) = ULPINV
00432                      WRITE( NOUT, FMT = 9999 )'SGGESX', LINFO, MPLUSN,
00433      $                  PRTYPE
00434                      INFO = LINFO
00435                      GO TO 30
00436                   END IF
00437 *
00438 *                 Compute the norm(A, B)
00439 *
00440                   CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
00441      $                         MPLUSN )
00442                   CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
00443      $                         WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
00444                   ABNRM = SLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
00445      $                    WORK )
00446 *
00447 *                 Do tests (1) to (4)
00448 *
00449                   CALL SGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
00450      $                         LDA, WORK, RESULT( 1 ) )
00451                   CALL SGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
00452      $                         LDA, WORK, RESULT( 2 ) )
00453                   CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
00454      $                         LDA, WORK, RESULT( 3 ) )
00455                   CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
00456      $                         LDA, WORK, RESULT( 4 ) )
00457                   NTEST = 4
00458 *
00459 *                 Do tests (5) and (6): check Schur form of A and
00460 *                 compare eigenvalues with diagonals.
00461 *
00462                   TEMP1 = ZERO
00463                   RESULT( 5 ) = ZERO
00464                   RESULT( 6 ) = ZERO
00465 *
00466                   DO 10 J = 1, MPLUSN
00467                      ILABAD = .FALSE.
00468                      IF( ALPHAI( J ).EQ.ZERO ) THEN
00469                         TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
00470      $                          MAX( SMLNUM, ABS( ALPHAR( J ) ),
00471      $                          ABS( AI( J, J ) ) )+
00472      $                          ABS( BETA( J )-BI( J, J ) ) /
00473      $                          MAX( SMLNUM, ABS( BETA( J ) ),
00474      $                          ABS( BI( J, J ) ) ) ) / ULP
00475                         IF( J.LT.MPLUSN ) THEN
00476                            IF( AI( J+1, J ).NE.ZERO ) THEN
00477                               ILABAD = .TRUE.
00478                               RESULT( 5 ) = ULPINV
00479                            END IF
00480                         END IF
00481                         IF( J.GT.1 ) THEN
00482                            IF( AI( J, J-1 ).NE.ZERO ) THEN
00483                               ILABAD = .TRUE.
00484                               RESULT( 5 ) = ULPINV
00485                            END IF
00486                         END IF
00487                      ELSE
00488                         IF( ALPHAI( J ).GT.ZERO ) THEN
00489                            I1 = J
00490                         ELSE
00491                            I1 = J - 1
00492                         END IF
00493                         IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
00494                            ILABAD = .TRUE.
00495                         ELSE IF( I1.LT.MPLUSN-1 ) THEN
00496                            IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
00497                               ILABAD = .TRUE.
00498                               RESULT( 5 ) = ULPINV
00499                            END IF
00500                         ELSE IF( I1.GT.1 ) THEN
00501                            IF( AI( I1, I1-1 ).NE.ZERO ) THEN
00502                               ILABAD = .TRUE.
00503                               RESULT( 5 ) = ULPINV
00504                            END IF
00505                         END IF
00506                         IF( .NOT.ILABAD ) THEN
00507                            CALL SGET53( AI( I1, I1 ), LDA, BI( I1, I1 ),
00508      $                                  LDA, BETA( J ), ALPHAR( J ),
00509      $                                  ALPHAI( J ), TEMP2, IINFO )
00510                            IF( IINFO.GE.3 ) THEN
00511                               WRITE( NOUT, FMT = 9997 )IINFO, J,
00512      $                           MPLUSN, PRTYPE
00513                               INFO = ABS( IINFO )
00514                            END IF
00515                         ELSE
00516                            TEMP2 = ULPINV
00517                         END IF
00518                      END IF
00519                      TEMP1 = MAX( TEMP1, TEMP2 )
00520                      IF( ILABAD ) THEN
00521                         WRITE( NOUT, FMT = 9996 )J, MPLUSN, PRTYPE
00522                      END IF
00523    10             CONTINUE
00524                   RESULT( 6 ) = TEMP1
00525                   NTEST = NTEST + 2
00526 *
00527 *                 Test (7) (if sorting worked)
00528 *
00529                   RESULT( 7 ) = ZERO
00530                   IF( LINFO.EQ.MPLUSN+3 ) THEN
00531                      RESULT( 7 ) = ULPINV
00532                   ELSE IF( MM.NE.N ) THEN
00533                      RESULT( 7 ) = ULPINV
00534                   END IF
00535                   NTEST = NTEST + 1
00536 *
00537 *                 Test (8): compare the estimated value DIF and its
00538 *                 value. first, compute the exact DIF.
00539 *
00540                   RESULT( 8 ) = ZERO
00541                   MN2 = MM*( MPLUSN-MM )*2
00542                   IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
00543 *
00544 *                    Note: for either following two causes, there are
00545 *                    almost same number of test cases fail the test.
00546 *
00547                      CALL SLAKF2( MM, MPLUSN-MM, AI, LDA,
00548      $                            AI( MM+1, MM+1 ), BI,
00549      $                            BI( MM+1, MM+1 ), C, LDC )
00550 *
00551                      CALL SGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
00552      $                            1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
00553      $                            INFO )
00554                      DIFTRU = S( MN2 )
00555 *
00556                      IF( DIFEST( 2 ).EQ.ZERO ) THEN
00557                         IF( DIFTRU.GT.ABNRM*ULP )
00558      $                     RESULT( 8 ) = ULPINV
00559                      ELSE IF( DIFTRU.EQ.ZERO ) THEN
00560                         IF( DIFEST( 2 ).GT.ABNRM*ULP )
00561      $                     RESULT( 8 ) = ULPINV
00562                      ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
00563      $                        ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
00564                         RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ),
00565      $                                DIFEST( 2 ) / DIFTRU )
00566                      END IF
00567                      NTEST = NTEST + 1
00568                   END IF
00569 *
00570 *                 Test (9)
00571 *
00572                   RESULT( 9 ) = ZERO
00573                   IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
00574                      IF( DIFTRU.GT.ABNRM*ULP )
00575      $                  RESULT( 9 ) = ULPINV
00576                      IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
00577      $                  RESULT( 9 ) = ULPINV
00578                      IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
00579      $                  RESULT( 9 ) = ULPINV
00580                      NTEST = NTEST + 1
00581                   END IF
00582 *
00583                   NTESTT = NTESTT + NTEST
00584 *
00585 *                 Print out tests which fail.
00586 *
00587                   DO 20 J = 1, 9
00588                      IF( RESULT( J ).GE.THRESH ) THEN
00589 *
00590 *                       If this is the first test to fail,
00591 *                       print a header to the data file.
00592 *
00593                         IF( NERRS.EQ.0 ) THEN
00594                            WRITE( NOUT, FMT = 9995 )'SGX'
00595 *
00596 *                          Matrix types
00597 *
00598                            WRITE( NOUT, FMT = 9993 )
00599 *
00600 *                          Tests performed
00601 *
00602                            WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
00603      $                        'transpose', ( '''', I = 1, 4 )
00604 *
00605                         END IF
00606                         NERRS = NERRS + 1
00607                         IF( RESULT( J ).LT.10000.0 ) THEN
00608                            WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
00609      $                        WEIGHT, M, J, RESULT( J )
00610                         ELSE
00611                            WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE,
00612      $                        WEIGHT, M, J, RESULT( J )
00613                         END IF
00614                      END IF
00615    20             CONTINUE
00616 *
00617    30          CONTINUE
00618    40       CONTINUE
00619    50    CONTINUE
00620    60 CONTINUE
00621 *
00622       GO TO 150
00623 *
00624    70 CONTINUE
00625 *
00626 *     Read in data from file to check accuracy of condition estimation
00627 *     Read input data until N=0
00628 *
00629       NPTKNT = 0
00630 *
00631    80 CONTINUE
00632       READ( NIN, FMT = *, END = 140 )MPLUSN
00633       IF( MPLUSN.EQ.0 )
00634      $   GO TO 140
00635       READ( NIN, FMT = *, END = 140 )N
00636       DO 90 I = 1, MPLUSN
00637          READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
00638    90 CONTINUE
00639       DO 100 I = 1, MPLUSN
00640          READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
00641   100 CONTINUE
00642       READ( NIN, FMT = * )PLTRU, DIFTRU
00643 *
00644       NPTKNT = NPTKNT + 1
00645       FS = .TRUE.
00646       K = 0
00647       M = MPLUSN - N
00648 *
00649       CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
00650       CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
00651 *
00652 *     Compute the Schur factorization while swaping the
00653 *     m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
00654 *
00655       CALL SGGESX( 'V', 'V', 'S', SLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
00656      $             MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST,
00657      $             WORK, LWORK, IWORK, LIWORK, BWORK, LINFO )
00658 *
00659       IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
00660          RESULT( 1 ) = ULPINV
00661          WRITE( NOUT, FMT = 9998 )'SGGESX', LINFO, MPLUSN, NPTKNT
00662          GO TO 130
00663       END IF
00664 *
00665 *     Compute the norm(A, B)
00666 *        (should this be norm of (A,B) or (AI,BI)?)
00667 *
00668       CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
00669       CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
00670      $             WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
00671       ABNRM = SLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK )
00672 *
00673 *     Do tests (1) to (4)
00674 *
00675       CALL SGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
00676      $             RESULT( 1 ) )
00677       CALL SGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
00678      $             RESULT( 2 ) )
00679       CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
00680      $             RESULT( 3 ) )
00681       CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
00682      $             RESULT( 4 ) )
00683 *
00684 *     Do tests (5) and (6): check Schur form of A and compare
00685 *     eigenvalues with diagonals.
00686 *
00687       NTEST = 6
00688       TEMP1 = ZERO
00689       RESULT( 5 ) = ZERO
00690       RESULT( 6 ) = ZERO
00691 *
00692       DO 110 J = 1, MPLUSN
00693          ILABAD = .FALSE.
00694          IF( ALPHAI( J ).EQ.ZERO ) THEN
00695             TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
00696      $              MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J,
00697      $              J ) ) )+ABS( BETA( J )-BI( J, J ) ) /
00698      $              MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) )
00699      $               / ULP
00700             IF( J.LT.MPLUSN ) THEN
00701                IF( AI( J+1, J ).NE.ZERO ) THEN
00702                   ILABAD = .TRUE.
00703                   RESULT( 5 ) = ULPINV
00704                END IF
00705             END IF
00706             IF( J.GT.1 ) THEN
00707                IF( AI( J, J-1 ).NE.ZERO ) THEN
00708                   ILABAD = .TRUE.
00709                   RESULT( 5 ) = ULPINV
00710                END IF
00711             END IF
00712          ELSE
00713             IF( ALPHAI( J ).GT.ZERO ) THEN
00714                I1 = J
00715             ELSE
00716                I1 = J - 1
00717             END IF
00718             IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
00719                ILABAD = .TRUE.
00720             ELSE IF( I1.LT.MPLUSN-1 ) THEN
00721                IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
00722                   ILABAD = .TRUE.
00723                   RESULT( 5 ) = ULPINV
00724                END IF
00725             ELSE IF( I1.GT.1 ) THEN
00726                IF( AI( I1, I1-1 ).NE.ZERO ) THEN
00727                   ILABAD = .TRUE.
00728                   RESULT( 5 ) = ULPINV
00729                END IF
00730             END IF
00731             IF( .NOT.ILABAD ) THEN
00732                CALL SGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA,
00733      $                      BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2,
00734      $                      IINFO )
00735                IF( IINFO.GE.3 ) THEN
00736                   WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT
00737                   INFO = ABS( IINFO )
00738                END IF
00739             ELSE
00740                TEMP2 = ULPINV
00741             END IF
00742          END IF
00743          TEMP1 = MAX( TEMP1, TEMP2 )
00744          IF( ILABAD ) THEN
00745             WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT
00746          END IF
00747   110 CONTINUE
00748       RESULT( 6 ) = TEMP1
00749 *
00750 *     Test (7) (if sorting worked)  <--------- need to be checked.
00751 *
00752       NTEST = 7
00753       RESULT( 7 ) = ZERO
00754       IF( LINFO.EQ.MPLUSN+3 )
00755      $   RESULT( 7 ) = ULPINV
00756 *
00757 *     Test (8): compare the estimated value of DIF and its true value.
00758 *
00759       NTEST = 8
00760       RESULT( 8 ) = ZERO
00761       IF( DIFEST( 2 ).EQ.ZERO ) THEN
00762          IF( DIFTRU.GT.ABNRM*ULP )
00763      $      RESULT( 8 ) = ULPINV
00764       ELSE IF( DIFTRU.EQ.ZERO ) THEN
00765          IF( DIFEST( 2 ).GT.ABNRM*ULP )
00766      $      RESULT( 8 ) = ULPINV
00767       ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
00768      $         ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
00769          RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
00770       END IF
00771 *
00772 *     Test (9)
00773 *
00774       NTEST = 9
00775       RESULT( 9 ) = ZERO
00776       IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
00777          IF( DIFTRU.GT.ABNRM*ULP )
00778      $      RESULT( 9 ) = ULPINV
00779          IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
00780      $      RESULT( 9 ) = ULPINV
00781          IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
00782      $      RESULT( 9 ) = ULPINV
00783       END IF
00784 *
00785 *     Test (10): compare the estimated value of PL and it true value.
00786 *
00787       NTEST = 10
00788       RESULT( 10 ) = ZERO
00789       IF( PL( 1 ).EQ.ZERO ) THEN
00790          IF( PLTRU.GT.ABNRM*ULP )
00791      $      RESULT( 10 ) = ULPINV
00792       ELSE IF( PLTRU.EQ.ZERO ) THEN
00793          IF( PL( 1 ).GT.ABNRM*ULP )
00794      $      RESULT( 10 ) = ULPINV
00795       ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
00796      $         ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
00797          RESULT( 10 ) = ULPINV
00798       END IF
00799 *
00800       NTESTT = NTESTT + NTEST
00801 *
00802 *     Print out tests which fail.
00803 *
00804       DO 120 J = 1, NTEST
00805          IF( RESULT( J ).GE.THRESH ) THEN
00806 *
00807 *           If this is the first test to fail,
00808 *           print a header to the data file.
00809 *
00810             IF( NERRS.EQ.0 ) THEN
00811                WRITE( NOUT, FMT = 9995 )'SGX'
00812 *
00813 *              Matrix types
00814 *
00815                WRITE( NOUT, FMT = 9994 )
00816 *
00817 *              Tests performed
00818 *
00819                WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
00820      $            'transpose', ( '''', I = 1, 4 )
00821 *
00822             END IF
00823             NERRS = NERRS + 1
00824             IF( RESULT( J ).LT.10000.0 ) THEN
00825                WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
00826             ELSE
00827                WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J )
00828             END IF
00829          END IF
00830 *
00831   120 CONTINUE
00832 *
00833   130 CONTINUE
00834       GO TO 80
00835   140 CONTINUE
00836 *
00837   150 CONTINUE
00838 *
00839 *     Summary
00840 *
00841       CALL ALASVM( 'SGX', NOUT, NERRS, NTESTT, 0 )
00842 *
00843       WORK( 1 ) = MAXWRK
00844 *
00845       RETURN
00846 *
00847  9999 FORMAT( ' SDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00848      $      I6, ', JTYPE=', I6, ')' )
00849 *
00850  9998 FORMAT( ' SDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00851      $      I6, ', Input Example #', I2, ')' )
00852 *
00853  9997 FORMAT( ' SDRGSX: SGET53 returned INFO=', I1, ' for eigenvalue ',
00854      $      I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
00855 *
00856  9996 FORMAT( ' SDRGSX: S not in Schur form at eigenvalue ', I6, '.',
00857      $      / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
00858 *
00859  9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form',
00860      $      ' problem driver' )
00861 *
00862  9994 FORMAT( 'Input Example' )
00863 *
00864  9993 FORMAT( ' Matrix types: ', /
00865      $      '  1:  A is a block diagonal matrix of Jordan blocks ',
00866      $      'and B is the identity ', / '      matrix, ',
00867      $      / '  2:  A and B are upper triangular matrices, ',
00868      $      / '  3:  A and B are as type 2, but each second diagonal ',
00869      $      'block in A_11 and ', /
00870      $      '      each third diaongal block in A_22 are 2x2 blocks,',
00871      $      / '  4:  A and B are block diagonal matrices, ',
00872      $      / '  5:  (A,B) has potentially close or common ',
00873      $      'eigenvalues.', / )
00874 *
00875  9992 FORMAT( / ' Tests performed:  (S is Schur, T is triangular, ',
00876      $      'Q and Z are ', A, ',', / 19X,
00877      $      ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
00878      $      / '  1 = | A - Q S Z', A,
00879      $      ' | / ( |A| n ulp )      2 = | B - Q T Z', A,
00880      $      ' | / ( |B| n ulp )', / '  3 = | I - QQ', A,
00881      $      ' | / ( n ulp )             4 = | I - ZZ', A,
00882      $      ' | / ( n ulp )', / '  5 = 1/ULP  if A is not in ',
00883      $      'Schur form S', / '  6 = difference between (alpha,beta)',
00884      $      ' and diagonals of (S,T)', /
00885      $      '  7 = 1/ULP  if SDIM is not the correct number of ',
00886      $      'selected eigenvalues', /
00887      $      '  8 = 1/ULP  if DIFEST/DIFTRU > 10*THRESH or ',
00888      $      'DIFTRU/DIFEST > 10*THRESH',
00889      $      / '  9 = 1/ULP  if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
00890      $      'when reordering fails', /
00891      $      ' 10 = 1/ULP  if PLEST/PLTRU > THRESH or ',
00892      $      'PLTRU/PLEST > THRESH', /
00893      $      '    ( Test 10 is only for input examples )', / )
00894  9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.4,
00895      $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
00896  9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.4,
00897      $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, E10.4 )
00898  9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00899      $      ' result ', I2, ' is', 0P, F8.2 )
00900  9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00901      $      ' result ', I2, ' is', 1P, E10.3 )
00902 *
00903 *     End of SDRGSX
00904 *
00905       END