LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sdrvsx.f
Go to the documentation of this file.
00001 *> \brief \b SDRVSX
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 SDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00012 *                          NIUNIT, NOUNIT, A, LDA, H, HT, WR, WI, WRT,
00013 *                          WIT, WRTMP, WITMP, VS, LDVS, VS1, RESULT, WORK,
00014 *                          LWORK, IWORK, BWORK, INFO )
00015 * 
00016 *       .. Scalar Arguments ..
00017 *       INTEGER            INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
00018 *      $                   NTYPES
00019 *       REAL               THRESH
00020 *       ..
00021 *       .. Array Arguments ..
00022 *       LOGICAL            BWORK( * ), DOTYPE( * )
00023 *       INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
00024 *       REAL               A( LDA, * ), H( LDA, * ), HT( LDA, * ),
00025 *      $                   RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
00026 *      $                   WI( * ), WIT( * ), WITMP( * ), WORK( * ),
00027 *      $                   WR( * ), WRT( * ), WRTMP( * )
00028 *       ..
00029 *  
00030 *
00031 *> \par Purpose:
00032 *  =============
00033 *>
00034 *> \verbatim
00035 *>
00036 *>    SDRVSX checks the nonsymmetric eigenvalue (Schur form) problem
00037 *>    expert driver SGEESX.
00038 *>
00039 *>    SDRVSX 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 SDRVSX 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 WR+sqrt(-1)*WI 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 WR+sqrt(-1)*WI are eigenvalues of T
00083 *>            1/ulp otherwise
00084 *>            If workspace sufficient, also compare WR, WI 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 signs.
00121 *>         (ULP = (first number larger than 1) - 1 )
00122 *>    (5)  A diagonal matrix with geometrically spaced entries
00123 *>         1, ..., ULP  and random signs.
00124 *>    (6)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
00125 *>         and random signs.
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 orthogonal and
00133 *>         T has evenly spaced entries 1, ..., ULP with random signs
00134 *>         on the diagonal and random O(1) entries in the upper
00135 *>         triangle.
00136 *>
00137 *>    (10) A matrix of the form  U' T U, where U is orthogonal and
00138 *>         T has geometrically spaced entries 1, ..., ULP with random
00139 *>         signs on the diagonal and random O(1) entries in the upper
00140 *>         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 *>         signs on the diagonal and random O(1) entries in the upper
00145 *>         triangle.
00146 *>
00147 *>    (12) A matrix of the form  U' T U, where U is orthogonal and
00148 *>         T has real or complex conjugate paired eigenvalues randomly
00149 *>         chosen from ( ULP, 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 signs on the diagonal and random O(1) entries
00155 *>         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 signs on the diagonal and random
00160 *>         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 signs on the diagonal and random O(1) entries
00165 *>         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 real or complex conjugate paired
00169 *>         eigenvalues randomly chosen from ( ULP, 1 ) and random
00170 *>         O(1) entries in the upper 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 SGEESX 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 SGEESX 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 SDRVSX to continue the same random number
00266 *>          sequence.
00267 *> \endverbatim
00268 *>
00269 *> \param[in] THRESH
00270 *> \verbatim
00271 *>          THRESH is REAL
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 REAL 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 REAL array, dimension (LDA, max(NN))
00311 *>          Another copy of the test matrix A, modified by SGEESX.
00312 *> \endverbatim
00313 *>
00314 *> \param[out] HT
00315 *> \verbatim
00316 *>          HT is REAL array, dimension (LDA, max(NN))
00317 *>          Yet another copy of the test matrix A, modified by SGEESX.
00318 *> \endverbatim
00319 *>
00320 *> \param[out] WR
00321 *> \verbatim
00322 *>          WR is REAL array, dimension (max(NN))
00323 *> \endverbatim
00324 *>
00325 *> \param[out] WI
00326 *> \verbatim
00327 *>          WI is REAL array, dimension (max(NN))
00328 *>
00329 *>          The real and imaginary parts of the eigenvalues of A.
00330 *>          On exit, WR + WI*i are the eigenvalues of the matrix in A.
00331 *> \endverbatim
00332 *>
00333 *> \param[out] WRT
00334 *> \verbatim
00335 *>          WRT is REAL array, dimension (max(NN))
00336 *> \endverbatim
00337 *>
00338 *> \param[out] WIT
00339 *> \verbatim
00340 *>          WIT is REAL array, dimension (max(NN))
00341 *>
00342 *>          Like WR, WI, these arrays contain the eigenvalues of A,
00343 *>          but those computed when SGEESX only computes a partial
00344 *>          eigendecomposition, i.e. not Schur vectors
00345 *> \endverbatim
00346 *>
00347 *> \param[out] WRTMP
00348 *> \verbatim
00349 *>          WRTMP is REAL array, dimension (max(NN))
00350 *> \endverbatim
00351 *>
00352 *> \param[out] WITMP
00353 *> \verbatim
00354 *>          WITMP is REAL array, dimension (max(NN))
00355 *>
00356 *>          More temporary storage for eigenvalues.
00357 *> \endverbatim
00358 *>
00359 *> \param[out] VS
00360 *> \verbatim
00361 *>          VS is REAL array, dimension (LDVS, max(NN))
00362 *>          VS holds the computed Schur vectors.
00363 *> \endverbatim
00364 *>
00365 *> \param[in] LDVS
00366 *> \verbatim
00367 *>          LDVS is INTEGER
00368 *>          Leading dimension of VS. Must be at least max(1,max(NN)).
00369 *> \endverbatim
00370 *>
00371 *> \param[out] VS1
00372 *> \verbatim
00373 *>          VS1 is REAL array, dimension (LDVS, max(NN))
00374 *>          VS1 holds another copy of the computed Schur vectors.
00375 *> \endverbatim
00376 *>
00377 *> \param[out] RESULT
00378 *> \verbatim
00379 *>          RESULT is REAL array, dimension (17)
00380 *>          The values computed by the 17 tests described above.
00381 *>          The values are currently limited to 1/ulp, to avoid overflow.
00382 *> \endverbatim
00383 *>
00384 *> \param[out] WORK
00385 *> \verbatim
00386 *>          WORK is REAL array, dimension (LWORK)
00387 *> \endverbatim
00388 *>
00389 *> \param[in] LWORK
00390 *> \verbatim
00391 *>          LWORK is INTEGER
00392 *>          The number of entries in WORK.  This must be at least
00393 *>          max(3*NN(j),2*NN(j)**2) for all j.
00394 *> \endverbatim
00395 *>
00396 *> \param[out] IWORK
00397 *> \verbatim
00398 *>          IWORK is INTEGER array, dimension (max(NN)*max(NN))
00399 *> \endverbatim
00400 *>
00401 *> \param[out] BWORK
00402 *> \verbatim
00403 *>          BWORK is LOGICAL array, dimension (max(NN))
00404 *> \endverbatim
00405 *>
00406 *> \param[out] INFO
00407 *> \verbatim
00408 *>          INFO is INTEGER
00409 *>          If 0,  successful exit.
00410 *>            <0,  input parameter -INFO is incorrect
00411 *>            >0,  SLATMR, SLATMS, SLATME or SGET24 returned an error
00412 *>                 code and INFO is its absolute value
00413 *>
00414 *>-----------------------------------------------------------------------
00415 *>
00416 *>     Some Local Variables and Parameters:
00417 *>     ---- ----- --------- --- ----------
00418 *>     ZERO, ONE       Real 0 and 1.
00419 *>     MAXTYP          The number of types defined.
00420 *>     NMAX            Largest value in NN.
00421 *>     NERRS           The number of tests which have exceeded THRESH
00422 *>     COND, CONDS,
00423 *>     IMODE           Values to be passed to the matrix generators.
00424 *>     ANORM           Norm of A; passed to matrix generators.
00425 *>
00426 *>     OVFL, UNFL      Overflow and underflow thresholds.
00427 *>     ULP, ULPINV     Finest relative precision and its inverse.
00428 *>     RTULP, RTULPI   Square roots of the previous 4 values.
00429 *>             The following four arrays decode JTYPE:
00430 *>     KTYPE(j)        The general type (1-10) for type "j".
00431 *>     KMODE(j)        The MODE value to be passed to the matrix
00432 *>                     generator for type "j".
00433 *>     KMAGN(j)        The order of magnitude ( O(1),
00434 *>                     O(overflow^(1/2) ), O(underflow^(1/2) )
00435 *>     KCONDS(j)       Selectw whether CONDS is to be 1 or
00436 *>                     1/sqrt(ulp).  (0 means irrelevant.)
00437 *> \endverbatim
00438 *
00439 *  Authors:
00440 *  ========
00441 *
00442 *> \author Univ. of Tennessee 
00443 *> \author Univ. of California Berkeley 
00444 *> \author Univ. of Colorado Denver 
00445 *> \author NAG Ltd. 
00446 *
00447 *> \date November 2011
00448 *
00449 *> \ingroup single_eig
00450 *
00451 *  =====================================================================
00452       SUBROUTINE SDRVSX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00453      $                   NIUNIT, NOUNIT, A, LDA, H, HT, WR, WI, WRT,
00454      $                   WIT, WRTMP, WITMP, VS, LDVS, VS1, RESULT, WORK,
00455      $                   LWORK, IWORK, BWORK, INFO )
00456 *
00457 *  -- LAPACK test routine (version 3.4.0) --
00458 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00459 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00460 *     November 2011
00461 *
00462 *     .. Scalar Arguments ..
00463       INTEGER            INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
00464      $                   NTYPES
00465       REAL               THRESH
00466 *     ..
00467 *     .. Array Arguments ..
00468       LOGICAL            BWORK( * ), DOTYPE( * )
00469       INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
00470       REAL               A( LDA, * ), H( LDA, * ), HT( LDA, * ),
00471      $                   RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
00472      $                   WI( * ), WIT( * ), WITMP( * ), WORK( * ),
00473      $                   WR( * ), WRT( * ), WRTMP( * )
00474 *     ..
00475 *
00476 *  =====================================================================
00477 *
00478 *     .. Parameters ..
00479       REAL               ZERO, ONE
00480       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00481       INTEGER            MAXTYP
00482       PARAMETER          ( MAXTYP = 21 )
00483 *     ..
00484 *     .. Local Scalars ..
00485       LOGICAL            BADNN
00486       CHARACTER*3        PATH
00487       INTEGER            I, IINFO, IMODE, ITYPE, IWK, J, JCOL, JSIZE,
00488      $                   JTYPE, MTYPES, N, NERRS, NFAIL, NMAX,
00489      $                   NNWORK, NSLCT, NTEST, NTESTF, NTESTT
00490       REAL               ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN,
00491      $                   RTULP, RTULPI, ULP, ULPINV, UNFL
00492 *     ..
00493 *     .. Local Arrays ..
00494       CHARACTER          ADUMMA( 1 )
00495       INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ),
00496      $                   KCONDS( MAXTYP ), KMAGN( MAXTYP ),
00497      $                   KMODE( MAXTYP ), KTYPE( MAXTYP )
00498 *     ..
00499 *     .. Arrays in Common ..
00500       LOGICAL            SELVAL( 20 )
00501       REAL               SELWI( 20 ), SELWR( 20 )
00502 *     ..
00503 *     .. Scalars in Common ..
00504       INTEGER            SELDIM, SELOPT
00505 *     ..
00506 *     .. Common blocks ..
00507       COMMON             / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
00508 *     ..
00509 *     .. External Functions ..
00510       REAL               SLAMCH
00511       EXTERNAL           SLAMCH
00512 *     ..
00513 *     .. External Subroutines ..
00514       EXTERNAL           SGET24, SLABAD, SLASUM, SLATME, SLATMR, SLATMS,
00515      $                   SLASET, XERBLA
00516 *     ..
00517 *     .. Intrinsic Functions ..
00518       INTRINSIC          ABS, MAX, MIN, SQRT
00519 *     ..
00520 *     .. Data statements ..
00521       DATA               KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
00522       DATA               KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
00523      $                   3, 1, 2, 3 /
00524       DATA               KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
00525      $                   1, 5, 5, 5, 4, 3, 1 /
00526       DATA               KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
00527 *     ..
00528 *     .. Executable Statements ..
00529 *
00530       PATH( 1: 1 ) = 'Single precision'
00531       PATH( 2: 3 ) = 'SX'
00532 *
00533 *     Check for errors
00534 *
00535       NTESTT = 0
00536       NTESTF = 0
00537       INFO = 0
00538 *
00539 *     Important constants
00540 *
00541       BADNN = .FALSE.
00542 *
00543 *     12 is the largest dimension in the input file of precomputed
00544 *     problems
00545 *
00546       NMAX = 12
00547       DO 10 J = 1, NSIZES
00548          NMAX = MAX( NMAX, NN( J ) )
00549          IF( NN( J ).LT.0 )
00550      $      BADNN = .TRUE.
00551    10 CONTINUE
00552 *
00553 *     Check for errors
00554 *
00555       IF( NSIZES.LT.0 ) THEN
00556          INFO = -1
00557       ELSE IF( BADNN ) THEN
00558          INFO = -2
00559       ELSE IF( NTYPES.LT.0 ) THEN
00560          INFO = -3
00561       ELSE IF( THRESH.LT.ZERO ) THEN
00562          INFO = -6
00563       ELSE IF( NIUNIT.LE.0 ) THEN
00564          INFO = -7
00565       ELSE IF( NOUNIT.LE.0 ) THEN
00566          INFO = -8
00567       ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
00568          INFO = -10
00569       ELSE IF( LDVS.LT.1 .OR. LDVS.LT.NMAX ) THEN
00570          INFO = -20
00571       ELSE IF( MAX( 3*NMAX, 2*NMAX**2 ).GT.LWORK ) THEN
00572          INFO = -24
00573       END IF
00574 *
00575       IF( INFO.NE.0 ) THEN
00576          CALL XERBLA( 'SDRVSX', -INFO )
00577          RETURN
00578       END IF
00579 *
00580 *     If nothing to do check on NIUNIT
00581 *
00582       IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
00583      $   GO TO 150
00584 *
00585 *     More Important constants
00586 *
00587       UNFL = SLAMCH( 'Safe minimum' )
00588       OVFL = ONE / UNFL
00589       CALL SLABAD( UNFL, OVFL )
00590       ULP = SLAMCH( 'Precision' )
00591       ULPINV = ONE / ULP
00592       RTULP = SQRT( ULP )
00593       RTULPI = ONE / RTULP
00594 *
00595 *     Loop over sizes, types
00596 *
00597       NERRS = 0
00598 *
00599       DO 140 JSIZE = 1, NSIZES
00600          N = NN( JSIZE )
00601          IF( NSIZES.NE.1 ) THEN
00602             MTYPES = MIN( MAXTYP, NTYPES )
00603          ELSE
00604             MTYPES = MIN( MAXTYP+1, NTYPES )
00605          END IF
00606 *
00607          DO 130 JTYPE = 1, MTYPES
00608             IF( .NOT.DOTYPE( JTYPE ) )
00609      $         GO TO 130
00610 *
00611 *           Save ISEED in case of an error.
00612 *
00613             DO 20 J = 1, 4
00614                IOLDSD( J ) = ISEED( J )
00615    20       CONTINUE
00616 *
00617 *           Compute "A"
00618 *
00619 *           Control parameters:
00620 *
00621 *           KMAGN  KCONDS  KMODE        KTYPE
00622 *       =1  O(1)   1       clustered 1  zero
00623 *       =2  large  large   clustered 2  identity
00624 *       =3  small          exponential  Jordan
00625 *       =4                 arithmetic   diagonal, (w/ eigenvalues)
00626 *       =5                 random log   symmetric, w/ eigenvalues
00627 *       =6                 random       general, w/ eigenvalues
00628 *       =7                              random diagonal
00629 *       =8                              random symmetric
00630 *       =9                              random general
00631 *       =10                             random triangular
00632 *
00633             IF( MTYPES.GT.MAXTYP )
00634      $         GO TO 90
00635 *
00636             ITYPE = KTYPE( JTYPE )
00637             IMODE = KMODE( JTYPE )
00638 *
00639 *           Compute norm
00640 *
00641             GO TO ( 30, 40, 50 )KMAGN( JTYPE )
00642 *
00643    30       CONTINUE
00644             ANORM = ONE
00645             GO TO 60
00646 *
00647    40       CONTINUE
00648             ANORM = OVFL*ULP
00649             GO TO 60
00650 *
00651    50       CONTINUE
00652             ANORM = UNFL*ULPINV
00653             GO TO 60
00654 *
00655    60       CONTINUE
00656 *
00657             CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
00658             IINFO = 0
00659             COND = ULPINV
00660 *
00661 *           Special Matrices -- Identity & Jordan block
00662 *
00663 *              Zero
00664 *
00665             IF( ITYPE.EQ.1 ) THEN
00666                IINFO = 0
00667 *
00668             ELSE IF( ITYPE.EQ.2 ) THEN
00669 *
00670 *              Identity
00671 *
00672                DO 70 JCOL = 1, N
00673                   A( JCOL, JCOL ) = ANORM
00674    70          CONTINUE
00675 *
00676             ELSE IF( ITYPE.EQ.3 ) THEN
00677 *
00678 *              Jordan Block
00679 *
00680                DO 80 JCOL = 1, N
00681                   A( JCOL, JCOL ) = ANORM
00682                   IF( JCOL.GT.1 )
00683      $               A( JCOL, JCOL-1 ) = ONE
00684    80          CONTINUE
00685 *
00686             ELSE IF( ITYPE.EQ.4 ) THEN
00687 *
00688 *              Diagonal Matrix, [Eigen]values Specified
00689 *
00690                CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
00691      $                      ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
00692      $                      IINFO )
00693 *
00694             ELSE IF( ITYPE.EQ.5 ) THEN
00695 *
00696 *              Symmetric, eigenvalues specified
00697 *
00698                CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
00699      $                      ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
00700      $                      IINFO )
00701 *
00702             ELSE IF( ITYPE.EQ.6 ) THEN
00703 *
00704 *              General, eigenvalues specified
00705 *
00706                IF( KCONDS( JTYPE ).EQ.1 ) THEN
00707                   CONDS = ONE
00708                ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
00709                   CONDS = RTULPI
00710                ELSE
00711                   CONDS = ZERO
00712                END IF
00713 *
00714                ADUMMA( 1 ) = ' '
00715                CALL SLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE,
00716      $                      ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4,
00717      $                      CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ),
00718      $                      IINFO )
00719 *
00720             ELSE IF( ITYPE.EQ.7 ) THEN
00721 *
00722 *              Diagonal, random eigenvalues
00723 *
00724                CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
00725      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00726      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
00727      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00728 *
00729             ELSE IF( ITYPE.EQ.8 ) THEN
00730 *
00731 *              Symmetric, random eigenvalues
00732 *
00733                CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
00734      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00735      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
00736      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00737 *
00738             ELSE IF( ITYPE.EQ.9 ) THEN
00739 *
00740 *              General, random eigenvalues
00741 *
00742                CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
00743      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00744      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
00745      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00746                IF( N.GE.4 ) THEN
00747                   CALL SLASET( 'Full', 2, N, ZERO, ZERO, A, LDA )
00748                   CALL SLASET( 'Full', N-3, 1, ZERO, ZERO, A( 3, 1 ),
00749      $                         LDA )
00750                   CALL SLASET( 'Full', N-3, 2, ZERO, ZERO, A( 3, N-1 ),
00751      $                         LDA )
00752                   CALL SLASET( 'Full', 1, N, ZERO, ZERO, A( N, 1 ),
00753      $                         LDA )
00754                END IF
00755 *
00756             ELSE IF( ITYPE.EQ.10 ) THEN
00757 *
00758 *              Triangular, random eigenvalues
00759 *
00760                CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
00761      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00762      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
00763      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00764 *
00765             ELSE
00766 *
00767                IINFO = 1
00768             END IF
00769 *
00770             IF( IINFO.NE.0 ) THEN
00771                WRITE( NOUNIT, FMT = 9991 )'Generator', IINFO, N, JTYPE,
00772      $            IOLDSD
00773                INFO = ABS( IINFO )
00774                RETURN
00775             END IF
00776 *
00777    90       CONTINUE
00778 *
00779 *           Test for minimal and generous workspace
00780 *
00781             DO 120 IWK = 1, 2
00782                IF( IWK.EQ.1 ) THEN
00783                   NNWORK = 3*N
00784                ELSE
00785                   NNWORK = MAX( 3*N, 2*N*N )
00786                END IF
00787                NNWORK = MAX( NNWORK, 1 )
00788 *
00789                CALL SGET24( .FALSE., JTYPE, THRESH, IOLDSD, NOUNIT, N,
00790      $                      A, LDA, H, HT, WR, WI, WRT, WIT, WRTMP,
00791      $                      WITMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT,
00792      $                      ISLCT, RESULT, WORK, NNWORK, IWORK, BWORK,
00793      $                      INFO )
00794 *
00795 *              Check for RESULT(j) > THRESH
00796 *
00797                NTEST = 0
00798                NFAIL = 0
00799                DO 100 J = 1, 15
00800                   IF( RESULT( J ).GE.ZERO )
00801      $               NTEST = NTEST + 1
00802                   IF( RESULT( J ).GE.THRESH )
00803      $               NFAIL = NFAIL + 1
00804   100          CONTINUE
00805 *
00806                IF( NFAIL.GT.0 )
00807      $            NTESTF = NTESTF + 1
00808                IF( NTESTF.EQ.1 ) THEN
00809                   WRITE( NOUNIT, FMT = 9999 )PATH
00810                   WRITE( NOUNIT, FMT = 9998 )
00811                   WRITE( NOUNIT, FMT = 9997 )
00812                   WRITE( NOUNIT, FMT = 9996 )
00813                   WRITE( NOUNIT, FMT = 9995 )THRESH
00814                   WRITE( NOUNIT, FMT = 9994 )
00815                   NTESTF = 2
00816                END IF
00817 *
00818                DO 110 J = 1, 15
00819                   IF( RESULT( J ).GE.THRESH ) THEN
00820                      WRITE( NOUNIT, FMT = 9993 )N, IWK, IOLDSD, JTYPE,
00821      $                  J, RESULT( J )
00822                   END IF
00823   110          CONTINUE
00824 *
00825                NERRS = NERRS + NFAIL
00826                NTESTT = NTESTT + NTEST
00827 *
00828   120       CONTINUE
00829   130    CONTINUE
00830   140 CONTINUE
00831 *
00832   150 CONTINUE
00833 *
00834 *     Read in data from file to check accuracy of condition estimation
00835 *     Read input data until N=0
00836 *
00837       JTYPE = 0
00838   160 CONTINUE
00839       READ( NIUNIT, FMT = *, END = 200 )N, NSLCT
00840       IF( N.EQ.0 )
00841      $   GO TO 200
00842       JTYPE = JTYPE + 1
00843       ISEED( 1 ) = JTYPE
00844       IF( NSLCT.GT.0 )
00845      $   READ( NIUNIT, FMT = * )( ISLCT( I ), I = 1, NSLCT )
00846       DO 170 I = 1, N
00847          READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N )
00848   170 CONTINUE
00849       READ( NIUNIT, FMT = * )RCDEIN, RCDVIN
00850 *
00851       CALL SGET24( .TRUE., 22, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT,
00852      $             WR, WI, WRT, WIT, WRTMP, WITMP, VS, LDVS, VS1,
00853      $             RCDEIN, RCDVIN, NSLCT, ISLCT, RESULT, WORK, LWORK,
00854      $             IWORK, BWORK, INFO )
00855 *
00856 *     Check for RESULT(j) > THRESH
00857 *
00858       NTEST = 0
00859       NFAIL = 0
00860       DO 180 J = 1, 17
00861          IF( RESULT( J ).GE.ZERO )
00862      $      NTEST = NTEST + 1
00863          IF( RESULT( J ).GE.THRESH )
00864      $      NFAIL = NFAIL + 1
00865   180 CONTINUE
00866 *
00867       IF( NFAIL.GT.0 )
00868      $   NTESTF = NTESTF + 1
00869       IF( NTESTF.EQ.1 ) THEN
00870          WRITE( NOUNIT, FMT = 9999 )PATH
00871          WRITE( NOUNIT, FMT = 9998 )
00872          WRITE( NOUNIT, FMT = 9997 )
00873          WRITE( NOUNIT, FMT = 9996 )
00874          WRITE( NOUNIT, FMT = 9995 )THRESH
00875          WRITE( NOUNIT, FMT = 9994 )
00876          NTESTF = 2
00877       END IF
00878       DO 190 J = 1, 17
00879          IF( RESULT( J ).GE.THRESH ) THEN
00880             WRITE( NOUNIT, FMT = 9992 )N, JTYPE, J, RESULT( J )
00881          END IF
00882   190 CONTINUE
00883 *
00884       NERRS = NERRS + NFAIL
00885       NTESTT = NTESTT + NTEST
00886       GO TO 160
00887   200 CONTINUE
00888 *
00889 *     Summary
00890 *
00891       CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT )
00892 *
00893  9999 FORMAT( / 1X, A3, ' -- Real Schur Form Decomposition Expert ',
00894      $      'Driver', / ' Matrix types (see SDRVSX for details):' )
00895 *
00896  9998 FORMAT( / ' Special Matrices:', / '  1=Zero matrix.             ',
00897      $      '           ', '  5=Diagonal: geometr. spaced entries.',
00898      $      / '  2=Identity matrix.                    ', '  6=Diagona',
00899      $      'l: clustered entries.', / '  3=Transposed Jordan block.  ',
00900      $      '          ', '  7=Diagonal: large, evenly spaced.', / '  ',
00901      $      '4=Diagonal: evenly spaced entries.    ', '  8=Diagonal: s',
00902      $      'mall, evenly spaced.' )
00903  9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / '  9=Well-cond., ev',
00904      $      'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
00905      $      'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
00906      $      ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
00907      $      'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
00908      $      'lex ', / ' 12=Well-cond., random complex ', '         ',
00909      $      ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
00910      $      'tioned, evenly spaced.     ', ' 18=Ill-cond., small rand.',
00911      $      ' complx ' )
00912  9996 FORMAT( ' 19=Matrix with random O(1) entries.    ', ' 21=Matrix ',
00913      $      'with small random entries.', / ' 20=Matrix with large ran',
00914      $      'dom entries.   ', / )
00915  9995 FORMAT( ' Tests performed with test threshold =', F8.2,
00916      $      / ' ( A denotes A on input and T denotes A on output)',
00917      $      / / ' 1 = 0 if T in Schur form (no sort), ',
00918      $      '  1/ulp otherwise', /
00919      $      ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
00920      $      / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
00921      $      ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
00922      $      '  1/ulp otherwise', /
00923      $      ' 5 = 0 if T same no matter if VS computed (no sort),',
00924      $      '  1/ulp otherwise', /
00925      $      ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
00926      $      ',  1/ulp otherwise' )
00927  9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', '  1/ulp otherwise',
00928      $      / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
00929      $      / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
00930      $      / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
00931      $      '  1/ulp otherwise', /
00932      $      ' 11 = 0 if T same no matter what else computed (sort),',
00933      $      '  1/ulp otherwise', /
00934      $      ' 12 = 0 if WR, WI same no matter what else computed ',
00935      $      '(sort), 1/ulp otherwise', /
00936      $      ' 13 = 0 if sorting succesful, 1/ulp otherwise',
00937      $      / ' 14 = 0 if RCONDE same no matter what else computed,',
00938      $      ' 1/ulp otherwise', /
00939      $      ' 15 = 0 if RCONDv same no matter what else computed,',
00940      $      ' 1/ulp otherwise', /
00941      $      ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
00942      $      / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
00943  9993 FORMAT( ' N=', I5, ', IWK=', I2, ', seed=', 4( I4, ',' ),
00944      $      ' type ', I2, ', test(', I2, ')=', G10.3 )
00945  9992 FORMAT( ' N=', I5, ', input example =', I3, ',  test(', I2, ')=',
00946      $      G10.3 )
00947  9991 FORMAT( ' SDRVSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00948      $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
00949 *
00950       RETURN
00951 *
00952 *     End of SDRVSX
00953 *
00954       END
 All Files Functions