LAPACK: TESTING/EIG/zdrgsx.f Source File

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