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