LAPACK: TESTING/EIG/cdrgsx.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
00002      $                   AI, BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK,
00003      $                   LWORK, RWORK, 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               RWORK( * ), S( * )
00018       COMPLEX            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 *  CDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
00028 *  problem expert driver CGGESX.
00029 *
00030 *  CGGES 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 CDRGSX is called with NSIZE > 0, five (5) types of built-in
00046 *  matrix pairs are used to test the routine CGGESX.
00047 *
00048 *  When CDRGSX is called with NSIZE = 0, it reads in test matrix data
00049 *  to test CGGESX.
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 CGGESX, 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 CGESVD) 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 SLATM5).
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 SGGESX.
00170 *
00171 *  NCMAX   (input) INTEGER
00172 *          Maximum allowable NMAX for generating Kroneker matrix
00173 *          in call to CLAKF2
00174 *
00175 *  THRESH  (input) REAL
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 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 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 array, dimension (LDA, NSIZE)
00205 *          Copy of A, modified by CGGESX.
00206 *
00207 *  BI      (workspace) COMPLEX array, dimension (LDA, NSIZE)
00208 *          Copy of B, modified by CGGESX.
00209 *
00210 *  Z       (workspace) COMPLEX array, dimension (LDA, NSIZE)
00211 *          Z holds the left Schur vectors computed by CGGESX.
00212 *
00213 *  Q       (workspace) COMPLEX array, dimension (LDA, NSIZE)
00214 *          Q holds the right Schur vectors computed by CGGESX.
00215 *
00216 *  ALPHA   (workspace) COMPLEX array, dimension (NSIZE)
00217 *  BETA    (workspace) COMPLEX array, dimension (NSIZE)
00218 *          On exit, ALPHA/BETA are the eigenvalues.
00219 *
00220 *  C       (workspace) COMPLEX array, dimension (LDC, LDC)
00221 *          Store the matrix generated by subroutine CLAKF2, 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) REAL array, dimension (LDC)
00229 *          Singular values of C
00230 *
00231 *  WORK    (workspace) COMPLEX array, dimension (LWORK)
00232 *
00233 *  LWORK   (input) INTEGER
00234 *          The dimension of the array WORK.  LWORK >= 3*NSIZE*NSIZE/2
00235 *
00236 *  RWORK   (workspace) REAL 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       REAL               ZERO, ONE, TEN
00255       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1 )
00256       COMPLEX            CZERO
00257       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+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       REAL               ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
00266      $                   TEMP2, THRSH2, ULP, ULPINV, WEIGHT
00267       COMPLEX            X
00268 *     ..
00269 *     .. Local Arrays ..
00270       REAL               DIFEST( 2 ), PL( 2 ), RESULT( 10 )
00271 *     ..
00272 *     .. External Functions ..
00273       LOGICAL            CLCTSX
00274       INTEGER            ILAENV
00275       REAL               CLANGE, SLAMCH
00276       EXTERNAL           CLCTSX, ILAENV, CLANGE, SLAMCH
00277 *     ..
00278 *     .. External Subroutines ..
00279       EXTERNAL           ALASVM, CGESVD, CGET51, CGGESX, CLACPY, CLAKF2,
00280      $                   CLASET, CLATM5, SLABAD, XERBLA
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, AIMAG, MAX, REAL, SQRT
00291 *     ..
00292 *     .. Statement Functions ..
00293       REAL               ABS1
00294 *     ..
00295 *     .. Statement Function definitions ..
00296       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
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 = -15
00314       ELSE IF( LIWORK.LT.NSIZE+2 ) 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 = 3*NSIZE*NSIZE / 2
00328 *
00329 *        workspace for cggesx
00330 *
00331          MAXWRK = NSIZE*( 1+ILAENV( 1, 'CGEQRF', ' ', NSIZE, 1, NSIZE,
00332      $            0 ) )
00333          MAXWRK = MAX( MAXWRK, NSIZE*( 1+ILAENV( 1, 'CUNGQR', ' ',
00334      $            NSIZE, 1, NSIZE, -1 ) ) )
00335 *
00336 *        workspace for cgesvd
00337 *
00338          BDSPAC = 3*NSIZE*NSIZE / 2
00339          MAXWRK = MAX( MAXWRK, NSIZE*NSIZE*
00340      $            ( 1+ILAENV( 1, 'CGEBRD', ' ', 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 = -18
00351 *
00352       IF( INFO.NE.0 ) THEN
00353          CALL XERBLA( 'CDRGSX', -INFO )
00354          RETURN
00355       END IF
00356 *
00357 *     Important constants
00358 *
00359       ULP = SLAMCH( 'P' )
00360       ULPINV = ONE / ULP
00361       SMLNUM = SLAMCH( 'S' ) / ULP
00362       BIGNUM = ONE / SMLNUM
00363       CALL SLABAD( 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 CGGESX, 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 CLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, AI,
00397      $                         LDA )
00398                   CALL CLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, BI,
00399      $                         LDA )
00400 *
00401                   CALL CLATM5( 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 CLCTSX
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 CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
00422                   CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
00423 *
00424                   CALL CGGESX( 'V', 'V', 'S', CLCTSX, SENSE, MPLUSN, AI,
00425      $                         LDA, BI, LDA, MM, ALPHA, BETA, Q, LDA, Z,
00426      $                         LDA, PL, DIFEST, WORK, LWORK, RWORK,
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 )'CGGESX', LINFO, MPLUSN,
00432      $                  PRTYPE
00433                      INFO = LINFO
00434                      GO TO 30
00435                   END IF
00436 *
00437 *                 Compute the norm(A, B)
00438 *
00439                   CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
00440      $                         MPLUSN )
00441                   CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
00442      $                         WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
00443                   ABNRM = CLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
00444      $                    RWORK )
00445 *
00446 *                 Do tests (1) to (4)
00447 *
00448                   RESULT( 2 ) = ZERO
00449                   CALL CGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
00450      $                         LDA, WORK, RWORK, RESULT( 1 ) )
00451                   CALL CGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
00452      $                         LDA, WORK, RWORK, RESULT( 2 ) )
00453                   CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
00454      $                         LDA, WORK, RWORK, RESULT( 3 ) )
00455                   CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
00456      $                         LDA, WORK, RWORK, 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                      TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) /
00469      $                       MAX( SMLNUM, ABS1( ALPHA( J ) ),
00470      $                       ABS1( AI( J, J ) ) )+
00471      $                       ABS1( BETA( J )-BI( J, J ) ) /
00472      $                       MAX( SMLNUM, ABS1( BETA( J ) ),
00473      $                       ABS1( 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                      TEMP1 = MAX( TEMP1, TEMP2 )
00487                      IF( ILABAD ) THEN
00488                         WRITE( NOUT, FMT = 9997 )J, MPLUSN, PRTYPE
00489                      END IF
00490    10             CONTINUE
00491                   RESULT( 6 ) = TEMP1
00492                   NTEST = NTEST + 2
00493 *
00494 *                 Test (7) (if sorting worked)
00495 *
00496                   RESULT( 7 ) = ZERO
00497                   IF( LINFO.EQ.MPLUSN+3 ) THEN
00498                      RESULT( 7 ) = ULPINV
00499                   ELSE IF( MM.NE.N ) THEN
00500                      RESULT( 7 ) = ULPINV
00501                   END IF
00502                   NTEST = NTEST + 1
00503 *
00504 *                 Test (8): compare the estimated value DIF and its
00505 *                 value. first, compute the exact DIF.
00506 *
00507                   RESULT( 8 ) = ZERO
00508                   MN2 = MM*( MPLUSN-MM )*2
00509                   IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
00510 *
00511 *                    Note: for either following two cases, there are
00512 *                    almost same number of test cases fail the test.
00513 *
00514                      CALL CLAKF2( MM, MPLUSN-MM, AI, LDA,
00515      $                            AI( MM+1, MM+1 ), BI,
00516      $                            BI( MM+1, MM+1 ), C, LDC )
00517 *
00518                      CALL CGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
00519      $                            1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
00520      $                            RWORK, INFO )
00521                      DIFTRU = S( MN2 )
00522 *
00523                      IF( DIFEST( 2 ).EQ.ZERO ) THEN
00524                         IF( DIFTRU.GT.ABNRM*ULP )
00525      $                     RESULT( 8 ) = ULPINV
00526                      ELSE IF( DIFTRU.EQ.ZERO ) THEN
00527                         IF( DIFEST( 2 ).GT.ABNRM*ULP )
00528      $                     RESULT( 8 ) = ULPINV
00529                      ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
00530      $                        ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
00531                         RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ),
00532      $                                DIFEST( 2 ) / DIFTRU )
00533                      END IF
00534                      NTEST = NTEST + 1
00535                   END IF
00536 *
00537 *                 Test (9)
00538 *
00539                   RESULT( 9 ) = ZERO
00540                   IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
00541                      IF( DIFTRU.GT.ABNRM*ULP )
00542      $                  RESULT( 9 ) = ULPINV
00543                      IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
00544      $                  RESULT( 9 ) = ULPINV
00545                      IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
00546      $                  RESULT( 9 ) = ULPINV
00547                      NTEST = NTEST + 1
00548                   END IF
00549 *
00550                   NTESTT = NTESTT + NTEST
00551 *
00552 *                 Print out tests which fail.
00553 *
00554                   DO 20 J = 1, 9
00555                      IF( RESULT( J ).GE.THRESH ) THEN
00556 *
00557 *                       If this is the first test to fail,
00558 *                       print a header to the data file.
00559 *
00560                         IF( NERRS.EQ.0 ) THEN
00561                            WRITE( NOUT, FMT = 9996 )'CGX'
00562 *
00563 *                          Matrix types
00564 *
00565                            WRITE( NOUT, FMT = 9994 )
00566 *
00567 *                          Tests performed
00568 *
00569                            WRITE( NOUT, FMT = 9993 )'unitary', '''',
00570      $                        'transpose', ( '''', I = 1, 4 )
00571 *
00572                         END IF
00573                         NERRS = NERRS + 1
00574                         IF( RESULT( J ).LT.10000.0 ) THEN
00575                            WRITE( NOUT, FMT = 9992 )MPLUSN, PRTYPE,
00576      $                        WEIGHT, M, J, RESULT( J )
00577                         ELSE
00578                            WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
00579      $                        WEIGHT, M, J, RESULT( J )
00580                         END IF
00581                      END IF
00582    20             CONTINUE
00583 *
00584    30          CONTINUE
00585    40       CONTINUE
00586    50    CONTINUE
00587    60 CONTINUE
00588 *
00589       GO TO 150
00590 *
00591    70 CONTINUE
00592 *
00593 *     Read in data from file to check accuracy of condition estimation
00594 *     Read input data until N=0
00595 *
00596       NPTKNT = 0
00597 *
00598    80 CONTINUE
00599       READ( NIN, FMT = *, END = 140 )MPLUSN
00600       IF( MPLUSN.EQ.0 )
00601      $   GO TO 140
00602       READ( NIN, FMT = *, END = 140 )N
00603       DO 90 I = 1, MPLUSN
00604          READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
00605    90 CONTINUE
00606       DO 100 I = 1, MPLUSN
00607          READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
00608   100 CONTINUE
00609       READ( NIN, FMT = * )PLTRU, DIFTRU
00610 *
00611       NPTKNT = NPTKNT + 1
00612       FS = .TRUE.
00613       K = 0
00614       M = MPLUSN - N
00615 *
00616       CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
00617       CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
00618 *
00619 *     Compute the Schur factorization while swaping the
00620 *     m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
00621 *
00622       CALL CGGESX( 'V', 'V', 'S', CLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
00623      $             MM, ALPHA, BETA, Q, LDA, Z, LDA, PL, DIFEST, WORK,
00624      $             LWORK, RWORK, IWORK, LIWORK, BWORK, LINFO )
00625 *
00626       IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
00627          RESULT( 1 ) = ULPINV
00628          WRITE( NOUT, FMT = 9998 )'CGGESX', LINFO, MPLUSN, NPTKNT
00629          GO TO 130
00630       END IF
00631 *
00632 *     Compute the norm(A, B)
00633 *        (should this be norm of (A,B) or (AI,BI)?)
00634 *
00635       CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
00636       CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
00637      $             WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
00638       ABNRM = CLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, RWORK )
00639 *
00640 *     Do tests (1) to (4)
00641 *
00642       CALL CGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
00643      $             RWORK, RESULT( 1 ) )
00644       CALL CGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
00645      $             RWORK, RESULT( 2 ) )
00646       CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
00647      $             RWORK, RESULT( 3 ) )
00648       CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
00649      $             RWORK, RESULT( 4 ) )
00650 *
00651 *     Do tests (5) and (6): check Schur form of A and compare
00652 *     eigenvalues with diagonals.
00653 *
00654       NTEST = 6
00655       TEMP1 = ZERO
00656       RESULT( 5 ) = ZERO
00657       RESULT( 6 ) = ZERO
00658 *
00659       DO 110 J = 1, MPLUSN
00660          ILABAD = .FALSE.
00661          TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) /
00662      $           MAX( SMLNUM, ABS1( ALPHA( J ) ), ABS1( AI( J, J ) ) )+
00663      $           ABS1( BETA( J )-BI( J, J ) ) /
00664      $           MAX( SMLNUM, ABS1( BETA( J ) ), ABS1( BI( J, J ) ) ) )
00665      $            / ULP
00666          IF( J.LT.MPLUSN ) THEN
00667             IF( AI( J+1, J ).NE.ZERO ) THEN
00668                ILABAD = .TRUE.
00669                RESULT( 5 ) = ULPINV
00670             END IF
00671          END IF
00672          IF( J.GT.1 ) THEN
00673             IF( AI( J, J-1 ).NE.ZERO ) THEN
00674                ILABAD = .TRUE.
00675                RESULT( 5 ) = ULPINV
00676             END IF
00677          END IF
00678          TEMP1 = MAX( TEMP1, TEMP2 )
00679          IF( ILABAD ) THEN
00680             WRITE( NOUT, FMT = 9997 )J, MPLUSN, NPTKNT
00681          END IF
00682   110 CONTINUE
00683       RESULT( 6 ) = TEMP1
00684 *
00685 *     Test (7) (if sorting worked)  <--------- need to be checked.
00686 *
00687       NTEST = 7
00688       RESULT( 7 ) = ZERO
00689       IF( LINFO.EQ.MPLUSN+3 )
00690      $   RESULT( 7 ) = ULPINV
00691 *
00692 *     Test (8): compare the estimated value of DIF and its true value.
00693 *
00694       NTEST = 8
00695       RESULT( 8 ) = ZERO
00696       IF( DIFEST( 2 ).EQ.ZERO ) THEN
00697          IF( DIFTRU.GT.ABNRM*ULP )
00698      $      RESULT( 8 ) = ULPINV
00699       ELSE IF( DIFTRU.EQ.ZERO ) THEN
00700          IF( DIFEST( 2 ).GT.ABNRM*ULP )
00701      $      RESULT( 8 ) = ULPINV
00702       ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
00703      $         ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
00704          RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
00705       END IF
00706 *
00707 *     Test (9)
00708 *
00709       NTEST = 9
00710       RESULT( 9 ) = ZERO
00711       IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
00712          IF( DIFTRU.GT.ABNRM*ULP )
00713      $      RESULT( 9 ) = ULPINV
00714          IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
00715      $      RESULT( 9 ) = ULPINV
00716          IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
00717      $      RESULT( 9 ) = ULPINV
00718       END IF
00719 *
00720 *     Test (10): compare the estimated value of PL and it true value.
00721 *
00722       NTEST = 10
00723       RESULT( 10 ) = ZERO
00724       IF( PL( 1 ).EQ.ZERO ) THEN
00725          IF( PLTRU.GT.ABNRM*ULP )
00726      $      RESULT( 10 ) = ULPINV
00727       ELSE IF( PLTRU.EQ.ZERO ) THEN
00728          IF( PL( 1 ).GT.ABNRM*ULP )
00729      $      RESULT( 10 ) = ULPINV
00730       ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
00731      $         ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
00732          RESULT( 10 ) = ULPINV
00733       END IF
00734 *
00735       NTESTT = NTESTT + NTEST
00736 *
00737 *     Print out tests which fail.
00738 *
00739       DO 120 J = 1, NTEST
00740          IF( RESULT( J ).GE.THRESH ) THEN
00741 *
00742 *           If this is the first test to fail,
00743 *           print a header to the data file.
00744 *
00745             IF( NERRS.EQ.0 ) THEN
00746                WRITE( NOUT, FMT = 9996 )'CGX'
00747 *
00748 *              Matrix types
00749 *
00750                WRITE( NOUT, FMT = 9995 )
00751 *
00752 *              Tests performed
00753 *
00754                WRITE( NOUT, FMT = 9993 )'unitary', '''', 'transpose',
00755      $            ( '''', I = 1, 4 )
00756 *
00757             END IF
00758             NERRS = NERRS + 1
00759             IF( RESULT( J ).LT.10000.0 ) THEN
00760                WRITE( NOUT, FMT = 9990 )NPTKNT, MPLUSN, J, RESULT( J )
00761             ELSE
00762                WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
00763             END IF
00764          END IF
00765 *
00766   120 CONTINUE
00767 *
00768   130 CONTINUE
00769       GO TO 80
00770   140 CONTINUE
00771 *
00772   150 CONTINUE
00773 *
00774 *     Summary
00775 *
00776       CALL ALASVM( 'CGX', NOUT, NERRS, NTESTT, 0 )
00777 *
00778       WORK( 1 ) = MAXWRK
00779 *
00780       RETURN
00781 *
00782  9999 FORMAT( ' CDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00783      $      I6, ', JTYPE=', I6, ')' )
00784 *
00785  9998 FORMAT( ' CDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00786      $      I6, ', Input Example #', I2, ')' )
00787 *
00788  9997 FORMAT( ' CDRGSX: S not in Schur form at eigenvalue ', I6, '.',
00789      $      / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
00790 *
00791  9996 FORMAT( / 1X, A3, ' -- Complex Expert Generalized Schur form',
00792      $      ' problem driver' )
00793 *
00794  9995 FORMAT( 'Input Example' )
00795 *
00796  9994 FORMAT( ' Matrix types: ', /
00797      $      '  1:  A is a block diagonal matrix of Jordan blocks ',
00798      $      'and B is the identity ', / '      matrix, ',
00799      $      / '  2:  A and B are upper triangular matrices, ',
00800      $      / '  3:  A and B are as type 2, but each second diagonal ',
00801      $      'block in A_11 and ', /
00802      $      '      each third diaongal block in A_22 are 2x2 blocks,',
00803      $      / '  4:  A and B are block diagonal matrices, ',
00804      $      / '  5:  (A,B) has potentially close or common ',
00805      $      'eigenvalues.', / )
00806 *
00807  9993 FORMAT( / ' Tests performed:  (S is Schur, T is triangular, ',
00808      $      'Q and Z are ', A, ',', / 19X,
00809      $      ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
00810      $      / '  1 = | A - Q S Z', A,
00811      $      ' | / ( |A| n ulp )      2 = | B - Q T Z', A,
00812      $      ' | / ( |B| n ulp )', / '  3 = | I - QQ', A,
00813      $      ' | / ( n ulp )             4 = | I - ZZ', A,
00814      $      ' | / ( n ulp )', / '  5 = 1/ULP  if A is not in ',
00815      $      'Schur form S', / '  6 = difference between (alpha,beta)',
00816      $      ' and diagonals of (S,T)', /
00817      $      '  7 = 1/ULP  if SDIM is not the correct number of ',
00818      $      'selected eigenvalues', /
00819      $      '  8 = 1/ULP  if DIFEST/DIFTRU > 10*THRESH or ',
00820      $      'DIFTRU/DIFEST > 10*THRESH',
00821      $      / '  9 = 1/ULP  if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
00822      $      'when reordering fails', /
00823      $      ' 10 = 1/ULP  if PLEST/PLTRU > THRESH or ',
00824      $      'PLTRU/PLEST > THRESH', /
00825      $      '    ( Test 10 is only for input examples )', / )
00826  9992 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.4,
00827      $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
00828  9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.4,
00829      $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, E10.4 )
00830  9990 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00831      $      ' result ', I2, ' is', 0P, F8.2 )
00832  9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00833      $      ' result ', I2, ' is', 1P, E10.3 )
00834 *
00835 *     End of CDRGSX
00836 *
00837       END