LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zdrvsx.f
Go to the documentation of this file.
00001 *> \brief \b ZDRVSX
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 ZDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00012 *                          NIUNIT, NOUNIT, A, LDA, H, HT, W, WT, WTMP, VS,
00013 *                          LDVS, VS1, RESULT, WORK, LWORK, RWORK, BWORK,
00014 *                          INFO )
00015 * 
00016 *       .. Scalar Arguments ..
00017 *       INTEGER            INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
00018 *      $                   NTYPES
00019 *       DOUBLE PRECISION   THRESH
00020 *       ..
00021 *       .. Array Arguments ..
00022 *       LOGICAL            BWORK( * ), DOTYPE( * )
00023 *       INTEGER            ISEED( 4 ), NN( * )
00024 *       DOUBLE PRECISION   RESULT( 17 ), RWORK( * )
00025 *       COMPLEX*16         A( LDA, * ), H( LDA, * ), HT( LDA, * ),
00026 *      $                   VS( LDVS, * ), VS1( LDVS, * ), W( * ),
00027 *      $                   WORK( * ), WT( * ), WTMP( * )
00028 *       ..
00029 *  
00030 *
00031 *> \par Purpose:
00032 *  =============
00033 *>
00034 *> \verbatim
00035 *>
00036 *>    ZDRVSX checks the nonsymmetric eigenvalue (Schur form) problem
00037 *>    expert driver ZGEESX.
00038 *>
00039 *>    ZDRVSX uses both test matrices generated randomly depending on
00040 *>    data supplied in the calling sequence, as well as on data
00041 *>    read from an input file and including precomputed condition
00042 *>    numbers to which it compares the ones it computes.
00043 *>
00044 *>    When ZDRVSX is called, a number of matrix "sizes" ("n's") and a
00045 *>    number of matrix "types" are specified.  For each size ("n")
00046 *>    and each type of matrix, one matrix will be generated and used
00047 *>    to test the nonsymmetric eigenroutines.  For each matrix, 15
00048 *>    tests will be performed:
00049 *>
00050 *>    (1)     0 if T is in Schur form, 1/ulp otherwise
00051 *>           (no sorting of eigenvalues)
00052 *>
00053 *>    (2)     | A - VS T VS' | / ( n |A| ulp )
00054 *>
00055 *>      Here VS is the matrix of Schur eigenvectors, and T is in Schur
00056 *>      form  (no sorting of eigenvalues).
00057 *>
00058 *>    (3)     | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
00059 *>
00060 *>    (4)     0     if W are eigenvalues of T
00061 *>            1/ulp otherwise
00062 *>            (no sorting of eigenvalues)
00063 *>
00064 *>    (5)     0     if T(with VS) = T(without VS),
00065 *>            1/ulp otherwise
00066 *>            (no sorting of eigenvalues)
00067 *>
00068 *>    (6)     0     if eigenvalues(with VS) = eigenvalues(without VS),
00069 *>            1/ulp otherwise
00070 *>            (no sorting of eigenvalues)
00071 *>
00072 *>    (7)     0 if T is in Schur form, 1/ulp otherwise
00073 *>            (with sorting of eigenvalues)
00074 *>
00075 *>    (8)     | A - VS T VS' | / ( n |A| ulp )
00076 *>
00077 *>      Here VS is the matrix of Schur eigenvectors, and T is in Schur
00078 *>      form  (with sorting of eigenvalues).
00079 *>
00080 *>    (9)     | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
00081 *>
00082 *>    (10)    0     if W are eigenvalues of T
00083 *>            1/ulp otherwise
00084 *>            If workspace sufficient, also compare W with and
00085 *>            without reciprocal condition numbers
00086 *>            (with sorting of eigenvalues)
00087 *>
00088 *>    (11)    0     if T(with VS) = T(without VS),
00089 *>            1/ulp otherwise
00090 *>            If workspace sufficient, also compare T with and without
00091 *>            reciprocal condition numbers
00092 *>            (with sorting of eigenvalues)
00093 *>
00094 *>    (12)    0     if eigenvalues(with VS) = eigenvalues(without VS),
00095 *>            1/ulp otherwise
00096 *>            If workspace sufficient, also compare VS with and without
00097 *>            reciprocal condition numbers
00098 *>            (with sorting of eigenvalues)
00099 *>
00100 *>    (13)    if sorting worked and SDIM is the number of
00101 *>            eigenvalues which were SELECTed
00102 *>            If workspace sufficient, also compare SDIM with and
00103 *>            without reciprocal condition numbers
00104 *>
00105 *>    (14)    if RCONDE the same no matter if VS and/or RCONDV computed
00106 *>
00107 *>    (15)    if RCONDV the same no matter if VS and/or RCONDE computed
00108 *>
00109 *>    The "sizes" are specified by an array NN(1:NSIZES); the value of
00110 *>    each element NN(j) specifies one size.
00111 *>    The "types" are specified by a logical array DOTYPE( 1:NTYPES );
00112 *>    if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
00113 *>    Currently, the list of possible types is:
00114 *>
00115 *>    (1)  The zero matrix.
00116 *>    (2)  The identity matrix.
00117 *>    (3)  A (transposed) Jordan block, with 1's on the diagonal.
00118 *>
00119 *>    (4)  A diagonal matrix with evenly spaced entries
00120 *>         1, ..., ULP  and random complex angles.
00121 *>         (ULP = (first number larger than 1) - 1 )
00122 *>    (5)  A diagonal matrix with geometrically spaced entries
00123 *>         1, ..., ULP  and random complex angles.
00124 *>    (6)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
00125 *>         and random complex angles.
00126 *>
00127 *>    (7)  Same as (4), but multiplied by a constant near
00128 *>         the overflow threshold
00129 *>    (8)  Same as (4), but multiplied by a constant near
00130 *>         the underflow threshold
00131 *>
00132 *>    (9)  A matrix of the form  U' T U, where U is unitary and
00133 *>         T has evenly spaced entries 1, ..., ULP with random
00134 *>         complex angles on the diagonal and random O(1) entries in
00135 *>         the upper triangle.
00136 *>
00137 *>    (10) A matrix of the form  U' T U, where U is unitary and
00138 *>         T has geometrically spaced entries 1, ..., ULP with random
00139 *>         complex angles on the diagonal and random O(1) entries in
00140 *>         the upper triangle.
00141 *>
00142 *>    (11) A matrix of the form  U' T U, where U is orthogonal and
00143 *>         T has "clustered" entries 1, ULP,..., ULP with random
00144 *>         complex angles on the diagonal and random O(1) entries in
00145 *>         the upper triangle.
00146 *>
00147 *>    (12) A matrix of the form  U' T U, where U is unitary and
00148 *>         T has complex eigenvalues randomly chosen from
00149 *>         ULP < |z| < 1   and random O(1) entries in the upper
00150 *>         triangle.
00151 *>
00152 *>    (13) A matrix of the form  X' T X, where X has condition
00153 *>         SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
00154 *>         with random complex angles on the diagonal and random O(1)
00155 *>         entries in the upper triangle.
00156 *>
00157 *>    (14) A matrix of the form  X' T X, where X has condition
00158 *>         SQRT( ULP ) and T has geometrically spaced entries
00159 *>         1, ..., ULP with random complex angles on the diagonal
00160 *>         and random O(1) entries in the upper triangle.
00161 *>
00162 *>    (15) A matrix of the form  X' T X, where X has condition
00163 *>         SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
00164 *>         with random complex angles on the diagonal and random O(1)
00165 *>         entries in the upper triangle.
00166 *>
00167 *>    (16) A matrix of the form  X' T X, where X has condition
00168 *>         SQRT( ULP ) and T has complex eigenvalues randomly chosen
00169 *>         from ULP < |z| < 1 and random O(1) entries in the upper
00170 *>         triangle.
00171 *>
00172 *>    (17) Same as (16), but multiplied by a constant
00173 *>         near the overflow threshold
00174 *>    (18) Same as (16), but multiplied by a constant
00175 *>         near the underflow threshold
00176 *>
00177 *>    (19) Nonsymmetric matrix with random entries chosen from (-1,1).
00178 *>         If N is at least 4, all entries in first two rows and last
00179 *>         row, and first column and last two columns are zero.
00180 *>    (20) Same as (19), but multiplied by a constant
00181 *>         near the overflow threshold
00182 *>    (21) Same as (19), but multiplied by a constant
00183 *>         near the underflow threshold
00184 *>
00185 *>    In addition, an input file will be read from logical unit number
00186 *>    NIUNIT. The file contains matrices along with precomputed
00187 *>    eigenvalues and reciprocal condition numbers for the eigenvalue
00188 *>    average and right invariant subspace. For these matrices, in
00189 *>    addition to tests (1) to (15) we will compute the following two
00190 *>    tests:
00191 *>
00192 *>   (16)  |RCONDE - RCDEIN| / cond(RCONDE)
00193 *>
00194 *>      RCONDE is the reciprocal average eigenvalue condition number
00195 *>      computed by ZGEESX and RCDEIN (the precomputed true value)
00196 *>      is supplied as input.  cond(RCONDE) is the condition number
00197 *>      of RCONDE, and takes errors in computing RCONDE into account,
00198 *>      so that the resulting quantity should be O(ULP). cond(RCONDE)
00199 *>      is essentially given by norm(A)/RCONDV.
00200 *>
00201 *>   (17)  |RCONDV - RCDVIN| / cond(RCONDV)
00202 *>
00203 *>      RCONDV is the reciprocal right invariant subspace condition
00204 *>      number computed by ZGEESX and RCDVIN (the precomputed true
00205 *>      value) is supplied as input. cond(RCONDV) is the condition
00206 *>      number of RCONDV, and takes errors in computing RCONDV into
00207 *>      account, so that the resulting quantity should be O(ULP).
00208 *>      cond(RCONDV) is essentially given by norm(A)/RCONDE.
00209 *> \endverbatim
00210 *
00211 *  Arguments:
00212 *  ==========
00213 *
00214 *> \param[in] NSIZES
00215 *> \verbatim
00216 *>          NSIZES is INTEGER
00217 *>          The number of sizes of matrices to use.  NSIZES must be at
00218 *>          least zero. If it is zero, no randomly generated matrices
00219 *>          are tested, but any test matrices read from NIUNIT will be
00220 *>          tested.
00221 *> \endverbatim
00222 *>
00223 *> \param[in] NN
00224 *> \verbatim
00225 *>          NN is INTEGER array, dimension (NSIZES)
00226 *>          An array containing the sizes to be used for the matrices.
00227 *>          Zero values will be skipped.  The values must be at least
00228 *>          zero.
00229 *> \endverbatim
00230 *>
00231 *> \param[in] NTYPES
00232 *> \verbatim
00233 *>          NTYPES is INTEGER
00234 *>          The number of elements in DOTYPE. NTYPES must be at least
00235 *>          zero. If it is zero, no randomly generated test matrices
00236 *>          are tested, but and test matrices read from NIUNIT will be
00237 *>          tested. If it is MAXTYP+1 and NSIZES is 1, then an
00238 *>          additional type, MAXTYP+1 is defined, which is to use
00239 *>          whatever matrix is in A.  This is only useful if
00240 *>          DOTYPE(1:MAXTYP) is .FALSE. and DOTYPE(MAXTYP+1) is .TRUE. .
00241 *> \endverbatim
00242 *>
00243 *> \param[in] DOTYPE
00244 *> \verbatim
00245 *>          DOTYPE is LOGICAL array, dimension (NTYPES)
00246 *>          If DOTYPE(j) is .TRUE., then for each size in NN a
00247 *>          matrix of that size and of type j will be generated.
00248 *>          If NTYPES is smaller than the maximum number of types
00249 *>          defined (PARAMETER MAXTYP), then types NTYPES+1 through
00250 *>          MAXTYP will not be generated.  If NTYPES is larger
00251 *>          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
00252 *>          will be ignored.
00253 *> \endverbatim
00254 *>
00255 *> \param[in,out] ISEED
00256 *> \verbatim
00257 *>          ISEED is INTEGER array, dimension (4)
00258 *>          On entry ISEED specifies the seed of the random number
00259 *>          generator. The array elements should be between 0 and 4095;
00260 *>          if not they will be reduced mod 4096.  Also, ISEED(4) must
00261 *>          be odd.  The random number generator uses a linear
00262 *>          congruential sequence limited to small integers, and so
00263 *>          should produce machine independent random numbers. The
00264 *>          values of ISEED are changed on exit, and can be used in the
00265 *>          next call to ZDRVSX to continue the same random number
00266 *>          sequence.
00267 *> \endverbatim
00268 *>
00269 *> \param[in] THRESH
00270 *> \verbatim
00271 *>          THRESH is DOUBLE PRECISION
00272 *>          A test will count as "failed" if the "error", computed as
00273 *>          described above, exceeds THRESH.  Note that the error
00274 *>          is scaled to be O(1), so THRESH should be a reasonably
00275 *>          small multiple of 1, e.g., 10 or 100.  In particular,
00276 *>          it should not depend on the precision (single vs. double)
00277 *>          or the size of the matrix.  It must be at least zero.
00278 *> \endverbatim
00279 *>
00280 *> \param[in] NIUNIT
00281 *> \verbatim
00282 *>          NIUNIT is INTEGER
00283 *>          The FORTRAN unit number for reading in the data file of
00284 *>          problems to solve.
00285 *> \endverbatim
00286 *>
00287 *> \param[in] NOUNIT
00288 *> \verbatim
00289 *>          NOUNIT is INTEGER
00290 *>          The FORTRAN unit number for printing out error messages
00291 *>          (e.g., if a routine returns INFO not equal to 0.)
00292 *> \endverbatim
00293 *>
00294 *> \param[out] A
00295 *> \verbatim
00296 *>          A is COMPLEX*16 array, dimension (LDA, max(NN))
00297 *>          Used to hold the matrix whose eigenvalues are to be
00298 *>          computed.  On exit, A contains the last matrix actually used.
00299 *> \endverbatim
00300 *>
00301 *> \param[in] LDA
00302 *> \verbatim
00303 *>          LDA is INTEGER
00304 *>          The leading dimension of A, and H. LDA must be at
00305 *>          least 1 and at least max( NN ).
00306 *> \endverbatim
00307 *>
00308 *> \param[out] H
00309 *> \verbatim
00310 *>          H is COMPLEX*16 array, dimension (LDA, max(NN))
00311 *>          Another copy of the test matrix A, modified by ZGEESX.
00312 *> \endverbatim
00313 *>
00314 *> \param[out] HT
00315 *> \verbatim
00316 *>          HT is COMPLEX*16 array, dimension (LDA, max(NN))
00317 *>          Yet another copy of the test matrix A, modified by ZGEESX.
00318 *> \endverbatim
00319 *>
00320 *> \param[out] W
00321 *> \verbatim
00322 *>          W is COMPLEX*16 array, dimension (max(NN))
00323 *>          The computed eigenvalues of A.
00324 *> \endverbatim
00325 *>
00326 *> \param[out] WT
00327 *> \verbatim
00328 *>          WT is COMPLEX*16 array, dimension (max(NN))
00329 *>          Like W, this array contains the eigenvalues of A,
00330 *>          but those computed when ZGEESX only computes a partial
00331 *>          eigendecomposition, i.e. not Schur vectors
00332 *> \endverbatim
00333 *>
00334 *> \param[out] WTMP
00335 *> \verbatim
00336 *>          WTMP is COMPLEX*16 array, dimension (max(NN))
00337 *>          More temporary storage for eigenvalues.
00338 *> \endverbatim
00339 *>
00340 *> \param[out] VS
00341 *> \verbatim
00342 *>          VS is COMPLEX*16 array, dimension (LDVS, max(NN))
00343 *>          VS holds the computed Schur vectors.
00344 *> \endverbatim
00345 *>
00346 *> \param[in] LDVS
00347 *> \verbatim
00348 *>          LDVS is INTEGER
00349 *>          Leading dimension of VS. Must be at least max(1,max(NN)).
00350 *> \endverbatim
00351 *>
00352 *> \param[out] VS1
00353 *> \verbatim
00354 *>          VS1 is COMPLEX*16 array, dimension (LDVS, max(NN))
00355 *>          VS1 holds another copy of the computed Schur vectors.
00356 *> \endverbatim
00357 *>
00358 *> \param[out] RESULT
00359 *> \verbatim
00360 *>          RESULT is DOUBLE PRECISION array, dimension (17)
00361 *>          The values computed by the 17 tests described above.
00362 *>          The values are currently limited to 1/ulp, to avoid overflow.
00363 *> \endverbatim
00364 *>
00365 *> \param[out] WORK
00366 *> \verbatim
00367 *>          WORK is COMPLEX*16 array, dimension (LWORK)
00368 *> \endverbatim
00369 *>
00370 *> \param[in] LWORK
00371 *> \verbatim
00372 *>          LWORK is INTEGER
00373 *>          The number of entries in WORK.  This must be at least
00374 *>          max(1,2*NN(j)**2) for all j.
00375 *> \endverbatim
00376 *>
00377 *> \param[out] RWORK
00378 *> \verbatim
00379 *>          RWORK is DOUBLE PRECISION array, dimension (max(NN))
00380 *> \endverbatim
00381 *>
00382 *> \param[out] BWORK
00383 *> \verbatim
00384 *>          BWORK is LOGICAL array, dimension (max(NN))
00385 *> \endverbatim
00386 *>
00387 *> \param[out] INFO
00388 *> \verbatim
00389 *>          INFO is INTEGER
00390 *>          If 0,  successful exit.
00391 *>            <0,  input parameter -INFO is incorrect
00392 *>            >0,  ZLATMR, CLATMS, CLATME or ZGET24 returned an error
00393 *>                 code and INFO is its absolute value
00394 *>
00395 *>-----------------------------------------------------------------------
00396 *>
00397 *>     Some Local Variables and Parameters:
00398 *>     ---- ----- --------- --- ----------
00399 *>     ZERO, ONE       Real 0 and 1.
00400 *>     MAXTYP          The number of types defined.
00401 *>     NMAX            Largest value in NN.
00402 *>     NERRS           The number of tests which have exceeded THRESH
00403 *>     COND, CONDS,
00404 *>     IMODE           Values to be passed to the matrix generators.
00405 *>     ANORM           Norm of A; passed to matrix generators.
00406 *>
00407 *>     OVFL, UNFL      Overflow and underflow thresholds.
00408 *>     ULP, ULPINV     Finest relative precision and its inverse.
00409 *>     RTULP, RTULPI   Square roots of the previous 4 values.
00410 *>             The following four arrays decode JTYPE:
00411 *>     KTYPE(j)        The general type (1-10) for type "j".
00412 *>     KMODE(j)        The MODE value to be passed to the matrix
00413 *>                     generator for type "j".
00414 *>     KMAGN(j)        The order of magnitude ( O(1),
00415 *>                     O(overflow^(1/2) ), O(underflow^(1/2) )
00416 *>     KCONDS(j)       Selectw whether CONDS is to be 1 or
00417 *>                     1/sqrt(ulp).  (0 means irrelevant.)
00418 *> \endverbatim
00419 *
00420 *  Authors:
00421 *  ========
00422 *
00423 *> \author Univ. of Tennessee 
00424 *> \author Univ. of California Berkeley 
00425 *> \author Univ. of Colorado Denver 
00426 *> \author NAG Ltd. 
00427 *
00428 *> \date November 2011
00429 *
00430 *> \ingroup complex16_eig
00431 *
00432 *  =====================================================================
00433       SUBROUTINE ZDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00434      $                   NIUNIT, NOUNIT, A, LDA, H, HT, W, WT, WTMP, VS,
00435      $                   LDVS, VS1, RESULT, WORK, LWORK, RWORK, BWORK,
00436      $                   INFO )
00437 *
00438 *  -- LAPACK test routine (version 3.4.0) --
00439 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00440 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00441 *     November 2011
00442 *
00443 *     .. Scalar Arguments ..
00444       INTEGER            INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
00445      $                   NTYPES
00446       DOUBLE PRECISION   THRESH
00447 *     ..
00448 *     .. Array Arguments ..
00449       LOGICAL            BWORK( * ), DOTYPE( * )
00450       INTEGER            ISEED( 4 ), NN( * )
00451       DOUBLE PRECISION   RESULT( 17 ), RWORK( * )
00452       COMPLEX*16         A( LDA, * ), H( LDA, * ), HT( LDA, * ),
00453      $                   VS( LDVS, * ), VS1( LDVS, * ), W( * ),
00454      $                   WORK( * ), WT( * ), WTMP( * )
00455 *     ..
00456 *
00457 *  =====================================================================
00458 *
00459 *     .. Parameters ..
00460       COMPLEX*16         CZERO
00461       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ) )
00462       COMPLEX*16         CONE
00463       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
00464       DOUBLE PRECISION   ZERO, ONE
00465       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00466       INTEGER            MAXTYP
00467       PARAMETER          ( MAXTYP = 21 )
00468 *     ..
00469 *     .. Local Scalars ..
00470       LOGICAL            BADNN
00471       CHARACTER*3        PATH
00472       INTEGER            I, IINFO, IMODE, ISRT, ITYPE, IWK, J, JCOL,
00473      $                   JSIZE, JTYPE, MTYPES, N, NERRS, NFAIL, NMAX,
00474      $                   NNWORK, NSLCT, NTEST, NTESTF, NTESTT
00475       DOUBLE PRECISION   ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN,
00476      $                   RTULP, RTULPI, ULP, ULPINV, UNFL
00477 *     ..
00478 *     .. Local Arrays ..
00479       INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ),
00480      $                   KCONDS( MAXTYP ), KMAGN( MAXTYP ),
00481      $                   KMODE( MAXTYP ), KTYPE( MAXTYP )
00482 *     ..
00483 *     .. Arrays in Common ..
00484       LOGICAL            SELVAL( 20 )
00485       DOUBLE PRECISION   SELWI( 20 ), SELWR( 20 )
00486 *     ..
00487 *     .. Scalars in Common ..
00488       INTEGER            SELDIM, SELOPT
00489 *     ..
00490 *     .. Common blocks ..
00491       COMMON             / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
00492 *     ..
00493 *     .. External Functions ..
00494       DOUBLE PRECISION   DLAMCH
00495       EXTERNAL           DLAMCH
00496 *     ..
00497 *     .. External Subroutines ..
00498       EXTERNAL           DLABAD, DLASUM, XERBLA, ZGET24, ZLASET, ZLATME,
00499      $                   ZLATMR, ZLATMS
00500 *     ..
00501 *     .. Intrinsic Functions ..
00502       INTRINSIC          ABS, MAX, MIN, SQRT
00503 *     ..
00504 *     .. Data statements ..
00505       DATA               KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
00506       DATA               KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
00507      $                   3, 1, 2, 3 /
00508       DATA               KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
00509      $                   1, 5, 5, 5, 4, 3, 1 /
00510       DATA               KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
00511 *     ..
00512 *     .. Executable Statements ..
00513 *
00514       PATH( 1: 1 ) = 'Zomplex precision'
00515       PATH( 2: 3 ) = 'SX'
00516 *
00517 *     Check for errors
00518 *
00519       NTESTT = 0
00520       NTESTF = 0
00521       INFO = 0
00522 *
00523 *     Important constants
00524 *
00525       BADNN = .FALSE.
00526 *
00527 *     8 is the largest dimension in the input file of precomputed
00528 *     problems
00529 *
00530       NMAX = 8
00531       DO 10 J = 1, NSIZES
00532          NMAX = MAX( NMAX, NN( J ) )
00533          IF( NN( J ).LT.0 )
00534      $      BADNN = .TRUE.
00535    10 CONTINUE
00536 *
00537 *     Check for errors
00538 *
00539       IF( NSIZES.LT.0 ) THEN
00540          INFO = -1
00541       ELSE IF( BADNN ) THEN
00542          INFO = -2
00543       ELSE IF( NTYPES.LT.0 ) THEN
00544          INFO = -3
00545       ELSE IF( THRESH.LT.ZERO ) THEN
00546          INFO = -6
00547       ELSE IF( NIUNIT.LE.0 ) THEN
00548          INFO = -7
00549       ELSE IF( NOUNIT.LE.0 ) THEN
00550          INFO = -8
00551       ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
00552          INFO = -10
00553       ELSE IF( LDVS.LT.1 .OR. LDVS.LT.NMAX ) THEN
00554          INFO = -20
00555       ELSE IF( MAX( 3*NMAX, 2*NMAX**2 ).GT.LWORK ) THEN
00556          INFO = -24
00557       END IF
00558 *
00559       IF( INFO.NE.0 ) THEN
00560          CALL XERBLA( 'ZDRVSX', -INFO )
00561          RETURN
00562       END IF
00563 *
00564 *     If nothing to do check on NIUNIT
00565 *
00566       IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
00567      $   GO TO 150
00568 *
00569 *     More Important constants
00570 *
00571       UNFL = DLAMCH( 'Safe minimum' )
00572       OVFL = ONE / UNFL
00573       CALL DLABAD( UNFL, OVFL )
00574       ULP = DLAMCH( 'Precision' )
00575       ULPINV = ONE / ULP
00576       RTULP = SQRT( ULP )
00577       RTULPI = ONE / RTULP
00578 *
00579 *     Loop over sizes, types
00580 *
00581       NERRS = 0
00582 *
00583       DO 140 JSIZE = 1, NSIZES
00584          N = NN( JSIZE )
00585          IF( NSIZES.NE.1 ) THEN
00586             MTYPES = MIN( MAXTYP, NTYPES )
00587          ELSE
00588             MTYPES = MIN( MAXTYP+1, NTYPES )
00589          END IF
00590 *
00591          DO 130 JTYPE = 1, MTYPES
00592             IF( .NOT.DOTYPE( JTYPE ) )
00593      $         GO TO 130
00594 *
00595 *           Save ISEED in case of an error.
00596 *
00597             DO 20 J = 1, 4
00598                IOLDSD( J ) = ISEED( J )
00599    20       CONTINUE
00600 *
00601 *           Compute "A"
00602 *
00603 *           Control parameters:
00604 *
00605 *           KMAGN  KCONDS  KMODE        KTYPE
00606 *       =1  O(1)   1       clustered 1  zero
00607 *       =2  large  large   clustered 2  identity
00608 *       =3  small          exponential  Jordan
00609 *       =4                 arithmetic   diagonal, (w/ eigenvalues)
00610 *       =5                 random log   symmetric, w/ eigenvalues
00611 *       =6                 random       general, w/ eigenvalues
00612 *       =7                              random diagonal
00613 *       =8                              random symmetric
00614 *       =9                              random general
00615 *       =10                             random triangular
00616 *
00617             IF( MTYPES.GT.MAXTYP )
00618      $         GO TO 90
00619 *
00620             ITYPE = KTYPE( JTYPE )
00621             IMODE = KMODE( JTYPE )
00622 *
00623 *           Compute norm
00624 *
00625             GO TO ( 30, 40, 50 )KMAGN( JTYPE )
00626 *
00627    30       CONTINUE
00628             ANORM = ONE
00629             GO TO 60
00630 *
00631    40       CONTINUE
00632             ANORM = OVFL*ULP
00633             GO TO 60
00634 *
00635    50       CONTINUE
00636             ANORM = UNFL*ULPINV
00637             GO TO 60
00638 *
00639    60       CONTINUE
00640 *
00641             CALL ZLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
00642             IINFO = 0
00643             COND = ULPINV
00644 *
00645 *           Special Matrices -- Identity & Jordan block
00646 *
00647             IF( ITYPE.EQ.1 ) THEN
00648 *
00649 *              Zero
00650 *
00651                IINFO = 0
00652 *
00653             ELSE IF( ITYPE.EQ.2 ) THEN
00654 *
00655 *              Identity
00656 *
00657                DO 70 JCOL = 1, N
00658                   A( JCOL, JCOL ) = ANORM
00659    70          CONTINUE
00660 *
00661             ELSE IF( ITYPE.EQ.3 ) THEN
00662 *
00663 *              Jordan Block
00664 *
00665                DO 80 JCOL = 1, N
00666                   A( JCOL, JCOL ) = ANORM
00667                   IF( JCOL.GT.1 )
00668      $               A( JCOL, JCOL-1 ) = CONE
00669    80          CONTINUE
00670 *
00671             ELSE IF( ITYPE.EQ.4 ) THEN
00672 *
00673 *              Diagonal Matrix, [Eigen]values Specified
00674 *
00675                CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
00676      $                      ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
00677      $                      IINFO )
00678 *
00679             ELSE IF( ITYPE.EQ.5 ) THEN
00680 *
00681 *              Symmetric, eigenvalues specified
00682 *
00683                CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
00684      $                      ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
00685      $                      IINFO )
00686 *
00687             ELSE IF( ITYPE.EQ.6 ) THEN
00688 *
00689 *              General, eigenvalues specified
00690 *
00691                IF( KCONDS( JTYPE ).EQ.1 ) THEN
00692                   CONDS = ONE
00693                ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
00694                   CONDS = RTULPI
00695                ELSE
00696                   CONDS = ZERO
00697                END IF
00698 *
00699                CALL ZLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE,
00700      $                      'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
00701      $                      A, LDA, WORK( 2*N+1 ), IINFO )
00702 *
00703             ELSE IF( ITYPE.EQ.7 ) THEN
00704 *
00705 *              Diagonal, random eigenvalues
00706 *
00707                CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
00708      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00709      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
00710      $                      ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
00711 *
00712             ELSE IF( ITYPE.EQ.8 ) THEN
00713 *
00714 *              Symmetric, random eigenvalues
00715 *
00716                CALL ZLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
00717      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00718      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
00719      $                      ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
00720 *
00721             ELSE IF( ITYPE.EQ.9 ) THEN
00722 *
00723 *              General, random eigenvalues
00724 *
00725                CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
00726      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00727      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
00728      $                      ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
00729                IF( N.GE.4 ) THEN
00730                   CALL ZLASET( 'Full', 2, N, CZERO, CZERO, A, LDA )
00731                   CALL ZLASET( 'Full', N-3, 1, CZERO, CZERO, A( 3, 1 ),
00732      $                         LDA )
00733                   CALL ZLASET( 'Full', N-3, 2, CZERO, CZERO,
00734      $                         A( 3, N-1 ), LDA )
00735                   CALL ZLASET( 'Full', 1, N, CZERO, CZERO, A( N, 1 ),
00736      $                         LDA )
00737                END IF
00738 *
00739             ELSE IF( ITYPE.EQ.10 ) THEN
00740 *
00741 *              Triangular, random eigenvalues
00742 *
00743                CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
00744      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00745      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
00746      $                      ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
00747 *
00748             ELSE
00749 *
00750                IINFO = 1
00751             END IF
00752 *
00753             IF( IINFO.NE.0 ) THEN
00754                WRITE( NOUNIT, FMT = 9991 )'Generator', IINFO, N, JTYPE,
00755      $            IOLDSD
00756                INFO = ABS( IINFO )
00757                RETURN
00758             END IF
00759 *
00760    90       CONTINUE
00761 *
00762 *           Test for minimal and generous workspace
00763 *
00764             DO 120 IWK = 1, 2
00765                IF( IWK.EQ.1 ) THEN
00766                   NNWORK = 2*N
00767                ELSE
00768                   NNWORK = MAX( 2*N, N*( N+1 ) / 2 )
00769                END IF
00770                NNWORK = MAX( NNWORK, 1 )
00771 *
00772                CALL ZGET24( .FALSE., JTYPE, THRESH, IOLDSD, NOUNIT, N,
00773      $                      A, LDA, H, HT, W, WT, WTMP, VS, LDVS, VS1,
00774      $                      RCDEIN, RCDVIN, NSLCT, ISLCT, 0, RESULT,
00775      $                      WORK, NNWORK, RWORK, BWORK, INFO )
00776 *
00777 *              Check for RESULT(j) > THRESH
00778 *
00779                NTEST = 0
00780                NFAIL = 0
00781                DO 100 J = 1, 15
00782                   IF( RESULT( J ).GE.ZERO )
00783      $               NTEST = NTEST + 1
00784                   IF( RESULT( J ).GE.THRESH )
00785      $               NFAIL = NFAIL + 1
00786   100          CONTINUE
00787 *
00788                IF( NFAIL.GT.0 )
00789      $            NTESTF = NTESTF + 1
00790                IF( NTESTF.EQ.1 ) THEN
00791                   WRITE( NOUNIT, FMT = 9999 )PATH
00792                   WRITE( NOUNIT, FMT = 9998 )
00793                   WRITE( NOUNIT, FMT = 9997 )
00794                   WRITE( NOUNIT, FMT = 9996 )
00795                   WRITE( NOUNIT, FMT = 9995 )THRESH
00796                   WRITE( NOUNIT, FMT = 9994 )
00797                   NTESTF = 2
00798                END IF
00799 *
00800                DO 110 J = 1, 15
00801                   IF( RESULT( J ).GE.THRESH ) THEN
00802                      WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE,
00803      $                  J, RESULT( J )
00804                   END IF
00805   110          CONTINUE
00806 *
00807                NERRS = NERRS + NFAIL
00808                NTESTT = NTESTT + NTEST
00809 *
00810   120       CONTINUE
00811   130    CONTINUE
00812   140 CONTINUE
00813 *
00814   150 CONTINUE
00815 *
00816 *     Read in data from file to check accuracy of condition estimation
00817 *     Read input data until N=0
00818 *
00819       JTYPE = 0
00820   160 CONTINUE
00821       READ( NIUNIT, FMT = *, END = 200 )N, NSLCT, ISRT
00822       IF( N.EQ.0 )
00823      $   GO TO 200
00824       JTYPE = JTYPE + 1
00825       ISEED( 1 ) = JTYPE
00826       READ( NIUNIT, FMT = * )( ISLCT( I ), I = 1, NSLCT )
00827       DO 170 I = 1, N
00828          READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N )
00829   170 CONTINUE
00830       READ( NIUNIT, FMT = * )RCDEIN, RCDVIN
00831 *
00832       CALL ZGET24( .TRUE., 22, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT,
00833      $             W, WT, WTMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT,
00834      $             ISLCT, ISRT, RESULT, WORK, LWORK, RWORK, BWORK,
00835      $             INFO )
00836 *
00837 *     Check for RESULT(j) > THRESH
00838 *
00839       NTEST = 0
00840       NFAIL = 0
00841       DO 180 J = 1, 17
00842          IF( RESULT( J ).GE.ZERO )
00843      $      NTEST = NTEST + 1
00844          IF( RESULT( J ).GE.THRESH )
00845      $      NFAIL = NFAIL + 1
00846   180 CONTINUE
00847 *
00848       IF( NFAIL.GT.0 )
00849      $   NTESTF = NTESTF + 1
00850       IF( NTESTF.EQ.1 ) THEN
00851          WRITE( NOUNIT, FMT = 9999 )PATH
00852          WRITE( NOUNIT, FMT = 9998 )
00853          WRITE( NOUNIT, FMT = 9997 )
00854          WRITE( NOUNIT, FMT = 9996 )
00855          WRITE( NOUNIT, FMT = 9995 )THRESH
00856          WRITE( NOUNIT, FMT = 9994 )
00857          NTESTF = 2
00858       END IF
00859       DO 190 J = 1, 17
00860          IF( RESULT( J ).GE.THRESH ) THEN
00861             WRITE( NOUNIT, FMT = 9992 )N, JTYPE, J, RESULT( J )
00862          END IF
00863   190 CONTINUE
00864 *
00865       NERRS = NERRS + NFAIL
00866       NTESTT = NTESTT + NTEST
00867       GO TO 160
00868   200 CONTINUE
00869 *
00870 *     Summary
00871 *
00872       CALL DLASUM( PATH, NOUNIT, NERRS, NTESTT )
00873 *
00874  9999 FORMAT( / 1X, A3, ' -- Complex Schur Form Decomposition Expert ',
00875      $      'Driver', / ' Matrix types (see ZDRVSX for details): ' )
00876 *
00877  9998 FORMAT( / ' Special Matrices:', / '  1=Zero matrix.             ',
00878      $      '           ', '  5=Diagonal: geometr. spaced entries.',
00879      $      / '  2=Identity matrix.                    ', '  6=Diagona',
00880      $      'l: clustered entries.', / '  3=Transposed Jordan block.  ',
00881      $      '          ', '  7=Diagonal: large, evenly spaced.', / '  ',
00882      $      '4=Diagonal: evenly spaced entries.    ', '  8=Diagonal: s',
00883      $      'mall, evenly spaced.' )
00884  9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / '  9=Well-cond., ev',
00885      $      'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
00886      $      'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
00887      $      ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
00888      $      'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
00889      $      'lex ', / ' 12=Well-cond., random complex ', '         ',
00890      $      ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
00891      $      'tioned, evenly spaced.     ', ' 18=Ill-cond., small rand.',
00892      $      ' complx ' )
00893  9996 FORMAT( ' 19=Matrix with random O(1) entries.    ', ' 21=Matrix ',
00894      $      'with small random entries.', / ' 20=Matrix with large ran',
00895      $      'dom entries.   ', / )
00896  9995 FORMAT( ' Tests performed with test threshold =', F8.2,
00897      $      / ' ( A denotes A on input and T denotes A on output)',
00898      $      / / ' 1 = 0 if T in Schur form (no sort), ',
00899      $      '  1/ulp otherwise', /
00900      $      ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
00901      $      / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ',
00902      $      / ' 4 = 0 if W are eigenvalues of T (no sort),',
00903      $      '  1/ulp otherwise', /
00904      $      ' 5 = 0 if T same no matter if VS computed (no sort),',
00905      $      '  1/ulp otherwise', /
00906      $      ' 6 = 0 if W same no matter if VS computed (no sort)',
00907      $      ',  1/ulp otherwise' )
00908  9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', '  1/ulp otherwise',
00909      $      / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
00910      $      / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
00911      $      / ' 10 = 0 if W are eigenvalues of T (sort),',
00912      $      '  1/ulp otherwise', /
00913      $      ' 11 = 0 if T same no matter what else computed (sort),',
00914      $      '  1/ulp otherwise', /
00915      $      ' 12 = 0 if W same no matter what else computed ',
00916      $      '(sort), 1/ulp otherwise', /
00917      $      ' 13 = 0 if sorting succesful, 1/ulp otherwise',
00918      $      / ' 14 = 0 if RCONDE same no matter what else computed,',
00919      $      ' 1/ulp otherwise', /
00920      $      ' 15 = 0 if RCONDv same no matter what else computed,',
00921      $      ' 1/ulp otherwise', /
00922      $      ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
00923      $      / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
00924  9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ),
00925      $      ' type ', I2, ', test(', I2, ')=', G10.3 )
00926  9992 FORMAT( ' N=', I5, ', input example =', I3, ',  test(', I2, ')=',
00927      $      G10.3 )
00928  9991 FORMAT( ' ZDRVSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00929      $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
00930 *
00931       RETURN
00932 *
00933 *     End of ZDRVSX
00934 *
00935       END
 All Files Functions