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