LAPACK: TESTING/EIG/ddrgsx.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
00002      $                   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       DOUBLE PRECISION   THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            BWORK( * )
00016       INTEGER            IWORK( * )
00017       DOUBLE PRECISION   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 *  DDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
00027 *  problem expert driver DGGESX.
00028 *
00029 *  DGGESX 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 DDRGSX is called with NSIZE > 0, five (5) types of built-in
00048 *  matrix pairs are used to test the routine DGGESX.
00049 *
00050 *  When DDRGSX is called with NSIZE = 0, it reads in test matrix data
00051 *  to test DGGESX.
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 DGGESX, 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 DGESVD) 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 DLATM5).
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 DGGESX.
00180 *
00181 *  NCMAX   (input) INTEGER
00182 *          Maximum allowable NMAX for generating Kroneker matrix
00183 *          in call to DLAKF2
00184 *
00185 *  THRESH  (input) DOUBLE PRECISION
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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDA, NSIZE)
00215 *          Copy of A, modified by DGGESX.
00216 *
00217 *  BI      (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
00218 *          Copy of B, modified by DGGESX.
00219 *
00220 *  Z       (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
00221 *          Z holds the left Schur vectors computed by DGGESX.
00222 *
00223 *  Q       (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
00224 *          Q holds the right Schur vectors computed by DGGESX.
00225 *
00226 *  ALPHAR  (workspace) DOUBLE PRECISION array, dimension (NSIZE)
00227 *  ALPHAI  (workspace) DOUBLE PRECISION array, dimension (NSIZE)
00228 *  BETA    (workspace) DOUBLE PRECISION array, dimension (NSIZE)
00229 *          On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
00230 *
00231 *  C       (workspace) DOUBLE PRECISION array, dimension (LDC, LDC)
00232 *          Store the matrix generated by subroutine DLAKF2, 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) DOUBLE PRECISION array, dimension (LDC)
00240 *          Singular values of C
00241 *
00242 *  WORK    (workspace) DOUBLE PRECISION 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       DOUBLE PRECISION   ZERO, ONE, TEN
00264       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+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       DOUBLE PRECISION   ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
00273      $                   TEMP2, THRSH2, ULP, ULPINV, WEIGHT
00274 *     ..
00275 *     .. Local Arrays ..
00276       DOUBLE PRECISION   DIFEST( 2 ), PL( 2 ), RESULT( 10 )
00277 *     ..
00278 *     .. External Functions ..
00279       LOGICAL            DLCTSX
00280       INTEGER            ILAENV
00281       DOUBLE PRECISION   DLAMCH, DLANGE
00282       EXTERNAL           DLCTSX, ILAENV, DLAMCH, DLANGE
00283 *     ..
00284 *     .. External Subroutines ..
00285       EXTERNAL           ALASVM, DGESVD, DGET51, DGET53, DGGESX, DLABAD,
00286      $                   DLACPY, DLAKF2, DLASET, DLATM5, 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          MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 )
00328 *
00329 *        workspace for sggesx
00330 *
00331          MAXWRK = 9*( NSIZE+1 ) + NSIZE*
00332      $            ILAENV( 1, 'DGEQRF', ' ', NSIZE, 1, NSIZE, 0 )
00333          MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE*
00334      $            ILAENV( 1, 'DORGQR', ' ', NSIZE, 1, NSIZE, -1 ) )
00335 *
00336 *        workspace for dgesvd
00337 *
00338          BDSPAC = 5*NSIZE*NSIZE / 2
00339          MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE*
00340      $            ILAENV( 1, 'DGEBRD', ' ', NSIZE*NSIZE / 2,
00341      $            NSIZE*NSIZE / 2, -1, -1 ) )
00342          MAXWRK = MAX( MAXWRK, BDSPAC )
00343 *
00344          MAXWRK = MAX( MAXWRK, MINWRK )
00345 *
00346          WORK( 1 ) = MAXWRK
00347       END IF
00348 *
00349       IF( LWORK.LT.MINWRK )
00350      $   INFO = -19
00351 *
00352       IF( INFO.NE.0 ) THEN
00353          CALL XERBLA( 'DDRGSX', -INFO )
00354          RETURN
00355       END IF
00356 *
00357 *     Important constants
00358 *
00359       ULP = DLAMCH( 'P' )
00360       ULPINV = ONE / ULP
00361       SMLNUM = DLAMCH( 'S' ) / ULP
00362       BIGNUM = ONE / SMLNUM
00363       CALL DLABAD( SMLNUM, BIGNUM )
00364       THRSH2 = TEN*THRESH
00365       NTESTT = 0
00366       NERRS = 0
00367 *
00368 *     Go to the tests for read-in matrix pairs
00369 *
00370       IFUNC = 0
00371       IF( NSIZE.EQ.0 )
00372      $   GO TO 70
00373 *
00374 *     Test the built-in matrix pairs.
00375 *     Loop over different functions (IFUNC) of DGGESX, types (PRTYPE)
00376 *     of test matrices, different size (M+N)
00377 *
00378       PRTYPE = 0
00379       QBA = 3
00380       QBB = 4
00381       WEIGHT = SQRT( ULP )
00382 *
00383       DO 60 IFUNC = 0, 3
00384          DO 50 PRTYPE = 1, 5
00385             DO 40 M = 1, NSIZE - 1
00386                DO 30 N = 1, NSIZE - M
00387 *
00388                   WEIGHT = ONE / WEIGHT
00389                   MPLUSN = M + N
00390 *
00391 *                 Generate test matrices
00392 *
00393                   FS = .TRUE.
00394                   K = 0
00395 *
00396                   CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI,
00397      $                         LDA )
00398                   CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI,
00399      $                         LDA )
00400 *
00401                   CALL DLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ),
00402      $                         LDA, AI( 1, M+1 ), LDA, BI, LDA,
00403      $                         BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA,
00404      $                         Q, LDA, Z, LDA, WEIGHT, QBA, QBB )
00405 *
00406 *                 Compute the Schur factorization and swapping the
00407 *                 m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
00408 *                 Swapping is accomplished via the function DLCTSX
00409 *                 which is supplied below.
00410 *
00411                   IF( IFUNC.EQ.0 ) THEN
00412                      SENSE = 'N'
00413                   ELSE IF( IFUNC.EQ.1 ) THEN
00414                      SENSE = 'E'
00415                   ELSE IF( IFUNC.EQ.2 ) THEN
00416                      SENSE = 'V'
00417                   ELSE IF( IFUNC.EQ.3 ) THEN
00418                      SENSE = 'B'
00419                   END IF
00420 *
00421                   CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
00422                   CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
00423 *
00424                   CALL DGGESX( 'V', 'V', 'S', DLCTSX, SENSE, MPLUSN, AI,
00425      $                         LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA,
00426      $                         Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK,
00427      $                         IWORK, LIWORK, BWORK, LINFO )
00428 *
00429                   IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
00430                      RESULT( 1 ) = ULPINV
00431                      WRITE( NOUT, FMT = 9999 )'DGGESX', LINFO, MPLUSN,
00432      $                  PRTYPE
00433                      INFO = LINFO
00434                      GO TO 30
00435                   END IF
00436 *
00437 *                 Compute the norm(A, B)
00438 *
00439                   CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
00440      $                         MPLUSN )
00441                   CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
00442      $                         WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
00443                   ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
00444      $                    WORK )
00445 *
00446 *                 Do tests (1) to (4)
00447 *
00448                   CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
00449      $                         LDA, WORK, RESULT( 1 ) )
00450                   CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
00451      $                         LDA, WORK, RESULT( 2 ) )
00452                   CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
00453      $                         LDA, WORK, RESULT( 3 ) )
00454                   CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
00455      $                         LDA, WORK, RESULT( 4 ) )
00456                   NTEST = 4
00457 *
00458 *                 Do tests (5) and (6): check Schur form of A and
00459 *                 compare eigenvalues with diagonals.
00460 *
00461                   TEMP1 = ZERO
00462                   RESULT( 5 ) = ZERO
00463                   RESULT( 6 ) = ZERO
00464 *
00465                   DO 10 J = 1, MPLUSN
00466                      ILABAD = .FALSE.
00467                      IF( ALPHAI( J ).EQ.ZERO ) THEN
00468                         TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
00469      $                          MAX( SMLNUM, ABS( ALPHAR( J ) ),
00470      $                          ABS( AI( J, J ) ) )+
00471      $                          ABS( BETA( J )-BI( J, J ) ) /
00472      $                          MAX( SMLNUM, ABS( BETA( J ) ),
00473      $                          ABS( BI( J, J ) ) ) ) / ULP
00474                         IF( J.LT.MPLUSN ) THEN
00475                            IF( AI( J+1, J ).NE.ZERO ) THEN
00476                               ILABAD = .TRUE.
00477                               RESULT( 5 ) = ULPINV
00478                            END IF
00479                         END IF
00480                         IF( J.GT.1 ) THEN
00481                            IF( AI( J, J-1 ).NE.ZERO ) THEN
00482                               ILABAD = .TRUE.
00483                               RESULT( 5 ) = ULPINV
00484                            END IF
00485                         END IF
00486                      ELSE
00487                         IF( ALPHAI( J ).GT.ZERO ) THEN
00488                            I1 = J
00489                         ELSE
00490                            I1 = J - 1
00491                         END IF
00492                         IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
00493                            ILABAD = .TRUE.
00494                         ELSE IF( I1.LT.MPLUSN-1 ) THEN
00495                            IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
00496                               ILABAD = .TRUE.
00497                               RESULT( 5 ) = ULPINV
00498                            END IF
00499                         ELSE IF( I1.GT.1 ) THEN
00500                            IF( AI( I1, I1-1 ).NE.ZERO ) THEN
00501                               ILABAD = .TRUE.
00502                               RESULT( 5 ) = ULPINV
00503                            END IF
00504                         END IF
00505                         IF( .NOT.ILABAD ) THEN
00506                            CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ),
00507      $                                  LDA, BETA( J ), ALPHAR( J ),
00508      $                                  ALPHAI( J ), TEMP2, IINFO )
00509                            IF( IINFO.GE.3 ) THEN
00510                               WRITE( NOUT, FMT = 9997 )IINFO, J,
00511      $                           MPLUSN, PRTYPE
00512                               INFO = ABS( IINFO )
00513                            END IF
00514                         ELSE
00515                            TEMP2 = ULPINV
00516                         END IF
00517                      END IF
00518                      TEMP1 = MAX( TEMP1, TEMP2 )
00519                      IF( ILABAD ) THEN
00520                         WRITE( NOUT, FMT = 9996 )J, MPLUSN, PRTYPE
00521                      END IF
00522    10             CONTINUE
00523                   RESULT( 6 ) = TEMP1
00524                   NTEST = NTEST + 2
00525 *
00526 *                 Test (7) (if sorting worked)
00527 *
00528                   RESULT( 7 ) = ZERO
00529                   IF( LINFO.EQ.MPLUSN+3 ) THEN
00530                      RESULT( 7 ) = ULPINV
00531                   ELSE IF( MM.NE.N ) THEN
00532                      RESULT( 7 ) = ULPINV
00533                   END IF
00534                   NTEST = NTEST + 1
00535 *
00536 *                 Test (8): compare the estimated value DIF and its
00537 *                 value. first, compute the exact DIF.
00538 *
00539                   RESULT( 8 ) = ZERO
00540                   MN2 = MM*( MPLUSN-MM )*2
00541                   IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
00542 *
00543 *                    Note: for either following two causes, there are
00544 *                    almost same number of test cases fail the test.
00545 *
00546                      CALL DLAKF2( MM, MPLUSN-MM, AI, LDA,
00547      $                            AI( MM+1, MM+1 ), BI,
00548      $                            BI( MM+1, MM+1 ), C, LDC )
00549 *
00550                      CALL DGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
00551      $                            1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
00552      $                            INFO )
00553                      DIFTRU = S( MN2 )
00554 *
00555                      IF( DIFEST( 2 ).EQ.ZERO ) THEN
00556                         IF( DIFTRU.GT.ABNRM*ULP )
00557      $                     RESULT( 8 ) = ULPINV
00558                      ELSE IF( DIFTRU.EQ.ZERO ) THEN
00559                         IF( DIFEST( 2 ).GT.ABNRM*ULP )
00560      $                     RESULT( 8 ) = ULPINV
00561                      ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
00562      $                        ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
00563                         RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ),
00564      $                                DIFEST( 2 ) / DIFTRU )
00565                      END IF
00566                      NTEST = NTEST + 1
00567                   END IF
00568 *
00569 *                 Test (9)
00570 *
00571                   RESULT( 9 ) = ZERO
00572                   IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
00573                      IF( DIFTRU.GT.ABNRM*ULP )
00574      $                  RESULT( 9 ) = ULPINV
00575                      IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
00576      $                  RESULT( 9 ) = ULPINV
00577                      IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
00578      $                  RESULT( 9 ) = ULPINV
00579                      NTEST = NTEST + 1
00580                   END IF
00581 *
00582                   NTESTT = NTESTT + NTEST
00583 *
00584 *                 Print out tests which fail.
00585 *
00586                   DO 20 J = 1, 9
00587                      IF( RESULT( J ).GE.THRESH ) THEN
00588 *
00589 *                       If this is the first test to fail,
00590 *                       print a header to the data file.
00591 *
00592                         IF( NERRS.EQ.0 ) THEN
00593                            WRITE( NOUT, FMT = 9995 )'SGX'
00594 *
00595 *                          Matrix types
00596 *
00597                            WRITE( NOUT, FMT = 9993 )
00598 *
00599 *                          Tests performed
00600 *
00601                            WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
00602      $                        'transpose', ( '''', I = 1, 4 )
00603 *
00604                         END IF
00605                         NERRS = NERRS + 1
00606                         IF( RESULT( J ).LT.10000.0D0 ) THEN
00607                            WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
00608      $                        WEIGHT, M, J, RESULT( J )
00609                         ELSE
00610                            WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE,
00611      $                        WEIGHT, M, J, RESULT( J )
00612                         END IF
00613                      END IF
00614    20             CONTINUE
00615 *
00616    30          CONTINUE
00617    40       CONTINUE
00618    50    CONTINUE
00619    60 CONTINUE
00620 *
00621       GO TO 150
00622 *
00623    70 CONTINUE
00624 *
00625 *     Read in data from file to check accuracy of condition estimation
00626 *     Read input data until N=0
00627 *
00628       NPTKNT = 0
00629 *
00630    80 CONTINUE
00631       READ( NIN, FMT = *, END = 140 )MPLUSN
00632       IF( MPLUSN.EQ.0 )
00633      $   GO TO 140
00634       READ( NIN, FMT = *, END = 140 )N
00635       DO 90 I = 1, MPLUSN
00636          READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
00637    90 CONTINUE
00638       DO 100 I = 1, MPLUSN
00639          READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
00640   100 CONTINUE
00641       READ( NIN, FMT = * )PLTRU, DIFTRU
00642 *
00643       NPTKNT = NPTKNT + 1
00644       FS = .TRUE.
00645       K = 0
00646       M = MPLUSN - N
00647 *
00648       CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
00649       CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
00650 *
00651 *     Compute the Schur factorization while swaping the
00652 *     m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
00653 *
00654       CALL DGGESX( 'V', 'V', 'S', DLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
00655      $             MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST,
00656      $             WORK, LWORK, IWORK, LIWORK, BWORK, LINFO )
00657 *
00658       IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
00659          RESULT( 1 ) = ULPINV
00660          WRITE( NOUT, FMT = 9998 )'DGGESX', LINFO, MPLUSN, NPTKNT
00661          GO TO 130
00662       END IF
00663 *
00664 *     Compute the norm(A, B)
00665 *        (should this be norm of (A,B) or (AI,BI)?)
00666 *
00667       CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
00668       CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
00669      $             WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
00670       ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK )
00671 *
00672 *     Do tests (1) to (4)
00673 *
00674       CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
00675      $             RESULT( 1 ) )
00676       CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
00677      $             RESULT( 2 ) )
00678       CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
00679      $             RESULT( 3 ) )
00680       CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
00681      $             RESULT( 4 ) )
00682 *
00683 *     Do tests (5) and (6): check Schur form of A and compare
00684 *     eigenvalues with diagonals.
00685 *
00686       NTEST = 6
00687       TEMP1 = ZERO
00688       RESULT( 5 ) = ZERO
00689       RESULT( 6 ) = ZERO
00690 *
00691       DO 110 J = 1, MPLUSN
00692          ILABAD = .FALSE.
00693          IF( ALPHAI( J ).EQ.ZERO ) THEN
00694             TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
00695      $              MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J,
00696      $              J ) ) )+ABS( BETA( J )-BI( J, J ) ) /
00697      $              MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) )
00698      $               / ULP
00699             IF( J.LT.MPLUSN ) THEN
00700                IF( AI( J+1, J ).NE.ZERO ) THEN
00701                   ILABAD = .TRUE.
00702                   RESULT( 5 ) = ULPINV
00703                END IF
00704             END IF
00705             IF( J.GT.1 ) THEN
00706                IF( AI( J, J-1 ).NE.ZERO ) THEN
00707                   ILABAD = .TRUE.
00708                   RESULT( 5 ) = ULPINV
00709                END IF
00710             END IF
00711          ELSE
00712             IF( ALPHAI( J ).GT.ZERO ) THEN
00713                I1 = J
00714             ELSE
00715                I1 = J - 1
00716             END IF
00717             IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
00718                ILABAD = .TRUE.
00719             ELSE IF( I1.LT.MPLUSN-1 ) THEN
00720                IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
00721                   ILABAD = .TRUE.
00722                   RESULT( 5 ) = ULPINV
00723                END IF
00724             ELSE IF( I1.GT.1 ) THEN
00725                IF( AI( I1, I1-1 ).NE.ZERO ) THEN
00726                   ILABAD = .TRUE.
00727                   RESULT( 5 ) = ULPINV
00728                END IF
00729             END IF
00730             IF( .NOT.ILABAD ) THEN
00731                CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA,
00732      $                      BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2,
00733      $                      IINFO )
00734                IF( IINFO.GE.3 ) THEN
00735                   WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT
00736                   INFO = ABS( IINFO )
00737                END IF
00738             ELSE
00739                TEMP2 = ULPINV
00740             END IF
00741          END IF
00742          TEMP1 = MAX( TEMP1, TEMP2 )
00743          IF( ILABAD ) THEN
00744             WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT
00745          END IF
00746   110 CONTINUE
00747       RESULT( 6 ) = TEMP1
00748 *
00749 *     Test (7) (if sorting worked)  <--------- need to be checked.
00750 *
00751       NTEST = 7
00752       RESULT( 7 ) = ZERO
00753       IF( LINFO.EQ.MPLUSN+3 )
00754      $   RESULT( 7 ) = ULPINV
00755 *
00756 *     Test (8): compare the estimated value of DIF and its true value.
00757 *
00758       NTEST = 8
00759       RESULT( 8 ) = ZERO
00760       IF( DIFEST( 2 ).EQ.ZERO ) THEN
00761          IF( DIFTRU.GT.ABNRM*ULP )
00762      $      RESULT( 8 ) = ULPINV
00763       ELSE IF( DIFTRU.EQ.ZERO ) THEN
00764          IF( DIFEST( 2 ).GT.ABNRM*ULP )
00765      $      RESULT( 8 ) = ULPINV
00766       ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
00767      $         ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
00768          RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
00769       END IF
00770 *
00771 *     Test (9)
00772 *
00773       NTEST = 9
00774       RESULT( 9 ) = ZERO
00775       IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
00776          IF( DIFTRU.GT.ABNRM*ULP )
00777      $      RESULT( 9 ) = ULPINV
00778          IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
00779      $      RESULT( 9 ) = ULPINV
00780          IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
00781      $      RESULT( 9 ) = ULPINV
00782       END IF
00783 *
00784 *     Test (10): compare the estimated value of PL and it true value.
00785 *
00786       NTEST = 10
00787       RESULT( 10 ) = ZERO
00788       IF( PL( 1 ).EQ.ZERO ) THEN
00789          IF( PLTRU.GT.ABNRM*ULP )
00790      $      RESULT( 10 ) = ULPINV
00791       ELSE IF( PLTRU.EQ.ZERO ) THEN
00792          IF( PL( 1 ).GT.ABNRM*ULP )
00793      $      RESULT( 10 ) = ULPINV
00794       ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
00795      $         ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
00796          RESULT( 10 ) = ULPINV
00797       END IF
00798 *
00799       NTESTT = NTESTT + NTEST
00800 *
00801 *     Print out tests which fail.
00802 *
00803       DO 120 J = 1, NTEST
00804          IF( RESULT( J ).GE.THRESH ) THEN
00805 *
00806 *           If this is the first test to fail,
00807 *           print a header to the data file.
00808 *
00809             IF( NERRS.EQ.0 ) THEN
00810                WRITE( NOUT, FMT = 9995 )'SGX'
00811 *
00812 *              Matrix types
00813 *
00814                WRITE( NOUT, FMT = 9994 )
00815 *
00816 *              Tests performed
00817 *
00818                WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
00819      $            'transpose', ( '''', I = 1, 4 )
00820 *
00821             END IF
00822             NERRS = NERRS + 1
00823             IF( RESULT( J ).LT.10000.0D0 ) THEN
00824                WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
00825             ELSE
00826                WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J )
00827             END IF
00828          END IF
00829 *
00830   120 CONTINUE
00831 *
00832   130 CONTINUE
00833       GO TO 80
00834   140 CONTINUE
00835 *
00836   150 CONTINUE
00837 *
00838 *     Summary
00839 *
00840       CALL ALASVM( 'SGX', NOUT, NERRS, NTESTT, 0 )
00841 *
00842       WORK( 1 ) = MAXWRK
00843 *
00844       RETURN
00845 *
00846  9999 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00847      $      I6, ', JTYPE=', I6, ')' )
00848 *
00849  9998 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00850      $      I6, ', Input Example #', I2, ')' )
00851 *
00852  9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', I1, ' for eigenvalue ',
00853      $      I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
00854 *
00855  9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', I6, '.',
00856      $      / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
00857 *
00858  9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form',
00859      $      ' problem driver' )
00860 *
00861  9994 FORMAT( 'Input Example' )
00862 *
00863  9993 FORMAT( ' Matrix types: ', /
00864      $      '  1:  A is a block diagonal matrix of Jordan blocks ',
00865      $      'and B is the identity ', / '      matrix, ',
00866      $      / '  2:  A and B are upper triangular matrices, ',
00867      $      / '  3:  A and B are as type 2, but each second diagonal ',
00868      $      'block in A_11 and ', /
00869      $      '      each third diaongal block in A_22 are 2x2 blocks,',
00870      $      / '  4:  A and B are block diagonal matrices, ',
00871      $      / '  5:  (A,B) has potentially close or common ',
00872      $      'eigenvalues.', / )
00873 *
00874  9992 FORMAT( / ' Tests performed:  (S is Schur, T is triangular, ',
00875      $      'Q and Z are ', A, ',', / 19X,
00876      $      ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
00877      $      / '  1 = | A - Q S Z', A,
00878      $      ' | / ( |A| n ulp )      2 = | B - Q T Z', A,
00879      $      ' | / ( |B| n ulp )', / '  3 = | I - QQ', A,
00880      $      ' | / ( n ulp )             4 = | I - ZZ', A,
00881      $      ' | / ( n ulp )', / '  5 = 1/ULP  if A is not in ',
00882      $      'Schur form S', / '  6 = difference between (alpha,beta)',
00883      $      ' and diagonals of (S,T)', /
00884      $      '  7 = 1/ULP  if SDIM is not the correct number of ',
00885      $      'selected eigenvalues', /
00886      $      '  8 = 1/ULP  if DIFEST/DIFTRU > 10*THRESH or ',
00887      $      'DIFTRU/DIFEST > 10*THRESH',
00888      $      / '  9 = 1/ULP  if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
00889      $      'when reordering fails', /
00890      $      ' 10 = 1/ULP  if PLEST/PLTRU > THRESH or ',
00891      $      'PLTRU/PLEST > THRESH', /
00892      $      '    ( Test 10 is only for input examples )', / )
00893  9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.4,
00894      $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
00895  9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.4,
00896      $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, D10.4 )
00897  9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00898      $      ' result ', I2, ' is', 0P, F8.2 )
00899  9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00900      $      ' result ', I2, ' is', 1P, D10.3 )
00901 *
00902 *     End of DDRGSX
00903 *
00904       END